数学分析实验

偏导数与全微分 格式 意义 D[f,x] 计算f 关于x的偏导数 D[f,x1,x2,…] 计算f 关于x1,x2,…的高阶偏导数 D[f,{xi,n}] 计算f 关于xi的n阶偏导数 D[f,x,NonConstants->{y}] 计算f 关于x的偏导数，y是x的函数 Dt[f,x] 计算f 关于x的全导数 Dt[f,x1,x2,…] 计算f 关于x1,x2,…的全导数 Dt[f,x,Constants-> {c1,c2,…}] 计算f 关于x的全导数，其中c1,c2,…为常数 Dt[f[x,y]] 求f 的全微分

D[Exp[x+y+z^2],z] D[x^2+y^2,x,y] D[Exp[x+2y],{y,5}] D[x[1]^2+x[2]^2,x[1]] D[g[x^2,y^2],x] D[g[x,y^2],x,x,y] D[x^2+y^2,x,NonConstants->{y}] D[x^2+y[x]^2,x] D[f[Sin[x],y[x^2]],x] Dt[x^2+y^2,x] Dt[x^2+y^2+z^2,x] Dt[x^2+y^2+z^2,x,Constants->{y}] Dt[ArcSin[x/y]]

t1={{"D[Exp[x+y+z^2],z]", "D[x^2+y^2,x,y]", "D[Exp[x+2y],{y,5}]}", "D[x[1]^2+x[2]^2,x[1]]", "D[g[x^2,y^2],x]", "D[g[x,y^2],x,x,y]", "D[x^2+y^2,x,NonConstants->{y}]", "D[x^2+y[x]^2,x]", "D[f[Sin[x],y[x^2]],x]", "Dt[x^2+y^2,x]", "Dt[x^2+y^2+z^2,x]", "Dt[x^2+y^2+z^2,x,Constants->{y}]", "Dt[ArcSin[x/y]]"}};

t2={{D[Exp[x+y+z^2],z], D[x^2+y^2,x,y], D[Exp[x+2y],{y,5}], D[x[1]^2+x[2]^2,x[1]], D[g[x^2,y^2],x], D[g[x,y^2],x,x,y], D[x^2+y^2,x,NonConstants->{y}], D[x^2+y[x]^2,x], D[f[Sin[x],y[x^2]],x], Dt[x^2+y^2,x], Dt[x^2+y^2+z^2,x], Dt[x^2+y^2+z^2,x,Constants->{y}], Dt[ArcSin[x/y]]}}; TableForm[Transpose[Join[t1,t2]]]

隐函数求导 例 已知x2+y2+z2-4z=0, 求 . eq=(x^2+y^2+z^2-4z==0); y/:Dt[y,x]=0; deq1=Dt[eq,x] deq2=Dt[eq,x,x] Solve[{deq1,deq2},{Dt[z,x],Dt[z,{x,2}]}]

例 设x u-y v=0, y u+x v=1,求 eq={x u-y v==0,y u+x v==1}; y/:Dt[y,x]=0; x/:Dt[x,y]=0; eqx=Dt[eq,x] {jie1}=Solve[eqx,{Dt[u,x],Dt[v,x]}] eqy=Dt[eq,y] {jie2}=Solve[eqy,{Dt[u,y],Dt[v,y]}] uv=Solve[eq,{u,v}]; u=u/.uv[[1]];v=v/.uv[[1]]; {D[u,x],D[u,y],D[v,x],D[v,y]}//Simplify {Dt[u,x],Dt[u,y],Dt[v,x],Dt[v,y]}/.jie1/.jie2//Simplify

多元函数的极值 例1 设z=x4-8xy+2y2-3, 求函数的极值点和极值. z=x^4-8x y+2y^2-3; dzx=D[z,x];dzy=D[z,y]; s0=Solve[{dzx==0,dzy==0},{x,y}] dzxx=D[z,x,x];dzxy=D[z,x,y];dzyy=D[z,y,y]; L=dzxx dzyy-dzxy^2 Module[{a,b,c,d,e},Do[a=L/.s0[[k]] ;b= dzxx/.s0[[k]];c=z/.s0[[k]];d=x/.s0[[k]];e=y/.s0[[k]]; If[a>0,If[b<0,Print["(",d,", ",e,") is a maximum point. "," z=",c], Print["(",d,", ",e,") is a minimum point."," z=",c]]], {k,1,Length[s0]}]]

