数学实验 Mathematic实验三 导数
高数项目-微积分(Mathematica)
附录Ⅰ 大学数学实验指导书（仅供交流学习）项目三 多元函数微积分实验1 多元函数微分学（基础实验）实验目的 掌握利用Mathematica 计算多元函数偏导数和全微分的方法, 掌握计算二元 函数极值和条件极值的方法. 理解和掌握曲面的切平面的作法. 通过作图和观察, 理解二元 函数的性质、方向导数、梯度和等高线的概念.基本命令1.求偏导数的命令D命令D 既可以用于求一元函数的导数, 也可以用于求多元函数的偏导数. 例如: 求),,(z y x f 对x 的偏导数, 则输入D[f[x,y,z],x] 求),,(z y x f 对y 的偏导数, 则输入D[f[x,y,z],y]求),,(z y x f 对x 的二阶偏导数, 则输入D[f[x,y,z],{x,2}] 求),,(z y x f 对y x ,的混合偏导数, 则输入D[f[x,y,z],x,y] …………2.求全微分的命令Dt该命令只用于求二元函数),(y x f 的全微分时, 其基本格式为Dt[f[x,y]]其输出的表达式中含有Dt[x],Dt[y], 它们分别表示自变量的微分d x ,d y . 若函数),(y x f 的表 达式中还含有其它用字符表示的常数, 例如a, 则Dt[f[x,y]]的输出中还会有Dt[a], 若采用选 项Constants->{a}, 就可以得到正确结果, 即只要输入Dt[f[x,y],Constants->{a}]3.在Oxy 平面上作二元函数),(y x f 的等高线的命令ContourPlot 命令的基本格式为ContourPlot[f[x,y],{x,x1,x2},{y,y1,y2}]例如,输入ContourPlot[x^2-y^2,{x,-2,2},{y,-2,2}]则输出函数22y x z -=的等高线图(图1.1). 该命令的选项比较多(详细的内容参见光盘中的实验案例库). 如选项Contours->15表示作15条等高线, 选项Contours->{0}表示只作函数值为0的等高线.实验举例求多元函数的偏导数与全微分例1.1 (教材 例1.1) 设),(cos )sin(2xy xy z +=求.,,,222yx z x z y z x z ∂∂∂∂∂∂∂∂∂ 输入Clear[z];z=Sin[x*y]+Cos[x*y]^2; D[z,x] D[z,y] D[z,{x,2}] D[z,x,y]则输出所求结果.例1.2 设,)1(y xy z +=求yzx z ∂∂∂∂,和全微分dz. 输入Clear[z];z=(1+x*y)^y; D[z,x] D[z,y]则有输出⎪⎪⎭⎫ ⎝⎛++++++-]1[1)1()1(12xy Log xy xy xy xy y y y再输入Dt[z]则得到输出⎪⎪⎭⎫⎝⎛+++++]1[][1])[][()1(xy Log y Dt xy y xDt x yDt y xy y 例1.3 (教材 例1.2) 设,)(y xy a z +=其中a 是常数, 求dz.输入Clear[z,a];z=(a+x*y)^y;wf=Dt[z,Constants->{a}]//Simplify则输出结果:(a+xy)-1+y (y 2Dt[x,Constants->{a}]+Dt[y,Constants->{a}](xy+(a+xy)Log[a+xy]))其中Dt[x,Constants->{a}]就是d x , Dt[y,Constants->{a}]就是d y . 可以用代换命令“/.”把它们 换掉. 输入wf/.{Dt[x,Constants->{a}]->dx,Dt[y,Constants->{a}]->dy}输出为(a+xy)-1+y (dxy 2+dy(xy+(a+xy)Log[a+xy]))例1.4 (教材 例1.3) 设v u e y v u e x u u cos ,sin -=+=,求.,,,yv x v y u x u ∂∂∂∂∂∂∂∂ 输入eq1=D[x==E^u+u*Sin[v],x,NonConstants->{u,v}](*第一个方程两边对x 求导数, 把u,v 看成x,y 的函数*) eq2=D[y==E^u-u*Cos[v],x,NonConstants->{u,v}](*第二个方程两边对x 求导数, 把u,v 看成x,y 的函数*) Solve[{eq1,eq2},{D[u,x,NonConstants->{u,v}],D[v,x,NonConstants->{u,v}]}]//Simplify(*解求导以后由eq1,eq2组成的方程组*)则输出}}][][1(][}],{tan ,,[,][][1][}],{tan ,,[{{v Sin E v Cos E u v Cos E v u ts NonCons x v D v Sin E v Cos E v Sin v u ts NonCons x u D u u u u u -+-->->-+->->-其中D[u,x,NonConstants->{u,v}]表示u 对x 的偏导数, 而D[v,x,NonCosnstants->{u,v}]表示v 对x 的偏导数. 类似地可求得u ,v 对y 的偏导数.微分学的几何应用例1.5 求出曲面222y x z +=在点(1,1)处的切平面、法线方程, 并画出图形. 解(1) 画出曲面的图形. 曲面的参数方程为⎪⎩⎪⎨⎧=∈∈==2]2,0[],2,0[,cos 2/sin rz r u u r y u f x π 输入命令Clear[f];f[x_,y_]=2x^2+y^2;p1=Plot3D[f[x,y],{x,-2,2},{y,-2,2}];g1=ParametricPlot3D[{r*Sin[u]/Sqrt[2.],r*Cos[u],r^2}, {u,0,2*Pi},{r,0,2}] 则输出相应图形(图1.2).(2) 画出切平面的图形. 输入命令a=D[f[x,y],x]/.{x->1,y->1}; b=D[f[x,y],y]/.{x->1,y->1}; p[x_,y_]=f[1,1]+a(x-1)+b(y-1);g2=Plot3D[p[x,y],{x,-2,2},{y,-2,2}];则输出切平面方程为,012=-+y x 及相应图形(图1.3).(3) 画出法线的图形. 输入命令ly[x_]=1+b(x-1)/a;lz[x_]=f[1,1]-(x-1)/a;g3=ParametricPlot3D[{x,ly[x],lz[x]},{x,-2,2}]; Show[p1,g2,g3,AspectRatio->Automatic,ViewPoint->{-2.530,-1.025,2.000}];则输出相应图形(图1.4).图1.4例1.6 (教材 例1.4) 求曲面14),(22++=y x y x k 在点⎪⎭⎫⎝⎛2164,21,41处的切平面方程, 并把曲面和它的切平面作在同一图形里.输入Clear[k,z];k[x_,y_]=4/(x^2+y^2+1); (*定义函数k(x,y)*)kx=D[k[x,y],x]/.{x->1/4,y->1/2};(*求函数k(x,y)对x 的偏导数, 并代入在指定点的值*) ky=D[k[x,y],y]/.{x->1/4,y->1/2};(*求函数k(x,y)对y 的偏导数, 并代入在指定的值*) z=kx*(x-1/4)+ky*(y-1/2)+k[1/4,1/2]; (*定义在指定点的切平面函数*)再输入qm=Plot3D[k[x,y],{x,-2,2},{y,-2,2},PlotRange->{0,4}, BoxRatios->{1,1,1},PlotPoints->30, DisplayFunction->Identity]; qpm=Plot3D[z,{x,-2,2},{y,-2,2}, DisplayFunction->Identity];Show[qm,qpm,DisplayFunction->$DisplayFunction]多元函数的极值例1.7 (教材 例1.5) 求x y x y x y x f 933),(2233-++-=的极值. 输入Clear[f];f[x_,y_]=x^3-y^3+3x^2+3y^2-9x; fx=D[f[x,y],x] fy=D[f[x,y],y]critpts=Solve[{fx==0,fy==0}]则分别输出所求偏导数和驻点:2236369y y x x -++-{{x->-3,y->0},{x->-3,y->2},{x->1,y->0},{x->1,y->2}}再输入求二阶偏导数和定义判别式的命令fxx=D[f[x,y],{x,2}]; fyy=D[f[x,y],{y,2}]; fxy=D[f[x,y],x,y]; disc=fxx*fyy-fxy^2输出为判别式函数2xy yy xx f f f -的形式:(6+6x)(6-6y)再输入data={x,y,fxx,disc,f[x,y]}/.critpts;TableForm[data,TableHeadings->{None,{ "x ", "y ", "fxx ", "disc ", "f "}}]最后我们得到了四个驻点处的判别式与xx f 的值并以表格形式列出.X y fxx disc f -3 0 -12 -72 27 -3 2 -12 72 31 1 0 12 72 -5 1 2 12 -72 -1易见,当2,3=-=y x 时,12-=xx f 判别式disc=72, 函数有极大值31;当0,1==y x 时,12=xx f 判别式disc=72, 函数有极小值-5;当0,3=-=y x 和2,1==y x 时, 判别式disc=-72, 函数在这些点没有极值. 最后,把函数的等高线和四个极值点用图形表示出来,输入d2={x,y}/.critpts;g4=ListPlot[d2,PlotStyle->PointSize[0.02],DisplayFunction->Identity]; g5=ContourPlot[f[x,y],{x,-5,3},{y,-3,5},Contours->40,PlotPoints->60,ContourShading->False,Frame->False,Axes->Automatic,AxesOrigin->{0,0},DisplayFunction->Identity];Show[g4,g5,DisplayFunction->$DisplayFunction]则输出图1.6.从上图可见, 在两个极值点附近, 函数的等高线为封闭的. 在非极值点附近, 等高线不 封闭. 这也是从图形上判断极值点的方法.注:在项目一的实验4中,我们曾用命令FindMinimum 来求一元函数的极值, 实际上,也可 以用它求多元函数的极值, 不过输入的初值要在极值点的附近. 对本例,可以输入以下命令FindMinimum[f[x,y],{x,-1},{y,1}]则输出{-5.,{x->1.,y->-2.36603×10-8}}从中看到在0,1==y x 的附近函数),(y x f 有极小值-5, 但y 的精度不够好.例1.8 求函数22y x z +=在条件0122=-+++y x y x 下的极值. 输入Clear[f,g,la]; f[x_,y_]=x^2+y^2;g[x_,y_]=x^2+y^2+x+y-1; la[x_,y_,r_]=f[x,y]+r*g[x,y]; extpts=Solve[{D[la[x,y,r],x]==0,D[la[x,y,r],y]==0,D[la[x,y,r],r]==0}]得到输出⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧+->-+->-+->-⎩⎨⎧⎭⎬⎫⎩⎨⎧-->--->--->-)31(21),31(21),33(31,)31(21),31(21),33(31y x r y x r再输入f[x,y]/.extpts//Simplify得到两个可能是条件极值的函数值}.32,32{-+但是否真的取到条件极值呢? 可利用等高线作图来判断.输入dian={x,y}/.Table[extpts[[s,j]],{s,1,2},{j,2,3}] g1=ListPlot[dian,PlotStyle->PointSize[0.03],DisplayFunction->Identity]cp1=ContourPlot[f[x,y],{x,-2,2},{y,-2,2},Contours->20,PlotPoints->60,ContourShading->False,Frame->False,Axes-> Automatic,AxesOrigin->{0,0},DisplayFunction->Identity]; cp2=ContourPlot[g[x,y],{x,-2,2},{y,-2,2},PlotPoints->60,Contours->{0},ContourShading-> False,Frame->False,Axes->Automatic,ContourStyle->Dashing[{0.01}],AxesOrigin->{0,0},DisplayFunction->Identity]; Show[g1,cp1,cp2,AspectRatio->1,DisplayFunction->$DisplayFunction]输出为⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧+-+-⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧----)31(21,2321,)31(21,2321及图1.7. 从图可见，在极值可疑点,2321,2321⎪⎪⎭⎫ ⎝⎛----⎪⎪⎭⎫ ⎝⎛+-+-2321,2321 处, 函数),(y x f z =的等高线与曲线0),(=y x g (虚线)相切. 函数),(y x f z =的等高线是一系列同心圆, 由里向外, 函数值在增大, 在)31(21),31(21--=--=y x 的附近观察, 可以得出),(y x f z =取条件极大的结论. 在),31(21+-=x)31(21+-=y 的附近观察, 可以得出),(y x f z =取条件极小的结论.梯度场例1.9 画出函数222),,(y x z z y x f --=的梯度向量. 解 输入命令<<Graphics`ContourPlot3D` <<Graphics`PlotField3D` <<Calculus`VectorAnalysis`SetCoordinates[Cartesian[x,y,z]];f=z^2-x^2-y^2;cp3d=ContourPlot3D[f,{x,-1.1,1.1},{y,-1.1,1.1},{z,-2,2},Contours->{1.0},Axes->True,AxesLabel->{"x","y","z"}];vecplot3d=PlotGradientField3D[f,{x,-1.1,1.1},{y,-1.1,1.1},{z,-2,2},PlotPoints->3,VectorHeads->True];Show[vecplot3d, cp3d];则输出相应图形(图1.8)例1.10 在同一坐标面上作出⎪⎪⎭⎫⎝⎛++=2211),(y x x y x u 和 ,11),(22⎪⎪⎭⎫⎝⎛+-=y x y y x v 的等高线图(0>x ), 并给出它们之间的关系.解 输入命令<<Calculus`VectorAnalysis` <<Graphics`PlotField`SetCoordinates[Cartesian[x,y,z]];check[u_,v_]:={Grad[u][[1]]-Grad[v][[2]],Grad[v][[1]]+Grad[u][[2]]} u=x(1+1/(x^2+y^2));v=y(1-1/(x^2+y^2)); check[u,v]//Simplifyugradplot=PlotGradientField[u,{x,-2,2},{y,-2,2},DisplayFunction->Identity];uplot=ContourPlot[u,{x,-2,2},{y,-2,2},ContourStyle->GrayLevel[0],ContourShading->False,DisplayFunction->Identity,Contours->40,PlotPoints->40]; g1=Show[uplot,ugradplot,DisplayFunction->$DisplayFunction];vgradplot=PlotGradientField[v,{x,-2,2},{y,-2,2},DisplayFunction->Identity];vplot=ContourPlot[v,{x,-2,2},{y,-2,2},ContourStyle->GrayLevel[0.7],ContourShading->False,DisplayFunction->Identity,Contours->40,PlotPoints->40]; g2=Show[vplot,vgradplot,DisplayFunction->$DisplayFunction]; g3=Show[uplot,vplot,DisplayFunction->$DisplayFunction];g4=Show[ugradplot,vgradplot,DisplayFunction->$DisplayFunction];则输出相应图形(图1.9)，其中(a) ),(y x u 的梯度与等高线图;(b) ),(y x v 的梯度与等高线图; (c) ),(y x u 与),(y x v 的等高线图; (d) ),(y x u 与),(y x v 的梯度图.图1.9从上述图中可以看出它们的等高线为一族正交曲线. 事实上, 有,,2222x v y x x y u y v y x x x u ∂∂-=+=∂∂∂∂=+=∂∂ 且,0=∇⋅∇v u 它们满足拉普拉斯方程022222222=∂∂+∂∂=∂∂+∂∂y vx v y u x u 例1.11 (教材 例1.6) 设,),()(22y x xe y x f +-=作出),(y x f 的图形和等高线, 再作出它的梯度向量gradf 的图形. 把上述等高线和梯度向量的图形叠加在一起, 观察它们之间的关系.输入调用作向量场图形的软件包命令<<Graphics\PlotField.m再输入Clear[f];f[x_,y_]=x*Exp[-x^2-y^2];dgx=ContourPlot[f[x,y],{x,-2,2},{y,-2,2},PlotPoints->60, Contours->25,ContourShading->False,Frame->False,Axes->Automatic,AxesOrigin->{0,0}] td=PlotGradientField[f[x,y],{x,-2,2},{y,-2,2},Frame->False] Show[dgx,td]输出为图1.10. 从图可以看到Oxy 平面上过每一点的等高线和梯度向量是垂直的, 且梯度的 方向是指向函数值增大的方向.图1.10例1.12 求出函数244),(y xy x y x f +-=的极值, 并画出函数),(y x f 的等高线、驻点以及),(y x f -的梯度向量的图形.输入命令<<Graphics`PlotField`f=x^4-4*x*y+y^2;FindMinimum[f,{x,1},{y,1}]conplot=ContourPlot[f,{x,-2,2},{y,-3,3},ContourShading->False,PlotPoints->100,Cont ours->{-4,-2,0,2,4,10,20}];fieldplot=PlotGradientField[-f,{x,-2,2},{y,-3,3},ScaleFunction->(Tanh[#/5]&)];critptplot=ListPlot[{{-Sqrt[2],-2*Sqrt[2]},{0,0},{Sqrt[2],2*Sqrt[2]}},PlotStyle->{Point Size[0.03]}];Show[conplot,fieldplot,critptplot];则得到),(y x f 的最小值.4)82843.2,41421.1(-=f 以及函数的图形(图1.11).实验习题 1.设,xy e z =求.dz2.设),,(y xy f z =求.,,22222y x zy z x z ∂∂∂∂∂∂∂3.设),sin (cos ),(228/)(22y x e y x g y x +=+-求.,,2yx z y z x z ∂∂∂∂∂∂∂4.试用例1.5的方法求265433051830120),(xy x x x x y x f +++--=的极值.5.求324y x z +=在01422=-+y x 条件下的极值.6.作出函数42210/)2(),(y x e y x f +-=的等高线和梯度线的图形, 并观察梯度线与等高线的 关系.实验2 多元函数积分学（基础实验）实验目的掌握用Mathematica 计算二重积分与三重积分的方法; 深入理解曲线积分、曲面积分的 概念和计算方法. 提高应用重积分和曲线、曲面积分解决各种问题的能力.基本命令1. 计算重积分的命令lntegrate 和NIntegrate 例如,计算dydx xy x ⎰⎰1002, 输入Integrate[x*y^2,{x,0,1},{y,0,x}]则输出 151又如,计算dydx xy )sin(1102⎰⎰的近似值, 输入NIntegrate[Sin[x*y^2],{x,0,1},{y,0,1}] 则输出 0.160839注: Integrate 命令先对后边的变量积分.计算三重积分时,命令Integrate 的使用格式与计算二重积分时类似. 由此可见, 利用 Mathematica 计算重积分, 关键是确定各个积分变量的积分限. 2. 柱坐标系中作三维图形的命令CylindricalPlot3D使用命令Cylindricalplot3D, 首先要调出作图软件包. 输入 <<Graphics`ParametricPlot3D` 执行成功后便可继续下面的工作.使用命令Cylindricalplot3D 时,一定要把z 表示成r ,θ的函数. 例如,在直角坐标系中方 程22y x z +=是一旋转抛物面, 在柱坐标系中它的方程为2r z =. 因此,输入 CylindricalPlot3D[r^2,{r,0,2},{t,0,2Pi}] 则在柱坐标系中作出了该旋转抛物面的图形.3. 球面坐标系中作三维图形命令SphericalPlot3D使用命令SphericalPlot3D, 首先要调出作图软件包. 输入 <<Graphics`ParametricPlot3D` 执行成功后便可继续下面的工作.命令SphericalPlot3D 的基本格式为SphericalPlot3D[r[],θϕ, {}],,{},,,2121θθθϕϕϕ其中r[],θϕ是曲面的球面坐标方程, 使用时一定要把球面坐标中的r 表示成ϕ、θ的函数. 例如,在球面坐标系中作出球面,22222=++z y x 输入Sphericalplot3D[2,{u,0,pi},|v,0,2,pi|,plotpoints->40]则在球面坐标系中作出了该球面的图形.4. 向量的内积用“.”表示两个向量的内积. 例如,输入 vecl={al,bl,cl} vec2={a2,b2,c2} 则定义了两个三维向量, 再输入 vec1. vec2 则得到它们的内积a1a2+b1b2+c1c2实验举例计算重积分 例2.1 (教材 例2.1) 计算,2dxdy xyD⎰⎰ 其中D 为由,,2y x y x ==+ 2=y 所围成的有界区域.先作出区域D 的草图, 易直接确定积分限,且应先对x 积分, 因此, 输入 Integrate[x*y^2,{y,1,2},{x,2-y,Sqrt[y]}] 则输出所求二重积分的计算结果.120193例2.2 (教材 例2.2) 计算,)(22dxdy e Dy x⎰⎰+- 其中D 为.122≤+y x如果用直角坐标计算, 输入Clear[f,r];f[x,y]=Exp [-(x^2+y^2)];Integrate[f[x,y],{x,-1,1},{y,-Sqrt[1-x^2],Sqrt[1-x^2]}]则输出为dx x 1Erf e 211x 2⎥⎦⎤⎢⎣⎡-π⎰--其中Erf 是误差函数. 显然积分遇到了困难.如果改用极坐标来计算, 也可用手工确定积分限. 输入Integrate[(f[x,y]/.{x->r*Cos[t],y->r*Sin[t]})*r,{t,0,2 Pi},{r,0,1}] 则输出所求二重积分的计算结果eπ-π 如果输入NIntegrate[(f[x,y]/.{x->r*Cos[t],y->r*Sin[t]})*r,{t,0,2 Pi},{r,0,1}] 则输出积分的近似值1.98587例2.3 (教材 例2.3) 计算dxdydz z y x)(22++⎰⎰⎰Ω, 其中Ω由曲面222y x z --=与22y x z +=围成.先作出区域Ω的图形. 输入g1=ParametricPlot3D[{Sqrt[2]*Sin[fi]*Cos[th],Sqrt[2]*Sin[fi]*Sin[th], Sqrt[2]*Cos[fi]},{fi,0,Pi/4},{th,0,2Pi}] g2=ParametricPlot3D[{z*Cos[t],z*Sin[t],z},{z,0,1},{t,0,2Pi}] Show[g1,g2,ViewPoint->{1.3,-2.4,1.0}]则分别输出三个图形（图2.1(a), (b), (c)）.考察上述图形, 可用手工确定积分限. 如果用直角坐标计算, 输入 g[x_,y_,z_]=x^2+y^2+z;Integrate[g[x,y,z],{x,-1,1},{y,-Sqrt[1-x^2], Sqrt[1-x^2]},{z,Sqrt[x^2+y^2],Sqrt[2-x^2-y^2]}] 执行后计算时间很长, 且未得到明确结果.现在改用柱面坐标和球面坐标来计算. 如果用柱坐标计算,输入Integrate[(g[x,y,z]/.{x->r*Cos[s],y->r*Sin[s]})*r,{r,0,1},{s,0,2Pi},{z,r,Sqrt[2-r^2]}]则输出π⎪⎪⎭⎫⎝⎛+-15281252 如果用球面坐标计算,输入Integrate[(g[x,y,z]/.{x->r*Sin[fi]*Cos[t],y->r*Sin[fi]*Sin[t],z->r*Cos[fi]})*r^2*Sin[fi],{s,0,2Pi},{fi,0,Pi/4},{r,0,Sqrt[2]}]则输出π⎪⎪⎭⎫ ⎝⎛+-321662551这与柱面坐标的结果相同.重积分的应用例2.4 求由曲面()y x y x f --=1,与()222,y x y x g --=所围成的空间区域Ω的体积.输入Clear[f,g];f[x_,y_]=1-x -y;g[x_,y_]=2-x^2-y^2;Plot3D[f[x,y],{x,-1,2},{y,-1,2}] Plot3D[g[x,y],{x,-1,2},{y,-1,2}] Show[%,%%]一共输出三个图形, 最后一个图形是图2.1.首先观察到Ω的形状. 为了确定积分限, 要把两曲面的交线投影到Oxy 平面上输入 jx=Solve[f[x,y]==g[x,y],y] 得到输出 ⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧⎪⎭⎫ ⎝⎛-++→⎭⎬⎫⎩⎨⎧⎪⎭⎫ ⎝⎛-+-→22445121,445121x x y x x y为了取出这两条曲线方程, 输入 y1=jx[[1,1,2]] y2=jx[[2,1,2]] 输出为⎪⎭⎫ ⎝⎛-+-2445121x x⎪⎭⎫ ⎝⎛-++2445121x x再输入tu1=Plot[y1,{x,-2,3},PlotStyle->{Dashing[{0.02}]},DisplayFunction->Identity];tu2=Plot[y2,{x,-2,3},DisplayFunction->Identity]; Show[tu1,tu2,AspectRatio->1, DisplayFunction-> $DisplayFunction]输出为图2.2, 由此可见,1y 是下半圆(虚线),2y 是上半圆,因此投影区域是一个圆.设21y y =的解为1x 与2x ,则21,x x 为x 的积分限. 输入 xvals=Solve[y1==y2,x]输出为 ()()⎭⎬⎫⎩⎨⎧⎭⎬⎫⎩⎨⎧+→⎭⎬⎫⎩⎨⎧-→6121,6121x x 为了取出21,x x , 输入x1=xvals[[1,1,2]]x2=xvals[[2,1,2]]输出为()6121- ()6121+ 这时可以作最后的计算了. 输入V olume=Integrate[g[x,y]-f[x,y],{x,x1,x2},{y,y1,y2}]//Simplify 输出结果为 89π例2.5 (教材 例2.4) 求旋转抛物面224y x z --=在Oxy 平面上部的面积.S 先调用软件包, 输入<<Graphics`ParametricPlot3D` 再输入CylindricalPlot3D[4-r^2,{r,0,2},{t,0,2 Pi}] 则输出图2.3.利用计算曲面面积的公式⎰⎰++=xyD y z dxdy z z S 221, 输入Clear[z,z1];z=4-x^2-y^2;z=Sqrt[D[z,x]^2+D[z,y]^2+1]输出为22441y x ++, 因此,利用极坐标计算. 再输入z1=Simplify[z/.{x->r*Cos[t],y->r*Sin[t]}]; Integrate[z1*r,{t,0,2 Pi},{r,0,2}]//Simplify则输出所求曲面的面积()π1717161+-例2.6 在Oxz 平面内有一个半径为2的圆, 它与z 轴在原点O 相切, 求它绕z 轴旋转一周所得旋转体体积.先作出这个旋转体的图形. 因为圆的方程是,422x z x =+它绕z 轴旋转所得的圆环面的方程为)(16)(222222y x z y x +=++,所以圆环面的球坐标方程是.sin 4φ=r 输入SphericalPlot3D[4 Sin[t],{t,0,Pi},{s,0,2 Pi},PlotPoints->30,ViewPoint->{4.0,0.54,2.0}]输出为图2.4.图2.4这是一个环面, 它的体积可以用三重积分计算(用球坐标). 输入 Integrate[r^2*Sin[t],{s,0,2 Pi},{t,0,Pi},{r,0,4 Sin[t]}] 得到这个旋转体的体积为216π计算曲线积分例2.7 (教材 例2.5) 求⎰Lds z y x f ),,(, 其中(),10301,,2y x z y x f ++=积分路径为L :,3,,22t z t y t x ===.20≤≤y注意到,弧长微元dt z y x ds t t t 222++=, 将曲线积分化为定积分,输入 Clear[x,y,z];luj={t,t^2,3t^2}; D[luj,t]则输出z y x ,,对t 的导数 }6,2,1{t t再输入ds=Sqrt[D[luj,t].D[luj,t]];Integrate[(Sqrt[1+30 x^2+10y]/.{x->t, y->t^2,z->3t^2})*ds,{t,0,2}]则输出所求曲线积分的结果:326/3.例2.8 (教材 例2.6) 求dr F L.⎰, 其中.20,sin cos 2)(,)2(356π≤≤+=++=t tj ti t r j xy x i xy F输入vecf={x*y^6,3x*(x*y^5+2)};vecr={2*Cos[t],Sin[t]};Integrate[(vecf.D[vecr,t])/.{x->2Cos[t],y->Sin[t]}, {t,0,2 Pi}]则输出所求积分的结果12π例2.9 求锥面0,222≥=+z z y x 与柱面x y x =+22的交线的长度.先画出锥面和柱面的交线的图形. 输入g1=ParametricPlot3D[{Sin[u]*Cos[v], Sin[u]*Sin[v], Sin[u]}, {u,0,Pi},{v,0,2Pi},DisplayFunction->Identity]; g2=ParametricPlot3D[{Cos[t]^2,Cos[t]*Sin[t],z}, {t,0,2Pi},{z,0,1.2}, DisplayFunction->Identity]; Show[g1,g2,ViewPoint->{1,-1,2},DisplayFunction->$DisplayFunction]输出为图2.5.输入直接作曲线的命令ParametricPlot3D[{Cos[t]^2,Cos[t]*Sin[t],Cos[t]},{t,-Pi/2,Pi/2}, ViewPoint->{1,-1,2},Ticks->False]输出为图2.6.为了用线积分计算曲线的弧长, 必须把曲线用参数方程表示出来. 因为空间曲线的投影曲线的方程为x y x =+22, 它可以化成t x 2cos =,,sin cos t t y =再代入锥面方程222z y x =+, 得[]().2/,2/cos ππ=∈=t t z因为空间曲线的弧长的计算公式是()()()⎰'+'+'=21222t t dt t z t y t x s ,因此输入Clear[x,y,z]; x=Cos[t]^2; y=Cos[t]*Sin[t]; z=Cos[t]; qx={x,y,z};95Integrate[Sqrt[D[qx,t]. D[qx,t]]//Simplify, {t,-Pi/2,Pi/2}]输出为 2Elliptice[-1]这是椭圆积分函数. 换算成近似值. 输入 %//N 输出为3.8202计算曲面积分例2.10 (教材 例2.7) 计算曲面积分⎰⎰∑++dS zx yz xy )(, 其中∑为锥面22y x z +=被柱面x y x 222=+所截得的有限部分.注意到,面积微元dxdy z z dS y x 221++=, 投影曲线x y x 222=+的极坐标方程为,22,cos 2ππ≤≤-=t t r将曲面积分化作二重积分,并采用极坐标计算重积分.输入Clear[f,g,r,t];f[x_,y_,z_]=x*y+y*z+z*x; g[x_,y_]=Sqrt[x^2+y^2];mj=Sqrt[1+D[g[x,y],x]^2+D[g[x,y],y]^2]//Simplify; Integrate[(f[x,y,g[x,y]]*mj/.{x->r*Cos[t],y->r* Sin[t]})*r,{t,-Pi/2,Pi/2},{r,0,2Cos[t]}]则输出所求曲面积分的计算结果15264例2.11 计算曲面积分,333dxdy z dzdx y dydz x ++⎰⎰∑其中∑为球面2222a x y x =++的外侧.可以利用两类曲面积分的关系, 化作对曲面面积的曲面积分⎰⎰∑nds A .. 这里{}{}a z y x n z y x A /,,,,,333==. 因为球坐标的体积元素,sin 2θϕϕd drd r dv =注意到在球面∑上a r =, 取1=dr 后得到面积元素的表示式:().20,0sin 2πθπϕθϕθ≤≤≤≤=d d a ds把对面积的曲面即直接化作对θϕ,的二重积分. 输入Clear[A,fa,ds]; A={x^3,y^3,z^3}; fa={x,y,z}/a; ds=a^2*Sin[u];Integrate[(A.fa/.{x->a*Sin[u]*Cos[v],y->a*Sin[u]*Sin[v], z->a*Cos[u]})*ds//Simplify,{u,0,Pi}, {v,0,2Pi}]输出为96 5122πa如果用高斯公式计算, 则化为三重积分()d v z y x ⎰⎰⎰Ω++2223, 其中Ω为2222a z y x ≤++.采用球坐标计算, 输入<<Calculus`VectorAnalysis` 执行后再输入SetCoordinates[Cartesian[x,y,z]]; (*设定坐标系*) diva=Div[A]; (*求向量场的散度*)Integrate[(diva/.{x->r*Sin[u]*Cos[v],y->r*Sin [u]*Sin[v],z->r*Cos[u]})*r^2Sin[u],{v,0,2Pi}, {u,0,Pi},{r,0,a}]输出结果相同.实验习题 1. 计算⎰⎰-6/02/0.sin sin ππydydx x x y2. 计算下列积分的近似值: (1)();cos 022dydx y x ⎰⎰-ππ(2)().sin 1010dydx e xy ⎰⎰3. 计算下列积分 (1)();23012dydzdx z y e x x z xz x -⎰⎰⎰+- (2)⎰⎰1010.)arctan(dydx xy4. 交换积分次序并计算下列积分 (1)()d ydx y x x⎰⎰30922cos . (2) .20422dxdy e yx ⎰⎰5. 用极坐标计算下列积分: (1) ;10122dydx y x yx ⎰⎰+ (2) .13/3/22dxdy yx y y y ⎰⎰-+6. 用适当方法计算下列积分:(1)(),2/3222dv zy x z⎰⎰⎰Ω++ 其中Ω是由22y x z +=与1=z 围成;(2),)(224dv z y x++⎰⎰⎰Ω其中Ω是.1222≤++z y x7. 求()ds z y x f L⎰,,的近似值. 其中(),51,,33y x z y x f ++=,路径L :3/,2t y t x ==,.20,≤≤=t t z8. 求⎰L dr F ., 其中().0,sin cos ,121322π≤≤+=+++=t tj ti t r j y i x F 9. 用柱面坐标作图命令作出xy z =被柱面122=+y x 所围部分的图形,并求出其面积.9710. 求曲面积分,22zdxdy y x⎰⎰∑其中∑为球面2222a z y x =++的下半部分的下侧.11. 求曲面积分⎰⎰∑++zdS y x ,其中∑为球面2222a z y x=++上)0(a h z <<≥的部分.实验3 最小二乘拟合（基础实验）实验目的 了解曲线拟合问题与最小二乘拟合原理. 学会观察给定数表的散点图, 选择 恰当的曲线拟合该数表.最小二乘拟合原理 给定平面上的一组点,,,2,1),,(n k y x k k = 寻求一条曲线),(x f y =使它较好的近似这组数据, 这就是曲线拟合. 最小二乘法是进行曲线拟合的常用方法.最小二乘拟合的原理是, 求),(x f 使∑=-=nk k kyx f 12])([δ达到最小. 拟合时, 选取适当的拟合函数形式),()()()(1100x c x c x c x f m m ϕϕϕ+++=其中)(,),(),(10x x x m ϕϕϕ 称为拟合函数的基底函数.为使δ取到极小值, 将)(x f 的表达式 代入, 对变量i c 求函数δ的偏导数, 令其等于零, 就得到由1+m 个方程组成的方程组, 从中 可解出).,,2,1,0(m i c i =基本命令1.求数据的拟合函数的命令Fit 拟合函数Fit[ ]的基本格式为Fit[data,funs,vars],其中,data 是数据, vars 为变量(可以是多个变量), funs 为1+m 个以vars 为变量的基底函数. 其 输出结果是以基底函数(funs)的线性组合形式为拟合函数的最佳拟合函数(最小二乘估计的 结果). Fit 命令既可以作曲线拟合, 也可以作曲面拟合. 这里只讨论曲线拟合问题.曲线拟合时的数据格式为}}.,{,},,{},,{{2211n n y x y x y x下面是作曲线拟合时常用的几种拟合函数的形式Fit[data,{1,x},x] 用线性函数bx a +拟合数据data.Fit[data,{1,x,x^2},x] 用二次函数2cx bx a ++拟合数据data.Fit[data,Table[x^i,Table[x^i,{i,0,n}],x] 用x 的n 次多项式拟合数据data. 2.多项式拟合函数PolynomialFitMathematica 在程序包NumericalMath 中提供了多项式拟合函数PolynomialFit, 其基本格 式为PolynomialFit[data,n]它按最小二乘法构造n 次多项式函数拟合数据data. 例如,输入<<NumericalMath`PolynomialFit`98 p=PolynomialFit[{1,4,9,16,25,36,49},3]则输出FittingPolynomial[< >, 3]这里虽然没有给出拟合多项式的解析表达式, 但在计算机中已经存在. 因此可以用来计算函 数的近似值. 输入p[10] (*计算)10(f 的近似值*) 就得到函数的近似值100. 如果要拟合多项式的解析表达式, 输入Expand[p[x]]则输出321515x .0x .1x 1077636.11010543.7++⨯+⨯---3.去掉矩阵中非数值列的命令DropNonNumericColumn 如果矩阵M 中有非数值的列, 可先输入调用软件包命令<<Statistics\DataManipulation.m执行以后, 再输入DropNonNumericColumn[M]则在输出的矩阵中已经把含有非数值的列去掉.4.在Mathematica 中作曲线拟合的一般步骤在Mathematica 中作曲线拟合, 可按以下步骤进行:(1)用ListPlot[数据]作散点图, 观察曲线的分布形状, 确定基底函数; (2)用Fit[ ]命令求拟合函数; (3)用Plot[ ]命令作拟合曲线图;(4)最后用Show[ ]命令把散点图与拟合曲线图放在同一坐标系内, 观察拟合效果.实验举例曲线拟合 例 3.1 (教材 例 3.1) 为研究某一化学反应过程中温度)(C x 对产品得率(%)y 的影响, 测得数据如下:x 100 110 120 130 140 150 160 170 180 190 y 45 51 54 61 66 70 74 78 85 89试求其拟合曲线.输入点的坐标, 作散点图, 即输入b2={{100,45},{110,51},{120,54},{130,61},{140,66},{150,70},{160,74},{170,78},{180,85},{190,89}};fp=ListPlot[b2]99通过观察发现散点基本位于一条直线附近, 可用直线拟合. 输入Fit[b2,{1,x},x] (*用Fit 作拟合, 这里是线性拟合*)则输出拟合直线-2.73939+0.48303x作图观察拟合效果. 输入gp=Plot[%,{x,100,190},PlotStyle->{RGBColor[1,0,0]},DisplayFunction->Identity]; (*作拟合曲线的图形*)Show[fp,gp,DisplayFunction->$DisplayFunction](*显示数据点与拟合曲线*)例3.2 (教材 例3.2) 给定平面上点的坐标如下表:3627.100253.99493.70978.74337.69378.55687.53057.51234.59.08.07.06.05.04.03.02.01.0y x 试求其拟合曲线.输入data={{0.1,5.1234},{0.2,5.3057},{0.3,5.5687},{0.4, 5.9378},{0.5,6.4337},{0.6,7.0978},{0.7,7.9493},{0.8,9.0253},{0.9,10.3627}};pd=ListPlot[data];观察发现这些点位于一条抛物线附近. 用抛物线拟合, 即取基底函数.,,12x x 输入f=Fit[data,{1,x,x^2},x]则输出100 5.30661-1.83196x+8.17149x 2再输入fd=Plot[f,{x,0,1},DisplayFunction->Identity]; Show[pd,fd,DisplayFunction->$DisplayFunction]则输出平面上的点与拟合抛物线的图形（图3.2）.下面的例子说明Fit 的第二个参数中可以使用复杂的函数, 而不限于2,,1x x 等. 例3.3 (教材 例3.3) 使用初等函数的组合进行拟合的例子. 先计算一个数表. 输入ft=Table[N[1+2Exp[-x/3]],{x,10}]则输出{2.43306,2.02683,1.73576,1.52719,1.37775,1.27067,1.19394,1.13897,1.09957,1.07135}然后用基函数)exp(),3/exp(,sin ,1x x x --来做曲线拟合. 输入Fit[ft,{1,Sin[x],Exp[-x/3],Exp[-x]},x]则输出拟合函数][1022045.2.21044089.4.1163/15x Sin e e x x ----⨯++⨯-其中有些基函数的系数非常小, 可将它们删除. 输入Chop[%]则输出3/.2.1x e -+实际上,我们正是用这个函数做的数表.注:命令Chop 的基本格式为Chop[expr,δ]其含义是去掉表达式expr 的系数中绝对值小于δ的项,δ的默认值为1010-.实验4 水箱的流量问题（综合实验）实验目的 掌握应用最小二乘拟合原理分析和解决实际问题的思想和方法,能通过观察测试数据的散点图,建立恰当的数学模型,并用所学知识分析和解决所给问题.问题 (1991年美国大学生数学建模竞赛的A 题. 问题中使用的长度单位为E(英尺, 1 E=30.24cm), 容积单位是G(加仑, 1 G=3.785L)).101某些州的用水管理机构需估计公众的用水速度(单位:G/h)和每天的总用水量. 许多供水单位由于没有测量流入或流出量的设备, 而只能测量水箱中的水位(误差不超过5%). 当水箱水位低于水位L 时, 水泵开始工作将水灌入水箱, 直至水位达到最高水位H 为止. 但是依然无法测量水泵灌水流量, 因此, 在水泵工作时无法立即将水箱中的水位和水量联系起来. 水泵一天灌水1~2次, 每次约2h. 试估计在任一时刻(包括水泵灌水期间) t 流出水箱的流量),(t f 并估计一天的总用水量.表1给出了某镇某一天的真实用水数据. 水箱是直径为57E, 高为40E 的正圆柱体. 当水位落到27E 以下, 水泵自动启动把水灌入水箱; 当水位回升至35.5E 时, 水泵停止工作.模型假设(1) 影响水箱流量的唯一因素是该区公众对水的普通需求. 所给数据反映该镇在通常情况下一天的用水量, 不包括任何非常情况, 如水泵故障、水管破裂、自然灾害等. 并且认为水位高度、大气情况、温度变化等物理因素对水的流速均无直接影响;(2) 水泵的灌水速度为常数;(3) 从水箱中流出水的最大流速小于水泵的灌水速度. 为了满足公众的用水需求不让水箱中的水用尽, 这是显然的要求;(4) 因为公众对水的消耗量是以全天的活动(诸如洗澡、做饭、洗衣服等)为基础的, 所以,可以认为每天的用水量分布都是相似的;(5) 水箱的水流量速度可用光滑曲线来近似.问题分析与模型建立为方便起见,记V 表示水的容积;i V 表示时刻i t (单位:h)水的容积;)(t f 表示流出水箱的水的流速(单位;G/h),它是时间的函数;p 表示水泵的灌水速度(G/h).先将表1中数据作变换, 时间单位用小时(h), 水位高转换成水的体积(,2h r V π=单位:G 481.7E 1,G 1033= ). 输入tt={0,3316,6635,10619,13937,17921,21240,25223, 28543,32284,35932, 39332,39435,43318,46636,49953,53936,57254,60574,64554,68535, 71854,75021,79254,82649,85968,89953,93270}/3600//Nvv=Pi*(57/2)^2*{3175,3110,3054,2994,2947,2892, 2850,2795,2752,2697, no_data,no_data,3550,3445,3350,3260,3167,3087,3012,2927,2842,2767,2697,no_data,no_data,3475,3397,3340}*10^(-2)*7.481/10^3//N则输出下表.由于要求的是水箱流量与时间的关系, 因此须由表2的数据计算出相邻时间区间的中点及在时间区间内水箱中流出的水的平均速度.平均流速=(区间左端点的水量-区间右端点的水量)/时间区间长度输入tt1=Table[(tt[[i+1]]+tt[[i]])/2,{i,27}]vv1=Table[(vv[[i]]-vv[[i+1]])/(tt[[i+1]]-tt[[i]]),{i,27}]则输出下表模型求解为了作出时间tt1与平均水流量vv1之间的散点图, 先输入调用统计软件包的命令<<Statistics\DataManipulation.m执行以后再输入102103Clear[L];L=Transpose[DropNonNumericColumn[{tt1,vv1*10^3}]](*命令中vv1*10^3,使平均水流量vv1的单位变为G/h*) g1=ListPlot[L]则输出图4.1图中空白区域为泵水时间. 从中可以看出数据分布不均匀. 我们采用8阶多项式进行拟合. 输入ft=Fit[L,Table[t^i,{i,0,8}],t]则输出8765432t 00024547.0t 0248438.0t 01085.1t 1138.21t 101.240t 79.1468t 62.4690t 96.78394.16281+-+-+-+-这就是流出水箱的水的流速关于时间t 的函数)(t f . 为作出其拟合曲线图, 输入fg=Plot[ft,{t,0,26},DisplayFunction->Identity]; Show[g1,fg,DisplayFunction->$DisplayFunction]则输出图4.2.求解结果将460556.0=t h 和460556.24=t h 代入到水的流速拟合函数),(t f 我们得到这两时刻的流速分别近似为13532.5G/h 和13196.1G/h,相差仅2.48587%, 从而可以认为)(t f 能近似表达一天的用水流量.于是, 一天里的用水总量近似地等于函数)(t f 在24小时周期内的积分.输入Integrate[ft,{t,0.46,24.46}]104 则输出 336013.G若按常规每1000人的用水量为105000G/d, 因此估计出这个地区大约有3200人.模型评价该模型数学概念简单, 并且容易实现, 任意时刻从水箱中流出水的速度都可通过该模型计算出来, 可以推测速度. 但数据太少, 只能参照一天的数据. 另外, 如果知道水泵的灌水速度, 就能更准确地估算水泵灌水期间水的流速.实验报告某装饰材料商店欲以每瓶2元的成本价购进一批彩漆. 一般来说, 随着彩漆售价的提高,预期销售量将减少, 对此进行了估算, 见下表.202225282932343841/00.650.500.550.400.450.300.350.200.2/万瓶预期销售量元售价为了尽快收回资金并获得较多的赢利, 装饰材料商店打算做广告. 投入一定的广告费后,销售量将有一个增长, 可由销售增长因子来表示. 例如, 投入4万元的广告费, 销售增长因子为1.95. 即销售量将是预期销售量的1.95倍. 根据经验, 广告费与销售增长因子的关系见下表.8.195.100.295.185.170.140.100.1706050403020100)(销售增长因子元广告费用试确定装饰材料商店的最佳营销策略, 即确定彩漆售价和广告费投入使得预期的利润最大?实验5 线性规划问题（综合实验）实验目的 通过建立投资收益和风险问题的线性规划模型, 掌握利用线性规划理论建立实际问题的数学模型的思想和方法. 掌握用Mathematica 求解线性规划问题的基本方法.基本命令1.约束最大与约束最小命令求解线性规划问题的命令为ConstrainedMax 与ConstrainedMin. 其的基本格式是:ConstrainedMax[f,{inequalities},{x,y,…}]在不等式或等式{inequalities}确定的可行区域上求线性目标函数f 的最大值, 约定变量{x ,y ,…}都大于或等于0;ConstrainedMin{f,{inequalities},{x,y,…}}在不等式或等式{inequalities}确定的可行区域上求线性目标函数f 的最小值, 约定变量{x ,y ,…}都大于或等于0.注：上面两个命令都有一个可选参数:Tolerance 允许误差 (默认值是610-).例如, 输入ConstrainedMin[1.5 x+2.5 y,{x+3 y>=3,x+y>=2},{x,y}]则输出{3.5,{x->1.5,y->0.5}}即当5.0y ,5.1x ==时, 函数取得最小值3.5. 在约束条件中可以使用等号, 但要用“= =”表105示. 例如,输入ConstrainedMax[5 x+3 y+2 z+4 t,{3 x+y+2 z+8 t==10,2 x+4 y+2 z+t==10},{x,y,z,t}]则输出{18,{x->3,y->1,z->0,t->0}}有时, 输出结果可能有些问题. 输入ConstrainedMax[3x+2y-1,{x<1,y<2},{x,y}]则输出{6,{x->1,y->2}}即当2,1==y x 时, 函数取最大值6.注: 约束条件使用严格不等号, 结果仍旧取在边界上. 输入ConstrainedMax[x+y,{x+y<=15},{x,y}]则输出{15,{x->15,y->10}}这个问题有无穷多最优解, 这里只给出其中之一, 而且没有给出任何提示信息.前面的例题总是给出一个最优解, 属于正常情况.下面的例子是非正常的情况. 例如, 输入ConstrainedMax[x+y,{x-y>=0,3x-y<=-3},{x,y}]则在输入行的下面给出提示ConstrainedMax::nsat:The specified constraints cannot be satisfied.并输出ConstrainedMax[x+y,{x-y>=0,3x-y<=-3},{x,y}]其含义是: 没有可行解, 因此没有最优解. 然后返回投资的收益和风险问题.输入ConstrainedMax[2x+y,{x-y>=-1,-0.5x+y<=2},{x,y}]在输入行的下面给出提示ConstrainedMax::nbddt:Specified domain appearsunbounded,with tolerance1.’*^-6.并输出{∞,{x->Indeterminate,y->Indeterminate}}其含义是: 可行区域无界, 问题没有最大值, 或说最大值是无穷大. 然后返回投资的收益和风险问题.2.线性规划命令LinearProgramming当自变量和约束不等式较多时, 用ConstrainedMax 或ConstrainedMin 求解就比较麻烦. 此时, 可将目标函数和约束条件用向量或矩阵表示, 然后使用LinearProgramming. 其基本格式为LinearProgramming[c,m,b]其中c 是行向量, b 是列向量, m 是矩阵, 自变量用x 表示, 使用该命令, 则在满足不等式b mx ≥且0≥x 的可行区域中, 求出函数cx 的最小值点x .注: 实际输入时, b 仍以行向量表示. 此外, 这个命令也有可选参数Tolerance, 其含义与前面的说明相同.例如, 用约束最小命令计算, 输入ConstrainedMin[2x-3y,{x+y<10,x-y>2,x>1},{x,y}]则输出{0,{x->6,y->4}}改为用线性规划命令计算, 输入LinearProgramming[{2,-3},{{-1,-1},{1,-1},{1,0}},{-10,2,1}]。
Mathematica数学实验报告 实验三
数学实验报告实验三学院：数学与统计学院班级：信息与计算科学(1)班姓名：郝玉霞学号：201171020107实验三一、实验名：最佳分数近似值二、实验目的：研究怎样用分数近似值去给定的无理数作最佳逼近。
“最佳”就是既要误差小，又要分母小。
我们首先需要对“最佳”定出具体而明确的标准，还要寻找一个求最佳分数近似值的简单易行的算法。
三、实验环境：学校机房，Mathematica 软件。
四、实验的基本理论和方法：1、根据高中数学及大学数学中所学内容，经过分析研究，得出基本结论，利用Mathematica 来进行验证，并寻找一个求最佳分数近似值的简单易行的算法。
2、计算圆周率π“连分数展开”方法，并且利用特定的函数来展开其他数。
五、实验的内容和步骤实验步骤： 1、计算对数值对给定的正实数b ，N 且b ≠1，要求对数值a=N b log ,也就是求实数a 使a b =N ，如果能找到整数p ，q 使q pN b≈,则N b qp ≈，N b log qp≈，以lg2为例：由102=1024≈1000=310可得lg2≈103=0.3，再要提高精确度，就要找出更大的q 使q2更接近10的某个幂q10，也就是使p q32更接近于1。
练习题1：让q 依次取遍1到10000的所有的正整数，对每一个q ，按如下的递推法则求出一个正整数p=p(q)使实数p qq 102)(=λ最接近于1：q=1时，p(1)=0,λ(1)=01102=2.设已对q 求出p(q)和λ(q)，计算2λ(q)，如果2λ(q)<10,则取p(q+1)=p(q),λ(q+1)=2λ(q),如果2λ（q ）≥10，则取p(q+1)=p(q)+1，λ(q+1)=10)(2q λ. 如果λ(q)比以前所有的λ(i)(11-≤≤q i )都更接近1，即|λ(q)-1|<|λ(i)-1|对所有3、Mathematica 中常用的展开数与多项式的函数的使用；的1≤i ≤q-1成立，就取qp都是最佳逼近lg2的的分数近似值，它们可以展开成小数近似值。
Mathematica实验报告
Mathematica 实验报告【实验名称】利用MA THEMA TICA 作图、运算及编程.【实验目的】1。
掌握用MA THEMATICA 作二维图形，熟练作图函数Plot 、ParametricPlot 等应用，对图形中曲线能做简单的修饰.2。
掌握用MATHEMA TICA 做三维图形，对于一些二元函数能做出其等高线图等，熟练函数Plot3D ，ParametricPlot 的用法。
3、掌握用MA THEMATICA 进行微积分基本运算：求极限、导数、积分等。
【实验原理】1．二维绘图命令：二维曲线作图：Plot[fx,{x ，xmin,xmax}],二维参数方程作图：ParametricPlot[｛fx ，fy}，{t ，tmin ，tmax}]2．三维绘图命令：三维作图plot3D ［f,{x ，xmin ，xmax}，｛y,ymin ，ymax}］，三维参数方程作图：ParameticaPlot3D[{fx,fy ，fz ｝，{t ，tmin,tmax ｝］【实验内容】（含基本步骤、主要程序清单及异常情况记录等)1。
作出函数)sin(22y x z +=π的图形. 步骤： z=Sin ［Pi Sqrt[x^2+y^2］];Plot3D ［z ，{x,-1，1},{y,—1，1}，PlotPoints →30,Lighting →True]2。
椭球面()⎪⎪⎩⎪⎪⎨⎧=∈⎪⎭⎫ ⎝⎛-∈==u z v u v u y v u x R R R R R R sin ,,,2,0,2,2,sin cos cos cos 332121πππ自行给定,作图. 步骤：ParametricPlot3D ［{4Cos[u ］Cos[v],3Cos ［u]Sin[v]，2Sin[u]｝，{u ，—Pi/2，Pi/2}，｛v,0，2Pi}]3.做出极坐标描绘的图形：)cos 1(4θ+=r步骤：r ［t_］:=4(1+Cos[t ］）；ParametricPlot ［{r ［t ］Cos[t]，r ［t ］Sin ［t]},｛t,0,2Pi}]【实验结果】结果1：结果2:结果3：【总结与思考】MATHEMATICA作图的常见错误：General::spell1: Possible spelling error，因为在MATHEMATICA中作图函数大小写有区别.由于拼写间要有空格，易导致错误。
用Mathematica进行求导运算
y ln ln x
In[2] : D Log[Log[x]],x
Out[2]
1
xLog[x ]
例：求下列函数的二阶导数
y x8
In[3] : D x 8,x,2
Out[3] 56x 6
y 1 x2 arctan x
In[4] :
f (x)
f (x0 )
f (x0 )(x x0 )
f
( x0 2
)
(x
x0
)2
L
f
(n) (x0 n!
)
(x
x0
)n
Rn
(x)
麦克劳林公式:
f (x)
f (0)
f (0)x
f (0) x2 L 2
f
(n) (0) n!
xn
Rn
(x)
用Mathematica进行级数运算
Out[1]
x
x2
2
x3
6
x4
12
O(x 5 )
练习7
将函数 y ex 在x=1处展开到x的4次幂 将函数 y sin x和y cos x 在x= 0处展开到x的y 5ex次幂
将函数 y ln(x 1) 在x=1处展开到x的3次幂 x
用Mathematica进行求导运算
在Mathematica 中，求函数的导数或偏导数的格式为：
D[ f , x]
表示f对x求一阶导数
D[ f ,{x, n}] 表示f对x求n阶导数
例：求下列函数的一阶导数
y x3 cos x
In[1] : D x 3 * Cos[x ],x
用Mathematica软件求函数偏导数与多元函数的极值
实验九 用Mathematica 软件求函数偏导数与多元函数的极值实验目的:掌握用Mathematica 软件求函数偏导数与全微分、多元函数的极值的语句和方法。
实验过程与要求:教师利用多媒体组织教学,边讲边操作示范。
实验的内容:一、求偏导数在Mathematica 系统中与求一元函数导数类似用D 函数求函数f 的偏导数，基本格式为：D[f ,{变量,n }] 给出对变量的n 阶偏导数.D[f ,变量1,变量2，…] 给出高阶混合偏导数.实验 求y x x z cos sin +=的两个一阶偏导数和四个二阶偏导数.解 In[1]:=Clear[x ,y ]In[2]:=f [x _,y _]:=Sin[x ]+x *Cos[y ]In[3]:=D[f [x,y ],x ]In[4]:=D[f [x,y ],y ]In[5]:=D[f [x,y ],{x ,2}]In[6]:=D[f [x,y ],{y ,2}]In[7]:=D[f [x,y ],x,y ]In[8]:=D[f [x,y ],y,x ]Out[3]=Out[4]=Out[5]=Out[6]=Out[7]=Out[8]=二、求全微分在Mathematica 系统中与求一元函数微分类似用Dt 函数求函数f 的全微分，基本格式为：Dt[f ]实验 求函数206933+-+-+=y x xy y x z 的全微分.解 In[9]:=Dt[x ^3+y ^3-x *y +9x -6y +20]Out[9]=三、求多元函数的极值在Mathematica 系统中与求一元函数极小值类似用FindMinimum 函数求多变量函数f 的极小值，基本格式为：FindMinimum [f ,{x,x 0}，{y, y 0},…]其中{ x 0，y 0，…}为初始值，表示求出的是f 在（x 0，y 0，…）附近的极小值.因此，一般需借助于Plot3D 函数先作出函数的图象，由图象确定初始值，再利用FindMinimum 求出f 在（x 0，y 0，…）附近的极小值.仍用FindMinimum 函数求函数的极大值，基本格式为：FindMinimum [-f,{x,x 0}，{y, y 0},…]其中{ x 0，y 0，…}为初始值，表示求出的是-f 在（x 0，y 0，…）附近的极小值，设为W ，实际上间接地求出了f 在（x 0，y 0，…）附近的极大值，为-W.实验 求函数206922+-+-+=y x xy y x z 的极值.解 In[10]:=Clear[f,x,y ]In[11]:=FindMinimum[x ^2+y ^2+9*x -x *y-6y +20,{x ,-4},{y,-4}]In[12]:=Plot3D[x ^2+y ^2+9*x -x*y-6y +20,{x ,-4,5},{y ,-4,5}]Out[11]=表示z 在x =-4,y =1处取得极小值-1该函数无极大值.图形如图实验)ln()()4()3()2()1(.122222y x y x z yx yx z y x z yx e z xy +-=-+=-=+=求下列函数的偏导数：)ln()()4()arcsin()3()2()1(.2222y x y x z xy z y x z e z yx ++====-求下列函数的全微分：x y x y x z y x y x z 933)2()4)1(.3223322-++-=---=（求二元函数的极值：。
mathematica实验报告(符号计算)
1.表达式的运算
（1）化简： ；
（2）展开多项式： ；
（3）分解因式： ；
2.求函数的极限：（1） ；（2） ；（3） .
3.求导数：（1） ，求 ；（2） ，求 .
4.求积分： .
5.将 在 点，展开至 。
6.求和式与积式：（1） ；（2） .
7.求解方程 .
8.求微分方程：
四、程序、命令与结果
2.运行结果（）A准确，表现效果好；B正确；C部分结果不准确；D有较严重错误.
3.其它问题______________________________________________________________________.
4.综合评定（）A优秀；B良好；C合格；D不合格；E有明显抄袭或雷同现象.
结果：
（2）命令：Limit[(Tan[x])^Tan[2*x],xPi/4]
结果：
（3）命令：Limit[Exp[1/x],Direction0]
结果：
三、（1）命令：
结果：
（2）命令：
结果：
成绩评定：1.程序、命令（）A准确、简洁、效率高；B命令基本准确，ቤተ መጻሕፍቲ ባይዱ有少量问题；C部分命令有问题；
D许多命令都有问题或错误.
一、
（1）命令：
P=(x-2)*(x^2+2*x+4)+(x+5)*(x^2-5*x+25);
Simplify[P]
结果：
（2）命令：P=(a+b)^3;
Expand[P]
结果：
（3）命令：P=x^5-x;
Factor[P]
结果：
二、
（1）命令：Limit[((x+m)/(x-n))^x,xInfinity]
数学实验三 软件Mathematica求导数全微分
Dxyt[y_,x_,t_]:=D[y,t]/D[x,t]
自定义函数用于求参数方程所确 定的导数
例：求下列函数的一阶导数
y x3 cos x
In[1] : D x 3 * Cos[x ],x
Out[1] 3x 2Cos[x ] x 3Sin[x ]
y ln ln x
In[2] : D Log[Log[x]],x
命令
D[f[x],x] D[f[x],{x,n}]
功能 计算一元函数导数df/dx 计算一元函数高阶导数f(n)(x)
D[f,{x,n},{y,m}]
求函数f对x的n阶，对y的m阶混 合偏导数
Dt[f]
求函数f的全微分
DFxy[f_,x_,y_]:=Solve[D[f,x]==0, 自定义函数用于隐函数求导 y′[x]]
学生实验
基础操作
用mathematica求下列函数的导数
y e4x
y axex
y x 1 x 1 x
y sin x2
y (x 1 x2 )n y ln tan x
应用部分
• 将一物体垂直上抛，其运动方 s 10t ，1 g试t 2 求： １）物体从t=1秒到t=2秒的平均速度；2 ２）物体从t=1秒到t=1+△t秒的平均速度 ２）物体在t=1时的瞬时速度； ３）物体从t秒到t+△t秒的平均速度； ４）物体在任意t秒时的瞬时速度。
某公司在推销一种产品个月后，每月销售额（千元）可表示为
S(t) 2t3 40t2 220t 160
1）分别求1个月，4个月，6个月，9个月，20个月后的每月销售额； 2）求变化率 S(t) 3）分别求在 t 1, 4,6,9,12 处的变化率； 4）解释该公司的CEO为什么不必为6月份的销售额下降而发愁。
Mathematica数学实验——极限和导数
Mathematica数学实验——极限和导数教师指导实验4实验名称：极限和导数的运算⼀、问题：求⼀元函数的极限和导数。
⼆、实验⽬的：学会使⽤Mathematica 求数列和⼀元函数的极限（包括左极限、右极限），会求⼀元函数的导数，及利⽤导函数求原函数的单调区间和极值。
三、预备知识：本实验所⽤的Mathematica 命令提⽰1、Limit[f,x →x 0] 求函数f(x) 在x →x 0时的极限；2、Limit[f,x →x 0,Direction →-1] 求函数f(x) 在x →x 0时的右极限；Limit[f,x →x 0,Direction →1] 求函数f(x) 在x →x 0时的左极限； 3、D[f, var] 求函数f(x) 对⾃变量var 的导数；SetAttributes[k,Constant] 设定k 为常数；4、FindMinimum[f, {x, x 0}] 从x 0出发求函数f(x)的⼀个极⼩值点和极⼩值。
四、实验的内容和要求：1、求数列的极限1lim 1nn n →∞??+ 、11lim (1)nn i i i →∞=+∑；2、求函数的极限0sin lim x xx →、/2lim tan x x π→+；1lim (1)x x x e →∞-3、求下列函数的导数；sin cos n x nx ?、2cos ln x x ?、2(sin )(cos2)f x f x +4、求函数2()2ln f x x x =-的导数，求其单调区间和极值。
五、操作提⽰1、求数列的极限1lim 1nn n →∞+ 、11lim (1)nn i i i →∞=+∑；In[1]:= Limit[?n11+n ,n->Infinity]Out[1]= e In[2]:= Limit[∑ni=11i(i+1),n->∞] Out[2]= 12、求函数的极限0sin lim x xx→、/2lim tan x x π→+；1lim (1)x x x e →∞-In[3]:= Limit[Sin[x]x,x->0]Out[3]= 1In[4]:= Limit[Tan[x],x->Pi/2,Direction->-1] Out[4]= -∞ In[5]:= Limit[x(E^1 x-1),x->Infinity] Out[5]= 13、求下列函数的导数；sin cos nx nx ?、2cos ln x x ?、2(sin )(cos2)f x f x +In[6]:= D[Sin[x]^n Cos[nx],x] Out[6]= nCos[nx]Cos[x]Sin[x]-1+nIn[7]:= ?x (Cos[x]^2 Log[x])（注:?x 可以在基本输⼊输出模板中输⼊）Out[7]=2Cos[x]-x2Cos[x]Log[x]Sin[x] In[8]:= D[f[Sin[x]^2]+f[Cos[2x]]]Out[8]= -2Sin[2x]f ’[Cos[2x]]+2Cos[x]Sin[x]f ’[Sin[x]2]4、求函数2()2ln f x x x =-的导数，求其单调区间和极值。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第五题 Clear[f]; D[f[x^2],x]
Clear[f]; D[f[x^2],{x,2}]
Clear[f]; D[f[x]^2,x]
Clear[f]; D[f[x]^2,{x,2}]
Clear[f]; D[Log[f[x]],x]
Clear[f]; D[Log[f[x]],{x,2}]
Clear[f]; D[f[E^x+E^f[x]],x]
【实验原理】
1.求导数命令D与求微分命令Dt. D[f,x]给出关于的导数,而将表达式中的其他变量看作常量,因此,
如果是多元函数,则给出关于的偏导数. D[f,{x,n}]给出关于的阶导数或者偏导数.
D[f,x,y,z]给出关于,,的混合偏导数. Dt[f,x]给出关于的全导数,将表达式中的其他变量都看作的函数. Dt[f]给出的微分.如果是多元函数,则给出的偏微分. 即使表达式是抽象函数,上述命令也可以给出相应的正确结果,当然 是一些抽象符号. 命令D的选项Nonconstants一>{…}指出{…}内的字母是的函数. 命令Dt的选项Constants一>{…}指出{…}内的字母是常数. 2.解方程或方程组的命令Solve. 解方程命令的格式为 Solve[f[x]＝＝0,x] 解方程组命令的格式为 Solve[{f[x,y]＝＝0,g[x,y]＝＝0},{x,y}] 执行命令后给出方程或方程组关于指定变量的解.方程中的等号要 用双等号“＝＝”.如果是方程组,要用大括号将所有方程括起来,各方 程之间用逗号隔开. Solve的输出形如{{x->a}}.因此,为了取出其中解的表达式a,要使 用取出集合中元素的命令[[ ]]. 例如,解方程组,输入 Solve[{x+y＝＝a,x-y＝＝b},{x,y}] 输出为 3.循环语句Do. 循环语句Do的基本形式是 Do[表达式,循环变量的范围] 表达式中一般有循环变量,有多种方法说明循环变量的取值范围.最 完整的形式是 Do[表达式,{循环变量名,最小值,最大值,增量}]. 当省略增量时,默认增量为1.省略最小值时,默认最小值为1. 例如输入
4．拉格朗日中值定理 例3．9 函数在区间[1,2]上满足拉格朗日中值定理的条件．因此
存在使．可以验证这个结论的正确性．输入 Clear[f]； f[x_]:=1/x^4； Solve[D[f[x],x]＝＝f[2]-f[1],x]//N
【实验结论】（结果）
1.利用求导命令和求微积分命令可以容易的求出导数和微积分; 2.每一个导数和微积分都可以算出来; 3.实验很成功.
{{x 0.5 (a+b)}} 第四题
Clear[f]; f[x_]:=Sin[x]; Clear[F] ; F[x_]:=x+Cos[x]; Solve[D[f[x],x]/D[F[x],x]==(f[Pi/2]f[0])/(F[Pi/2]-F[0]),x]//N
Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found. {{Null (x 0.533458)}}
{x,-1,1},PlotStyle{RGBColor[0,0,1],RGBColor[1,0,0]}] 例3.2 作函数的图形和在处的切线. 输入 Clear[f]； f[x_]=2x^3+3x^2-12x+7； plotf=Plot[f[x],{x,-4,3},DisplayFunction Identity]； plot2=Plot[f'[-1]*(x+1)+f[-1],
第七题 deq1=D[Log[x]+Exp[-y[x]/x]
==E,x]
Solve[deq1,y ' [x]]
deq1=D[ArcTan[y[x]/x]==Log[Sqrt[x^2+y[x]^2]] ,x]
Solve[deq1,y ' [x]]
第八题 D[Sin[t]^3,t]/D[Cos[t]^3,t]
Do[Print[Sin[n*x],{n,1,10}] 则在屏幕上显示sin[x],sin[2x],……,sin[10x]等10个函数.
【实验环境】
Mathematic 4
二、实验内容：
【实验方案】
1.导数概念与导数的几何意义; 2.求函数的导数与微分;
3.求隐函数的导数,由参数方程定义的函数的导数; 4.拉格朗日中值定理.
{x,-4,3},PlotStyle GrayLevel[0.5],DisplayFunction Identity]； Show[plotf,plot2,DisplayFunction $DisplayFunction].
2．求函数的导数与微分. 例3.3 求函数的一阶导数 D[x^n,x] 例3.4 求函数 的一阶导数,并求 diff[x_]=D[Sin[a*x]*Cos[b*x],x] diff[1/(a+b)] 例3.5 求函数的1阶到11阶导数. 输入 Clear[f]; f[x_]=x^10+2(x-10)^9; D[f[x],{x,2}] Do[Print[D[f[x],{x,n}]],{n,1,11}] 或输入 Table[D[f[x],{x,n},{x,11}] 例3.6 求函数与例3.4中函数的微分.
D[%,t]/D[Cos[t]^3,t]//Simplify
D[6t^2/(1+t^3),t]/D[6t/(1+t^3),t]
D[%,t]/D[6t/(1+t^3),t]//Simplify
附录2：实验报告填Байду номын сангаас说明
1.实验项目名称：要求与实验教学大纲一致. 2.实验目的：目的要明确,要抓住重点,符合实验教学大纲要求. 3.实验原理：简要说明本实验项目所涉及的理论知识. 4.实验环境：实验用的软、硬件环境. 5.实验方案（思路、步骤和方法等）：这是实验报告极其重要的内容.概括整个实验过程.
4实验结论正确.
成 绩：
指导教师签名： 批阅日期：
附录1：源 程 序
第一题
实验三 导 数
Clear[f];f[x_]:=Log[Sin[x]]; Solve[D[f[x],x]==0,x]//N Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found. {{Null (x 1.5708)}} 第二题
7.实验结论（结果）：根据实验过程中得到的结果,做出结论. 8.实验小结：本次实验心得体会、思考和建议. 9.指导教师评语及成绩：指导教师依据学生的实际报告内容,给出本次实验报告的评价.
输入 Dt[Sin[2*x]] Dt[Sin[a*x]*Cos[b*x],Constants->{a,b}]//Simplify Dt[Sin[a*x]*Cos[b*x]]
3．求隐函数的导数,由参数方程定义的函数的导数. 例3.7 求由方程确定的隐函数的导数. 输入 deq1=D[2x^2-2x*y[x]+y[x]^2+x+2y[x]+1＝＝0,x]； Solve[deq1,y'[x]] deq2=Dt[2x^2-2x*y[x]+y[x]^2+x+2y[x]+1＝＝0,x] Solve[deq2,Dt[y,x]] deq3=D[deq1,x]； Solve[{deq1,deq3},{y'[x],y''[x]}]//Simplify 例3.8 求由参数方程确定的函数的导数 D[E^t*Sin[t],t]/D[E^t*Cos[t],t] D[%,t]/D[E^t*Cos[t],t]//Simplify
Clear[f]; D[f[E^x+E^f[x]],{x,2}]
第六题 Clear[f]; f[x_]=x*Sinh[x]; D[f[x],{x,100}]
Clear[f]; f[x_]=x^2*Cos[x]; D[f[x],{x,10}]
Clear[f]; f[x_]=x^2*Sin[2x]; D[f[x],{x,50}]
【实验小结】（收获体会）
1.用Mathematic 4求导数和求微积分很方便; 2.可以把实验内容和实际的计算联系起来.
三、指导教师评语及成绩：
评语
评语等级 优 良 中 及格 不及格
1.实验报告按时完成,字迹清楚,文字叙述流畅,逻辑 性强
2.实验方案设计合理
3.实验过程（实验步骤详细,记录完整,数据合理,分 析透彻）
【实验过程】（实验步骤、记录、数据、分析）
1．导数概念与导数的几何意义. 例3.1 用定义求函数的导数. 输入 Clear[g]； g[x_]=x^3-3x^2+x+1； Simplify[(g[x+h]-g[x])/h] %/.h 0 Limit[%%,h 0] Plot[{g[x],g'[x]},
对于验证性实验,要写明依据何种原理、操作方法进行实验,要写明需要经过哪几个步骤来实 现其操作.对于设计性和综合性实验,在上述内容基础上还应该画出流程图、设计思路和设计 方法,再配以相应的文字说明.对于创新性实验,应注明其创新点、特色.
6.实验过程（实验中涉及的记录、数据、分析）：写明具体实验方案的具体实施步骤,包括实验 过程中的记录、数据和相应的分析.
天水师范学院数学与统计学院
实验报告
实验项目名称 导 数
所属课程名称
数学实验
实 验 类 型 微积分实验
实 验 日 期 2011.10.05
班级 学号 姓名 成绩
一、实验概述：
【实验目的】
1. 深入理解导数与微分的概念,导数的几何意义. 2. 掌握用Mathematica求导数与高阶导数的方法. 3.深入理解和掌握求隐函数的导数,以及求由参数方程定义的函数 的导数的方法. 4. 完成数学实验报告,总结方法,增强数学思维能力.