04项目三 多元函数微积分学

实验五 矩阵运算与方程组求解 附录Ⅰ 大学数学实验指导书 项目三 多元函数微积分 实验1 多元函数微分学（基础实验） 实验目的 掌握利用Mathematica计算多元函数偏导数和全微分的 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 , 掌握计算二元 函数极值和条件极值的方法. 理解和掌握曲面的切平面的作法. 通过作图和观察, 理解二元 函数的性质、方向导数、梯度和等高线的概念. 基本命令 1.求偏导数的命令D 命令D既可以用于求一元函数的导数, 也可以用于求多元函数的偏导数. 例如: 求 对x的偏导数, 则输入D[f[x,y,z],x] 求 对y的偏导数, 则输入D[f[x,y,z],y] 求 对x的二阶偏导数, 则输入D[f[x,y,z],{x,2}] 求 对 的混合偏导数, 则输入D[f[x,y,z],x,y] ………… 2.求全微分的命令Dt 该命令只用于求二元函数 的全微分时, 其基本格式为 Dt[f[x,y]] 其输出的表达式中含有Dt[x],Dt[y], 它们分别表示自变量的微分dx,dy. 若函数 的表 达式中还含有其它用字符表示的常数, 例如a, 则Dt[f[x,y]]的输出中还会有Dt[a], 若采用选 项Constants->{a}, 就可以得到正确结果, 即只要输入 Dt[f[x,y],Constants->{a}] 3.在 平面上作二元函数 的等高线的命令ContourPlot 命令的基本格式为 ContourPlot[f[x,y],{x,x1,x2},{y,y1,y2}] 例如,输入 ContourPlot[x^2-y^2,{x,-2,2},{y,-2,2}] 则输出函数 的等高线图(图1.1). 该命令的选项比较多(详细的内容参见光盘中的实验案例库). 如选项Contours->15表示作15条等高线, 选项Contours->{0}表示只作函数值为0的等高线. 图1.1 实验举例 求多元函数的偏导数与全微分 例1.1 (教材 例1.1) 设 求 输入 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 设 求 和全微分dz. 输入 Clear[z];z=(1+x*y)^y; D[z,x] D[z,y] 则有输出 再输入 Dt[z] 则得到输出 例1.3 (教材 例1.2) 设 其中a是常数, 求dz. 输入 Clear[z,a];z=(a+x*y)^y; wf=Dt[z,Constants->{a}]//Simplify 则输出结果: (a+xy)-1+y(y2Dt[x,Constants->{a}]+ Dt[y,Constants->{a}](xy+(a+xy)Log[a+xy])) 其中Dt[x,Constants->{a}]就是dx, Dt[y,Constants->{a}]就是dy. 可以用代换命令“/.”把它们 换掉. 输入 wf/.{Dt[x,Constants->{a}]->dx,Dt[y,Constants->{a}]->dy} 输出为 (a+xy)-1+y(dxy2+dy(xy+(a+xy)Log[a+xy])) 例1.4 (教材 例1.3) 设 ,求 输入 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组成的方程组*) 则输出 其中D[u,x,NonConstants->{u,v}]表示u对x的偏导数, 而D[v,x,NonCosnstants->{u,v}]表示v 对x的偏导数. 类似地可求得u,v对y的偏导数. 微分学的几何应用 例1.5 求出曲面 在点(1,1)处的切平面、法线方程, 并画出图形. 解(1) 画出曲面的图形. 曲面的参数方程为 输入命令 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). 图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}]; 则输出切平面方程为 及相应图形(图1.3). 图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) 求曲面 在点 处的切平面方程, 并把曲面和它的切平面作在同一图形里. 输入 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.5）. 图1.5 多元函数的极值 例1.7 (教材 例1.5) 求 的极值. 输入 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}] 则分别输出所求偏导数和驻点: {{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 输出为判别式函数 的形式: (6+6x)(6-6y) 再输入 data={x,y,fxx,disc,f[x,y]}/.critpts; TableForm[data,TableHeadings->{None,{ "x ", "y ", "fxx ", "disc ", "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 易见,当 时 判别式disc=72, 函数有极大值31; 当 时 判别式disc=72, 函数有极小值-5; 当 和 时, 判别式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. 图1.6 从上图可见, 在两个极值点附近, 函数的等高线为封闭的. 在非极值点附近, 等高线不 封闭. 这也是从图形上判断极值点的方法. 注:在项目一的实验4中,我们曾用命令FindMinimum来求一元函数的极值, 实际上,也可 以用它求多元函数的极值, 不过输入的初值要在极值点的附近. 对本例,可以输入以下命令 FindMinimum[f[x,y],{x,-1},{y,1}] 则输出 {-5.,{x->1.,y->-2.36603×10-8}} 从中看到在 的附近函数 有极小值-5, 但y的精度不够好. 例1.8 求函数 在条件 下的极值. 输入 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}] 得到输出 再输入 f[x,y]/.extpts//Simplify 得到两个可能是条件极值的函数值 但是否真的取到条件极值呢? 可利用等高线作图来判断. 输入 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] 输出为 及图1.7. 从图可见，在极值可疑点 处, 函数 的等高线与曲线 (虚线)相切. 函数 的等高线是一系列同心圆, 由里向外, 函数值在增大, 在 的附近观察, 可以得出 取条件极大的结论. 在 的附近观察, 可以得出 取条件极小的结论. 图1.7 梯度场 例1.9 画出函数 的梯度向量. 解 输入命令 <{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.8 例1.10 在同一坐标面上作出 和 的等高线图( ), 并给出它们之间的关系. 解 输入命令 <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) 的梯度与等高线图; (b) 的梯度与等高线图; (c) 与 的等高线图; (d) 与 的梯度图. (a) (b) (c) (d) 图1.9 从上述图中可以看出它们的等高线为一族正交曲线. 事实上, 有 且 它们满足拉普拉斯方程 例1.11 (教材 例1.6) 设 作出 的图形和等高线, 再作出它的梯度向量gradf的图形. 把上述等高线和梯度向量的图形叠加在一起, 观察它们之间的关系. 输入调用作向量场图形的软件包命令 <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. 从图可以看到 平面上过每一点的等高线和梯度向量是垂直的, 且梯度的 方向是指向函数值增大的方向. 图1.10 例1.12 求出函数 的极值, 并画出函数 的等高线、驻点以及 的梯度向量的图形. 输入命令 <False,PlotPoints->100,Contours->{-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->{PointSize[0.03]}]; Show[conplot,fieldplot,critptplot]; 则得到 的最小值 以及函数的图形(图1.11). 图1.11 实验习题 1.设 求 2.设 求 3.设 求 4.试用例1.5的方法求 的极值. 5.求 在 条件下的极值. 6.作出函数 的等高线和梯度线的图形, 并观察梯度线与等高线的 关系. 实验2 多元函数积分学（基础实验） 实验目的 掌握用Mathematica计算二重积分与三重积分的方法; 深入理解曲线积分、曲面积分的 概念和计算方法. 提高应用重积分和曲线、曲面积分解决各种问题的能力. 基本命令 1. 计算重积分的命令lntegrate和NIntegrate 例如,计算 , 输入 Integrate[x*y^2,{x,0,1},{y,0,x}] 则输出 又如,计算 的近似值, 输入 NIntegrate[Sin[x*y^2],{x,0,1},{y,0,1}] 则输出 0.160839 注: Integrate命令先对后边的变量积分. 计算三重积分时,命令Integrate的使用格式与计算二重积分时类似. 由此可见, 利用 Mathematica计算重积分, 关键是确定各个积分变量的积分限. 2. 柱坐标系中作三维图形的命令CylindricalPlot3D 使用命令Cylindricalplot3D, 首先要调出作图软件包. 输入 <40] 则在球面坐标系中作出了该球面的图形. 4. 向量的内积 用“.”表示两个向量的内积. 例如,输入 vecl={al,bl,cl} vec2={a2,b2,c2} 则定义了两个三维向量, 再输入 vec1. vec2 则得到它们的内积 a1a2+b1b2+c1c2 实验举例 计算重积分 例2.1 (教材 例2.1) 计算 其中 为由 所围成的有界区域.先作出区域 的草图, 易直接确定积分限,且应先对 积分, 因此, 输入 Integrate[x*y^2,{y,1,2},{x,2-y,Sqrt[y]}] 则输出所求二重积分的计算结果 例2.2 (教材 例2.2) 计算 其中 为 如果用直角坐标计算, 输入 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]}] 则输出为 其中Erf是误差函数. 显然积分遇到了困难. 如果改用极坐标来计算, 也可用手工确定积分限. 输入 Integrate[(f[x,y]/.{x->r*Cos[t],y->r*Sin[t]})*r,{t,0,2 Pi},{r,0,1}] 则输出所求二重积分的计算结果 如果输入 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) 计算 , 其中 由曲面 与 围成. 先作出区域 的图形. 输入 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)）. (a) (b) 图2.1 考察上述图形, 可用手工确定积分限. 如果用直角坐标计算, 输入 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]}] 则输出 如果用球面坐标计算,输入 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]}] 则输出 这与柱面坐标的结果相同. 重积分的应用 例2.4 求由曲面 与 所围成的空间区域 的体积. 输入 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. 图2.2 首先观察到 的形状. 为了确定积分限, 要把两曲面的交线投影到 平面上输入 jx=Solve[f[x,y]==g[x,y],y] 得到输出 为了取出这两条曲线方程, 输入 y1=jx[[1,1,2]] y2=jx[[2,1,2]] 输出为 再输入 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, 由此可见, 是下半圆(虚线), 是上半圆,因此投影区域是一个圆. 图2.2 设 的解为 与 ,则 为 的积分限. 输入 xvals=Solve[y1==y2,x] 输出为 为了取出 , 输入 x1=xvals[[1,1,2]] x2=xvals[[2,1,2]] 输出为 这时可以作最后的计算了. 输入 Volume=Integrate[g[x,y]-f[x,y],{x,x1,x2}, {y,y1,y2}]//Simplify 输出结果为 例2.5 (教材 例2.4) 求旋转抛物面 在 平面上部的面积 先调用软件包, 输入 <r*Cos[t],y->r*Sin[t]}]; Integrate[z1*r,{t,0,2 Pi},{r,0,2}]//Simplify 则输出所求曲面的面积 例2.6 在 平面内有一个半径为2的圆, 它与 轴在原点 相切, 求它绕 轴旋转一周所得旋转体体积. 先作出这个旋转体的图形. 因为圆的方程是 它绕 轴旋转所得的圆环面的方程为 , 所以圆环面的球坐标方程是 输入 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]}] 得到这个旋转体的体积为 计算曲线积分 例2.7 (教材 例2.5) 求 , 其中 积分路径为 : 注意到,弧长微元 , 将曲线积分化为定积分,输入 Clear[x,y,z]; luj={t,t^2,3t^2}; D[luj,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) 求 , 其中 输入 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 求锥面 与柱面 的交线的长度. 先画出锥面和柱面的交线的图形. 输入 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. 图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. 图2.6 为了用线积分计算曲线的弧长, 必须把曲线用参数方程表示出来. 因为空间曲线的投影曲线的方程为 , 它可以化成 , 再代入锥面方程 , 得 因为空间曲线的弧长的计算公式是 , 因此输入 Clear[x,y,z]; x=Cos[t]^2; y=Cos[t]*Sin[t]; z=Cos[t]; qx={x,y,z}; Integrate[Sqrt[D[qx,t]. D[qx,t]]//Simplify, {t,-Pi/2,Pi/2}] 输出为 2Elliptice[-1] 这是椭圆积分函数. 换算成近似值. 输入 %//N 输出为 3.8202 计算曲面积分 例2.10 (教材 例2.7) 计算曲面积分 , 其中 为锥面 被柱面 所截得的有限部分. 注意到,面积微元 , 投影曲线 的极坐标方程为 将曲面积分化作二重积分,并采用极坐标计算重积分. 输入 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]}] 则输出所求曲面积分的计算结果 例2.11 计算曲面积分 其中 为球面 的外侧. 可以利用两类曲面积分的关系, 化作对曲面面积的曲面积分 . 这里 . 因为球坐标的体积元素 注意到在球面 上 , 取 后得到面积元素的表示式: 把对面积的曲面即直接化作对 的二重积分. 输入 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}] 输出为 如果用高斯公式计算, 则化为三重积分 , 其中 为 . 采用球坐标计算, 输入 <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. 计算 2. 计算下列积分的近似值: (1) (2) 3. 计算下列积分 (1) (2) 4. 交换积分次序并计算下列积分 (1) . (2) 5. 用极坐标计算下列积分: (1) (2) 6. 用适当方法计算下列积分: (1) 其中 是由 与 围成; (2) 其中 是 7. 求 的近似值. 其中 ,路径 : , 8. 求 , 其中 9. 用柱面坐标作图命令作出 被柱面 所围部分的图形,并求出其面积. 10. 求曲面积分 其中 为球面 的下半部分的下侧. 11. 求曲面积分 ,其中 为球面 上 的部分. 实验3 最小二乘拟合（基础实验） 实验目的 了解曲线拟合问题与最小二乘拟合原理. 学会观察给定数表的散点图, 选择 恰当的曲线拟合该数表. 最小二乘拟合原理 给定平面上的一组点 寻求一条曲线 使它较好的近似这组数据, 这就是曲线拟合. 最小二乘法是进行曲线拟合的常用方法. 最小二乘拟合的原理是, 求 使 达到最小. 拟合时, 选取适当的拟合函数形式 其中 称为拟合函数的基底函数.为使 取到极小值, 将 的表达式 代入, 对变量 求函数 的偏导数, 令其等于零, 就得到由 个方程组成的方程组, 从中 可解出 基本命令 1.求数据的拟合函数的命令Fit 拟合函数Fit[ ]的基本格式为 Fit[data,funs,vars], 其中,data是数据, vars为变量(可以是多个变量), funs为 个以vars为变量的基底函数. 其 输出结果是以基底函数(funs)的线性组合形式为拟合函数的最佳拟合函数(最小二乘估计的 结果). Fit命令既可以作曲线拟合, 也可以作曲面拟合. 这里只讨论曲线拟合问题. 曲线拟合时的数据格式为 下面是作曲线拟合时常用的几种拟合函数的形式 Fit[data,{1,x},x] 用线性函数 拟合数据data. Fit[data,{1,x,x^2},x] 用二次函数 拟合数据data. Fit[data,Table[x^i,Table[x^i,{i,0,n}],x] 用x的n次多项式拟合数据data. 2.多项式拟合函数PolynomialFit Mathematica在程序包NumericalMath中提供了多项式拟合函数PolynomialFit, 其基本格 式为 PolynomialFit[data,n] 它按最小二乘法构造n次多项式函数拟合数据data. 例如,输入 <, 3] 这里虽然没有给出拟合多项式的解析表达式, 但在计算机中已经存在. 因此可以用来计算函 数的近似值. 输入 p[10] (*计算 的近似值*) 就得到函数的近似值100. 如果要拟合多项式的解析表达式, 输入 Expand[p[x]] 则输出 3.去掉矩阵中非数值列的命令DropNonNumericColumn 如果矩阵M中有非数值的列, 可先输入调用软件包命令 <{RGBColor[1,0,0]}, DisplayFunction->Identity]; (*作拟合曲线的图形*) Show[fp,gp,DisplayFunction->$DisplayFunction] (*显示数据点与拟合曲线*) 则输出平面上的点与拟合抛物线的图形（图3.1）. 图3.1 例3.2 (教材 例3.2) 给定平面上点的坐标如下表: 试求其拟合曲线. 输入 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]; 则输出题设数据的散点图. 观察发现这些点位于一条抛物线附近. 用抛物线拟合, 即取基底函数 输入 f=Fit[data,{1,x,x^2},x] 则输出 5.30661-1.83196x+8.17149x2 再输入 fd=Plot[f,{x,0,1},DisplayFunction->Identity]; Show[pd,fd,DisplayFunction->$DisplayFunction] 则输出平面上的点与拟合抛物线的图形（图3.2）. 图3.2 下面的例子说明Fit的第二个参数中可以使用复杂的函数, 而不限于 等. 例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} 然后用基函数 来做曲线拟合. 输入 Fit[ft,{1,Sin[x],Exp[-x/3],Exp[-x]},x] 则输出拟合函数 其中有些基函数的系数非常小, 可将它们删除. 输入 Chop[%] 则输出 实际上,我们正是用这个函数做的数表. 注:命令Chop的基本格式为 Chop[expr, ] 其含义是去掉表达式expr的系数中绝对值小于 的项, 的默认值为 . 实验4 水箱的流量问题（综合实验） 实验目的 掌握应用最小二乘拟合原理分析和解决实际问题的思想和方法,能通过观察测试数据的散点图,建立恰当的数学模型,并用所学知识分析和解决所给问题. 问题 (1991年美国大学生数学建模竞赛的A题. 问题中使用的长度单位为E(英尺, 1 E=30.24cm), 容积单位是G(加仑, 1 G=3.785L)). 某些州的用水管理机构需估计公众的用水速度(单位:G/h)和每天的总用水量. 许多供水单位由于没有测量流入或流出量的设备, 而只能测量水箱中的水位(误差不超过5%). 当水箱水位低于水位L时, 水泵开始工作将水灌入水箱, 直至水位达到最高水位H为止. 但是依然无法测量水泵灌水流量, 因此, 在水泵工作时无法立即将水箱中的水位和水量联系起来. 水泵一天灌水1~2次, 每次约2h. 试估计在任一时刻(包括水泵灌水期间) t流出水箱的流量 并估计一天的总用水量. 表1给出了某镇某一天的真实用水数据. 水箱是直径为57E, 高为40E的正圆柱体. 当水位落到27E以下, 水泵自动启动把水灌入水箱; 当水位回升至35.5E时, 水泵停止工作. 表1 时间/s 水位 E 时间/s 水位 E 0 3316 6635 10619 13937 17921 21240 25223 28543 32284 35932 39332 39435 43318 3175 3110 3054 2994 2947 2892 2850 2795 2752 2697 泵水 泵水 3550 3445 46636 49953 53936 57254 60574 64554 68535 71854 75021 79254 82649 85968 89953 93270 3350 3260 3167 3087 3012 2927 2842 2767 2697 泵水 泵水 3475 3397 3340 模型假设 (1) 影响水箱流量的唯一因素是该区公众对水的普通需求. 所给数据反映该镇在通常情况下一天的用水量, 不包括任何非常情况, 如水泵故障、水管破裂、自然灾害等. 并且认为水位高度、大气情况、温度变化等物理因素对水的流速均无直接影响; (2) 水泵的灌水速度为常数; (3) 从水箱中流出水的最大流速小于水泵的灌水速度. 为了满足公众的用水需求不让水箱中的水用尽, 这是显然的要求; (4) 因为公众对水的消耗量是以全天的活动(诸如洗澡、做饭、洗衣服等)为基础的, 所以,可以认为每天的用水量分布都是相似的; (5) 水箱的水流量速度可用光滑曲线来近似. 问题分析与模型建立 为方便起见,记V表示水的容积; 表示时刻 (单位:h)水的容积; 表示流出水箱的水的流速(单位;G/h),它是时间的函数;p表示水泵的灌水速度(G/h). 先将表1中数据作变换, 时间单位用小时(h), 水位高转换成水的体积( 单位: ). 输入 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//N vv=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 时间/h 水量/ G 时间/h 水量/ G 0. 0.921111 1.84306 2.94972 3.87139 4.97806 5.9 7.00639 7.92861 8.96778 9.98111 10.9256 10.9542 12.0328 606.098 593.69 583. 571.546 562.574 552.074 544.057 533.557 525.349 514.849 no_data no_data 677.685 657.64 12.9544 13.87558 14.9822 15.9039 16.8261 17.9317 19.0375 19.9594 20.8392 22.015 22.9581 23.88 24.9869 25.9083 639.505 622.324 604.571 598.299 574.982 558.756 542.529 528.212 514.849 no_data no_data 663.367 648.477 637.593 由于要求的是水箱流量与时间的关系, 因此须由表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}] 则输出下表 表3 时间区间的中点值/h 平均水流量/ G/h 时间区间的中点值/h 平均水流量/ G/h 0.460556 1.38208 2.39639 3.41056 4.42472 5.43903 6.45319 7.4675 8.44819 9.47444 10.4533 10.9399 11.4935 12.4936 13.471 11.5953 10.3498 9.73471 9.48735 8.69649 9.48974 8.90086 10.1036 no_data no_data no_data 18.5833 19.6766 13.4151 14.429 15.4431 16.365 17.3789 18.4846 19.4985 20.3993 21.4271 22.4865 23.419 24.4335 25.4476 18.6466 16.0463 16.5697 15.5248 14.677 14.6733 15.5294 15.1898 no_data no_data no_data 13.4514 11.8095 模型求解 为了作出时间tt1与平均水流量vv1之间的散点图, 先输入调用统计软件包的命令 <Identity]; Show[g1,fg,DisplayFunction->$DisplayFunction] 则输出图4.2. 图4.2 求解结果 将 h和 h代入到水的流速拟合函数 我们得到这两时刻的流速分别近似为13532.5G/h和13196.1G/h,相差仅2.48587%, 从而可以认为 能近似表达一天的用水流量.于是, 一天里的用水总量近似地等于函数 在24小时周期内的积分. 输入 Integrate[ft,{t,0.46,24.46}] 则输出 336013.G 若按常规每1000人的用水量为105000G/d, 因此估计出这个地区大约有3200人. 模型 评价 LEC评价法下载LEC评价法下载评价量规免费下载学院评价表文档下载学院评价表文档下载 该模型数学概念简单, 并且容易实现, 任意时刻从水箱中流出水的速度都可通过该模型计算出来, 可以推测速度. 但数据太少, 只能参照一天的数据. 另外, 如果知道水泵的灌水速度, 就能更准确地估算水泵灌水期间水的流速. 实验 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 某装饰材料商店欲以每瓶2元的成本价购进一批彩漆. 一般来说, 随着彩漆售价的提高,预期销售量将减少, 对此进行了估算, 见下表. 为了尽快收回资金并获得较多的赢利, 装饰材料商店打算做广告. 投入一定的广告费后,销售量将有一个增长, 可由销售增长因子来表示. 例如, 投入4万元的广告费, 销售增长因子为1.95. 即销售量将是预期销售量的1.95倍. 根据 经验 班主任工作经验交流宣传工作经验交流材料优秀班主任经验交流小学课改经验典型材料房地产总经理管理经验 , 广告费与销售增长因子的关系见下表. 试确定装饰材料商店的最佳营销策略, 即确定彩漆售价和广告费投入使得预期的利润最大? 实验5 线性规划问题（综合实验） 实验目的 通过建立投资收益和风险问题的线性规划模型, 掌握利用线性规划理论建立实际问题的数学模型的思想和方法. 掌握用Mathematica求解线性规划问题的基本方法. 基本命令 1.约束最大与约束最小命令 求解线性规划问题的命令为ConstrainedMax