也可以通过图形来观察极值点与驻点 Plot3D[z,{x,-4,4},{y,-6,6},PlotPoints->100,Mesh->False] 缩小值域再作观察: small=Plot3D[z,{x,-4,4},{y,-6,6},PlotPoints->100,PlotRange-> {-20,3},Mesh->False] Show[small,ClipFill->None] ContourPlot[z,{x,-4,4},{y,-6,6},PlotRange->{-30,3}]

条件极值与Lagrange乘数法 例2 求表面积为a2且体积最大的长方体的体积. Clear[x,y,z,r,g] g[x_,y_,z_]:=x y z-k(2x y+2y z+2 x z-a^2) eq1=D[g[x,y,z],x]==0; eq2=D[g[x,y,z],y]==0; eq3=D[g[x,y,z],z]==0; eq4=D[g[x,y,z],k]==0; r=Solve[{eq1,eq2,eq3,eq4},{x,y,z,k}] g[x,y,z]/.r[[2]] Clear[g,eq1,eq2,eq3,eq4,r]

例3 已知两条直线的方程分别为： 求这两条直线间的最短距离. Clear[x,y,z,s,t,r,k1,k2,m1,m2,L] L=(x-s)^2+(y-t)^2+(z-r)^2+k1(x+y-z-1)+m1(2x+y-z-2)+ k2(s+2t-r-2)+m2(s+2t+2r+4); rr=Solve[{D[L,x]==0,D[L,y]==0,D[L,z]==0,D[L,s]==0, D[L,t]==0,D[L,r]==0,D[L,k1]==0,D[L,k2]==0, D[L,m1]==0,D[L,m2]==0},{x,y,z,s,t,r,k1,k2,m1,m2}] d=Sqrt[(x-s)^2+(y-t)^2+(z-r)^2]/.rr[[1]] Print[p1={x,y,z}/.rr[[1]],",",p2={s,t,r}/.rr[[1]],",",p2-p1]

二元函数在区域D内的最大值和最小值 例4 求函数z=x2+y2-4x-2y+7在上半圆域x2+y216, y0上的最大值和最小值. (1)画出函数的等高线和区域D的边界线. Clear[x,y,z,r,g] g[x_,y_]:=x^2+y^2-4x-2y+7 pts=Table[N[{4Cos[t],4Sin[t]}],{t,0,Pi,Pi/20}]; ContourPlot[g[x,y],{x,-5,5},{y,-5,5},Epilog->{RGBColor[1,1,0], Thickness[0.003],Line[pts],Line[{{-4,0},{4,0}}]}]; (2)求函数在区域D内的驻点，计算相应的函数值. r=Solve[{D[g[x,y],x]==0,D[g[x,y],y]==0},{x,y}] g[x,y]/.r[[1]]

(3)求函数在直线边界y=0, -4x4上的最值. Plot[g[x,0],{x,-4,4},PlotRange->All]; r=Solve[D[g[x,0],x]==0,x] Print[g[x,0]/.r[[1]], " , ",g[-4,0], " , ",g[4,0]] (4)求函数在圆弧边界x2+y2=16, y0上的最值. Plot[g[4Cos[t],4Sin[t]],{t,0,Pi}]; eqt=D[g[4Cos[t],4Sin[t]],t]==0; rt=FindRoot[eqt,{t,0.5}] {4Cos[t],4Sin[t]}/.rt g[4Cos[t],4Sin[t]]/.rt 根据上述计算可知：函数在区域D上的点(2,1)处取得最小值2；在边界点(-4,0)处取得最大值39.

平面上的几何变换 函数图形的平移变换 例1 设 ，在同一坐标系中作出y=f(x), y=f(x+1)和y=f(x-2)的图象. f[x_]:=Exp[-x^2/2] Plot[{f[x],f[x+1],f[x-2]},{x,-5,5},PlotStyle->{RGBColor[1,0,0], RGBColor[0,1,0], RGBColor[0,0,1]}] ParametricPlot[{{Cos[t],Sin[t]}, {Cos[t]+3,Sin[t]-2}},{t,0,2Pi}, PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]},AspectRatio ->Automatic] Plot[{f[x],f[x]+2,f[x-2]-3},{x,-5,5},PlotStyle->{RGBColor[1,0,0], RGBColor[0,1,0], RGBColor[0,0,1]}]

函数图形的伸缩变换 例2 设 ，在同一坐标系中作出y=f(x), y=f(2x)和y=f(x/2)的图象. f[x_]:=Exp[-x^2/2] Plot[{f[x],f[2x],f[x/2]},{x,-5,5},PlotStyle->{RGBColor[1,0,0], RGBColor[0,1,0], RGBColor[0,0,1]}] 通过动画观察图象的变化情况： Do[Plot[f[s x],{x,-5,5},PlotStyle->RGBColor[0,0,1], PlotRange->{0,1},PlotLabel->"s="<>ToString[s]],{s,0.2,5,0.2}] 进一步观察t f(s x),f(-x),-f(x)及参数方程图象的变化情况.

函数图形的旋转变换 例3 将抛物线 绕原点旋转30度. z[t_]:=(t+t^2 I)Exp[Pi/6 I] ParametricPlot[{Re[z[t]],Im[z[t]]},{t,-2,2},AspectRatio-> Automatic,PlotRange->{{-5,5},{-1,5}},PlotStyle-> RGBColor[0,0,1]] 动画观察图象旋转的效果: z1[t_,s_]:= (t+t^2 I)Exp[s I] Do[ParametricPlot[{Re[z1[t,s]],Im[z1[t,s]]},{t,-2,2}, AspectRatio-> Automatic,PlotRange->{{-5,5},{-5,5}}, PlotStyle-> RGBColor[0,0,1]],{s,0,2Pi,Pi/12}]

函数图形的对称变换 例4 已知 ，试扩展f(x)的定义域，使f(x)分别成为奇函数和偶函数，作出延拓后的函数图象. (1) 偶延拓 f=.;f[x_]:=Sqrt[2x-x^2]/;x>=0 f[x_]:=f[-x]/;x<0 Plot[f[x],{x,-2,2},AspectRatio->Automatic] (2) 奇延拓 f1[x_]:=Sqrt[2x-x^2]/;x>=0 f1[x_]:=-f1[-x]/;x<0 Plot[f1[x],{x,-2,2},AspectRatio->Automatic]

例5 将例4中的函数f(x)延拓为以2 为周期的周期函数，作出延拓后的函数图象. 周期延拓 Clear[f]; f[x_]:=Which[0<=x<=2,Sqrt[2x-x^2],x>2,f[x-2],x<0,f[x+2]] Plot[f[x],{x,-6,6},AspectRatio->Automatic]

例6 作出例4中的函数y=f(x)关于直线y=ax的镜像. ang[x_]:=ArcTan[x];a=1/2; zf[t_]:=(t+Sqrt[2t-t^2] I); zf1[t_,u_]:=Conjugate[zf[t] Exp[-ang[u]I]]Exp[ang[u]I]; ParametricPlot[{{Re[zf1[t,a]],Im[zf1[t,a]]},{t,t/2},{Re[zf[t]], Im[zf[t]]}},{t,0,2},AspectRatio->Automatic,PlotStyle-> {RGBColor[1,0,0], RGBColor[0,1,0], RGBColor[0,0,1]}] 动画观察图象镜像的效果: Do[ParametricPlot[{{Re[zf1[t,b]],Im[zf1[t,b]]},{t,b t}, {Re[zf[t]],Im[zf[t]]}},{t,0,2},PlotRange->{{-1,2.01},{-1.01,2}}, AspectRatio->Automatic, PlotStyle->{RGBColor[1,0,0], RGBColor[0,1,0], RGBColor[0,0,1]}],{b,0,2,1/12}]

例7 作出例4中的函数y=f(x)关于直线y=3x+5的镜像. a=3;b=5I; zf2[t_,u_]:=Conjugate[(zf[t]-b)Exp[-ang[u]I]]Exp[ang[u]I]+b; ParametricPlot[{{Re[zf2[t,a]],Im[zf2[t,a]]},{t-2,3(t-2)+5}, {Re[zf[t]],Im[zf[t]]}},{t,0,2}, AspectRatio->Automatic, PlotStyle->{RGBColor[1,0,0], RGBColor[0,1,0], RGBColor[0,0,1]}]

几个练习 例1 P143，练习五 z=.;x=.;y=.;point={x->1,y->1,z->1}; eqns={2x-3y+5z-4==0,x^2+y^2+z^2-3x==0}; jie=Solve[Dt[eqns,x],{Dt[y,x],Dt[z,x]}] dzx=Dt[z,x]/.jie[[1]]/.point dyx=Dt[y,x]/.jie[[1]]/.point v0={1,dyx,dzx} (x-1)/1==(y-1)/dyx==(z-1)/dzx Expand[(x-1)+(y-1)dyx+(z-1)dzx]==0

例2 P143，练习六 Module[{a,b,b0,m,u,x,y,z}, m[v_]:=Sqrt[v.v]; u=x^2+y^2+z^2;z[t_]:={t,t^2,t^2}; b=D[z[t],t]/.{t->1};b0=b/m[b]; grad={D[u,x],D[u,y],D[u,z]}/.{x->1,y->1,z->1} ; db0=grad.b0; Print[db0," , ",grad]]

例3 P143，练习七 (2) fff[x_,y_]:=Evaluate[Normal[Series[x^y,{x,1,5},{y,1,5}]]]; N[fff[1101/1000,1021/1000],10]

例4 P163，练习四 Module[{a,s,data,xydata,T1,T2,TT,g1,i}, a=Pi/30;s=Sin[a]+Cos[a]; data={{-1,-1,1},{1,-1,1},{1,1,1},{-1,1,1}}; T1={{Cos[a],Sin[a],0},{-Sin[a],Cos[a],0},{0,0,1}}; T2={{s,0,0},{0,s,0},{0,0,1}};TT=N[T1.T2]; g1={}; For[i=0,i<21,i++, xydata=Transpose[Drop[Transpose[data],-1]]; xydata=Append[xydata,xydata[[1]]]; g1=Append[g1,Graphics[Line[xydata]]]; Show[g1,AspectRatio->1,PlotRange->{{-10,10},{-10,10}}]; data=data.TT]]

复数法 Module[{a,s,data,xydata,g1,i}, a=Pi/30;s=Sin[a]+Cos[a]; data={-1-I,1-I,1+I,-1+I,-1-I}; g1={}; For[i=0,i<21,i++, xydata=Transpose[{Re[data],Im[data]}]; g1=Append[g1, Graphics[{RGBColor[1,0,1],Line[xydata]}]]; Show[g1,AspectRatio->1,PlotRange->{{-10,10},{-10,10}}]; data=N[data Exp[a I] s]] ]

机翼轮廓线的加工 问题的提出 在飞机制造业中，机翼的加工是一项关键技术。由于机翼的尺寸很大，通常在图纸中只能标出某些关键点的数据。下表给出的是某型号飞机的机翼上缘轮廓线的部分数据。 x 0. 4.74 9.50 19.00 38.00 57.00 76.00 y 5.32 8.10 11.97 16.15 17.10 16.34 95.00 114.0 133.0 152.0 171.0 190.0 14.63 12.16 9.69 7.03 3.99 0.0

t={{0.,0.},{4.74,5.32},{9.50,8.10},{19.00,11.97},{38.00,16.15}, {57.00,17.10},{76.00,16.34},{95.00,14.63},{114.00,12.16}, {133.00,9.69},{152.00,7.03},{171.00,3.99},{190.00,0.}}; TableForm[Join[{{x,y}},t]] ListPlot[t,PlotJoined->True,AspectRatio->1/4,PlotStyle-> RGBColor[1,0,0]]

多项式插值 设多项式为 ，根据已知条件，对应每个节点 ，应有 所以系数a0,a1,…,an应满足如下方程组： 当x0,x1,…,xn互不相同时，由n+1个节点可确定唯一的n次插值多项式Pn(x).

例1（P.172) 略. s={{0,853.},{30,921.},{60,934.}}; n=3;eqns={};vars={}; Pn[x_]:=a[0]+Sum[a[k]x^k,{k,1,n-1}] For[i=1,i<=n,i++,eqn[i]=(s[[i,2]]==Pn[s[[i,1]]])] For[i=1,i<=n,i++,AppendTo[eqns,eqn[i]]; AppendTo[vars,a[i-1]]] xishu=Solve[eqns,vars] fy=Pn[x]/.xishu[[1]] Plot[fy,{x,0,60}] FindMinimum[-fy,{x,30}] 待定系数法 通过计算可以看出，当x=52，即6月22日的白天最长，约为935.11分钟。

线性插值 考虑过两个节点(x0,y0),(x1,y1)的线性插值. 抛物线插值 考虑过节点(x0,y0),(x1,y1), (x2,y2)的抛物线插值. x0=0;x1=30;x2=60;y0=853.;y1=921.;y2=934.; l[x_,a_,b_,c_]:=(x-b)/(a-b) (x-c)/(a-c) L0[x_]:=l[x,x0,x1,x2];L1[x_]:=l[x,x1,x0,x2]; L2[x_]:=l[x,x2,x0,x1] p12[x_]:=L0[x] y0+L1[x]y1+L2[x]y2 FindMinimum[-p12[x],{x,30}] 基函数法

Lagrange插值 InterpolatingPolynomial[data,x] InterpolatingPolynomial[{1,4,5,8},x]//Simplify s={{0,853.},{30,921.},{60,934.}}; InterpolatingPolynomial[s,x]//Expand

三次样条插值 三次样条函数 满足以下： (a) 在每一区间段[xi-1,xi]上都是三次多项式，记为Si(x). (b) S(x)在区间[a,b]上有二阶连续导数. (c) 在插值节点处满足插值条件S(xi)=yi, i=0,1,2,…,n. 使用方法 <<NumericalMath`SplineFit` SplineFit[Data,Cubic] sindata=Table[{t,Sin[t]},{t,0,6Pi,0.3Pi}]; splinesin=SplineFit[sindata,Cubic] ParametricPlot[splinesin[t],{t,0,20},Compiled->False] fff=SplineFit[t,Cubic]; ParametricPlot[fff[t],{t,0,12}, AspectRatio->1/3]

3.4 数控机床的刀具补偿 实际问题：刀具补偿的估计 在数控机床加工零件时，由于刀具磨损会影响加工精度，要对刀具的磨损进行补偿. 为测定刀具的磨损速度，实验室每隔一小时测量一次刀具的厚度（单位: mm），得到实验数据如下： 时间 1 2 3 4 5 6 7 厚度 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8

t={0,1,2,3,4,5,6,7}; y={27.0,26.8,26.5,26.3,26.1,25.7,25.3,24.8}; data=Transpose[{t,y}] ph0=ListPlot[data,PlotStyle-> {RGBColor[1,0,0],PointSize[0.02]]

最小二乘法 从数据图可以看出：这些点大致位于一条直线的附近。因此设y与t近似满足线性关系，设y=a+b t. 可以使用最小二乘法来确定a,b. M[a_,b_]:=Sum[(y[[i]]-(a+b t[[i]]))^2,{i,1,8}] eqt={D[M[a,b],a]==0,D[M[a,b],b]==0}//Simplify r=Solve[eqt,{a,b}] a=.;b=.; ss1=a+b t/.r[[1]] Ph1=Plot[ss1,{x,0,7},PlotStyle->RGBColor[0,1,0]]

Mathematica系统提供了Fit专用函数来实现在最小二乘原则下的线性拟合运算. Fit[Data,funcs,var] ss1=Fit[data,{1,x},x] ss2=Fit[data,{1,x,x^2},x] ss3=Fit[data,{1,x,x^2,x^3},x] ss4=Fit[data,{1,Sin[x],x^2},x] ph2=Plot[{ss1,ss2,ss3,ss4},{x,0,7},PlotStyle-> {RGBColor[0,1,0], RGBColor[0,0,1], RGBColor[1,1,0], RGBColor[0,1,1]}] Show[ph0,ph1]

曲线拟合 例 在研究某单质分子的化学反应速度时，得到如下数据： t 3 6 9 12 15 18 21 24 y 57.6 41.9 31.0 22.7 16.6 12.2 8.9 6.5 试根据实验数据确定经验分布函数y=f(x). t1={3,6,9,12,15,18,21,24}; y1={57.6,41.9,31.0,22.7,16.6,12.2,8.9,6.5}; data1=Transpose[{t1,y1}]; d2=ListPlot[data1,PlotStyle-> {RGBColor[0,0,1],PointSize[.015]}]

根据化学反应速度理论，y=f(t)应该是指数函数y=keat，其中a, k是待定常数。取对数得，lny=lnk+at=b+at，即lny 与t满足线性关系。 logy=Log[y1]; data2=Transpose[{t1,logy}]; d3=ListPlot[data2, PlotStyle-> {RGBColor[0,0,1],PointSize[.015]}] ly=Fit[data2,{1,x},x] y=Exp[ly]//Factor g=Plot[y,{x,1,25}, PlotStyle->RGBColor[1,0,0]] Show[d2,g] sss=Fit[data1,{1,x,x^2,x^3},x];g1=Plot[sss,{x,1,25}]; Show[d2,g1]