1、 附录 大学数学实验指导书项目三 多元函数微积分实验 1 多元函数微分学(基础实验)实验目的 掌握利用 Mathematica 计算多元函数偏导数和全微分的方法, 掌握计算二元函数极值和条件极值的方法. 理解和掌握曲面的切平面的作法. 通过作图和观察, 理解二元函数的性质、方向导数、梯度和等高线的概念.求多元函数的偏导数与全微分例 1.1 (教材 例 1.1) 设 求),(cos)sin(2xyz.,2yxzz输入Clearz;z=Sinx*y+Cosx*y2;Dz,xDz,yDz,x,2Dz,x,y则输出所求结果.例 1.2 设 求 和全微分 dz.,)1(yxzz,输入Clearz;z=(
2、1+x*y)y;Dz,xDz,y则有输出 1)1(12xyLogxyyy再输入Dtz则得到输出 11)()1( xyLogDtxyttxy例 1.3 (教材 例 1.2) 设 其中 a 是常数, 求 dz.,(az输入Clearz,a;z=(a+x*y)y;wf=Dtz,Constants-a/Simplify则输出结果:(a+xy)-1+y(y2Dtx,Constants-a+Dty,Constants-a(xy+(a+xy)Loga+xy)其中 Dtx,Constants-a就是 dx, Dty,Constants-a就是 dy. 可以用代换命令“/.”把它们换掉. 输入wf/.Dtx,C
3、onstants-a-dx,Dty,Constants-a-dy输出为(a+xy)-1+y(dxy2+dy(xy+(a+xy)Loga+xy)例 1.4 (教材 例 1.3) 设 ,求vueyvuexcos,sin .,yvxu输入eq1=Dx=Eu+u*Sinv,x,NonConstants-u,v(*第一个方程两边对 x 求导数, 把 u,v 看成 x,y 的函数*)eq2=Dy=Eu-u*Cosv,x,NonConstants-u,v(*第二个方程两边对 x 求导数, 把 u,v 看成 x,y 的函数*)Solveeq1,eq2,Du,x,NonConstants-u,v,Dv,x,No
4、nConstants-u,v/Simplify(*解求导以后由 eq1,eq2 组成的方程组*)则输出1(,tan,t, vSinECosuvsNoCxvDiviu uu其中 Du,x,NonConstants-u,v表示 u 对 x 的偏导数, 而 Dv,x,NonCosnstants-u,v表示 v对 x 的偏导数 . 类似地可求得 u,v 对 y 的偏导数.微分学的几何应用例 1.5 求出曲面 在点(1,1)处的切平面、法线方程, 并画出图形.2xz解(1) 画出曲面的图形. 曲面的参数方程为22,0,cos/inrzruyf输入命令Clearf;fx_,y_=2x2+y2;p1=Plo
5、t3Dfx,y,x,-2,2,y,-2,2;g1=ParametricPlot3Dr*Sinu/Sqrt2.,r*Cosu,r2,u,0,2*Pi,r,0,2则输出相应图形(图 1.2).-1 01-2 -10 1201234图 1.2(2) 画出切平面的图形. 输入命令a=Dfx,y,x/.x-1,y-1;b=Dfx,y,y/.x-1,y-1;px_,y_=f1,1+a(x-1)+b(y-1);g2=Plot3Dpx,y,x,-2,2,y,-2,2;则输出切平面方程为 及相应图形(图 1.3).,012yx-2 -10 12-2-1012-100图 1.3(3) 画出法线的图形. 输入命令l
6、yx_=1+b(x-1)/a;lzx_=f1,1-(x-1)/a;g3=ParametricPlot3Dx,lyx,lzx,x,-2,2;Showp1,g2,g3,AspectRatio-Automatic,ViewPoint-2.530,-1.025,2.000;则输出相应图形(图 1.4).-2-1 012-2-1012-10010图 1.4例 1.6 (教材 例 1.4) 求曲面 在点 处的切平面方程, 并把14),(2yxk 2164,曲面和它的切平面作在同一图形里.输入Cleark,z;kx_,y_=4/(x2+y2+1);(*定义函数 k(x,y)*)kx=Dkx,y,x/.x-1
7、/4,y-1/2;(*求函数 k(x,y)对 x 的偏导数, 并代入在指定点的值*)ky=Dkx,y,y/.x-1/4,y-1/2;(*求函数 k(x,y)对 y 的偏导数, 并代入在指定的值*)z=kx*(x-1/4)+ky*(y-1/2)+k1/4,1/2;(*定义在指定点的切平面函数*)再输入qm=Plot3Dkx,y,x,-2,2,y,-2,2,PlotRange-0,4,BoxRatios-1,1,1,PlotPoints-30,DisplayFunction-Identity;qpm=Plot3Dz,x,-2,2,y,-2,2,DisplayFunction-Identity;Sh
8、owqm,qpm,DisplayFunction-$DisplayFunction则输出所求曲面与切平面的图形(图 1.5).-2 -10 12-2-101 201234图 1.5多元函数的极值例 1.7 (教材 例 1.5) 求 的极值.xyxyxf 93),(23输入Clearf;fx_,y_=x3-y3+3x2+3y2-9x;fx=Dfx,y,xfy=Dfx,y,ycritpts=Solvefx=0,fy=0则分别输出所求偏导数和驻点: 2369yxx-3,y-0,x-3,y-2,x-1,y-0,x-1,y-2再输入求二阶偏导数和定义判别式的命令fxx=Dfx,y,x,2;fyy=Dfx
9、,y,y,2;fxy=Dfx,y,x,y;disc=fxx*fyy-fxy2输出为判别式函数 的形式:2xyxff(6+6x)(6-6y)再输入data=x,y,fxx,disc,fx,y/.critpts;TableFormdata,TableHeadings-None, “x “, “y “, “fxx “, “disc “, “f “最后我们得到了四个驻点处的判别式与 的值并以表格形式列出.xfX y fxx disc f-3 0 -12 -72 27-3 2 -12 72 311 0 12 72 -51 2 12 -72 -1易见,当 时 判别式 disc=72, 函数有极大值 31;
10、,3yx,1xf当 时 判别式 disc=72, 函数有极小值-5;,f当 和 时, 判别式 disc=-72, 函数在这些点没有极值.0最后,把函数的等高线和四个极值点用图形表示出来,输入d2=x,y/.critpts;g4=ListPlotd2,PlotStyle-PointSize0.02,DisplayFunction-Identity;g5=ContourPlotfx,y,x,-5,3,y,-3,5,Contours-40,PlotPoints-60,ContourShading-False,Frame-False,Axes-Automatic,AxesOrigin-0,0,Disp
11、layFunction-Identity;Showg4,g5,DisplayFunction-$DisplayFunction则输出图 1.6.-4 -2 2-224图 1.6从上图可见, 在两个极值点附近, 函数的等高线为封闭的. 在非极值点附近, 等高线不封闭. 这也是从图形上判断极值点的方法.注:在项目一的实验 4 中,我们曾用命令 FindMinimum 来求一元函数的极值, 实际上,也可以用它求多元函数的极值, 不过输入的初值要在极值点的附近. 对本例,可以输入以下命令FindMinimumfx,y,x,-1,y,1则输出-5.,x-1.,y-2.3660310-8从中看到在 的附近
12、函数 有极小值-5, 但 y 的精度不够好 .0,1yx),(yxf例 1.8 求函数 在条件 下的极值.2xz012x输入Clearf,g,la; fx_,y_=x2+y2;gx_,y_=x2+y2+x+y-1;lax_,y_,r_=fx,y+r*gx,y;extpts=SolveDlax,y,r,x=0,Dlax,y,r,y=0,Dlax,y,r,r=0得到输出 )31(2),31(2),3(1 ,yxr再输入fx,y/.extpts/Simplify得到两个可能是条件极值的函数值 但是否真的取到条件极值呢? .,可利用等高线作图来判断.输入dian=x,y/.Tableextptss,j
13、,s,1,2,j,2,3g1=ListPlotdian,PlotStyle-PointSize0.03,DisplayFunction-Identitycp1=ContourPlotfx,y,x,-2,2,y,-2,2,Contours-20,PlotPoints-60,ContourShading-False,Frame-False,Axes-Automatic,AxesOrigin-0,0,DisplayFunction-Identity;cp2=ContourPlotgx,y,x,-2,2,y,-2,2,PlotPoints-60,Contours-0,ContourShading-Fa
14、lse,Frame-False,Axes-Automatic,ContourStyle-Dashing0.01,AxesOrigin-0,0,DisplayFunction-Identity;Showg1,cp1,cp2,AspectRatio-1,DisplayFunction-$DisplayFunction输出为 )31(2,)31(2,及图 1.7. 从图可见,在极值可疑点 ,1 ,处, 函数 的等高线与曲线 (虚线) 相切. 函数 的等高)(yxfz0(yxg ),(yxfz线是一系列同心圆, 由里向外, 函数值在增大, 在 的312),31(2x附近观察, 可以得出 取条件极大的结
15、论. 在 )(fz ),(的附近观察, 可以得出 取条件极小的结论.)31(2y ),(yxfz-2 -1 1 2-2-112图 1.7梯度场例 1.9 画出函数 的梯度向量.22),(yxzyxf解 输入命令1.0,Axes-True,AxesLabel-“x“,“y“,“z“;vecplot3d=PlotGradientField3Df,x,-1.1,1.1,y,-1.1,1.1,z,-2,2,PlotPoints-3,VectorHeads-True;Showvecplot3d, cp3d;则输出相应图形(图 1.8)图 1.8例 1.10 在同一坐标面上作出和 21),(yxyxu ,
16、1),(2yxyxv的等高线图( ), 并给出它们之间的关系.0解 输入命令Identity;uplot=ContourPlotu,x,-2,2,y,-2,2,ContourStyle-GrayLevel0,ContourShading-False,DisplayFunction-Identity,Contours-40,PlotPoints-40;g1=Showuplot,ugradplot,DisplayFunction-$DisplayFunction;vgradplot=PlotGradientFieldv,x,-2,2,y,-2,2,DisplayFunction-Identity;
17、vplot=ContourPlotv,x,-2,2,y,-2,2,ContourStyle-GrayLevel0.7,ContourShading-False,DisplayFunction-Identity,Contours-40,PlotPoints-40;g2=Showvplot,vgradplot,DisplayFunction-$DisplayFunction;g3=Showuplot,vplot,DisplayFunction-$DisplayFunction;g4=Showugradplot,vgradplot,DisplayFunction-$DisplayFunction;则
18、输出相应图形(图 1.9),其中(a) 的梯度与等高线图;),(yxu(b) 的梯度与等高线图;v(c) 与 的等高线图;,),(d) 与 的梯度图.)(-2 -1 0 1 2-2-1012-2 -1 0 1 2-2-1012(a) (b)-2 -1 0 1 2-2-1012(c) (d)图 1.9从上述图中可以看出它们的等高线为一族正交曲线. 事实上, 有 ,22 xvyxuyvxu 且 它们满足拉普拉斯方程,0vu 022yvxyx例 1.11 (教材 例 1.6) 设 作出 的图形和等高线, 再作出它的,),()(ef)(f梯度向量 gradf 的图形. 把上述等高线和梯度向量的图形叠加
19、在一起 , 观察它们之间的关系.输入调用作向量场图形的软件包命令60, Contours-25,ContourShading-False,Frame-False,Axes-Automatic,AxesOrigin-0,0td=PlotGradientFieldfx,y,x,-2,2,y,-2,2,Frame-FalseShowdgx,td输出为图 1.10. 从图可以看到 平面上过每一点的等高线和梯度向量是垂直的, 且梯度的Oxy方向是指向函数值增大的方向.-2 -1 1 2-2-112图 1.10例 1.12 求出函数 的极值, 并画出函数 的等高线、驻点4),(yxyxf),(yxf以及
20、的梯度向量的图形.),(yxf输入命令False,PlotPoints-100,Contours-4,-2,0,2,4,10,20;fieldplot=PlotGradientField-f,x,-2,2,y,-3,3,ScaleFunction-(Tanh#/5critptplot=ListPlot-Sqrt2,-2*Sqrt2,0,0,Sqrt2,2*Sqrt2,PlotStyle-PointSize0.03;Showconplot,fieldplot,critptplot;则得到 的最小值 以及函数的图形(图 1.11).),(yxf .4)823.,41.(f-2 -1 0 1 2-3
21、-2-10123图 1.11实验 2 多元函数积分学(基础实验)实验目的掌握用 Mathematica 计算二重积分与三重积分的方法; 深入理解曲线积分、曲面积分的概念和计算方法. 提高应用重积分和曲线、曲面积分解决各种问题的能力.计算重积分例 2.1 (教材 例 2.1) 计算 其中 为由 所围成的有,2dxyD ,2yxy2界区域.先作出区域 的草图, 易直接确定积分限,且应先对 积分, 因此, 输入Integratex*y2,y,1,2,x,2-y,Sqrty则输出所求二重积分的计算结果 .12093例 2.2 (教材 例 2.2) 计算 其中 为,)(2dxyeDD.12yx如果用直角
22、坐标计算, 输入Clearf,r;fx,y=Exp-(x2+y2);Integratefx,y,x,-1,1,y,-Sqrt1-x2,Sqrt1-x2则输出为dx1Erfe21x2其中 Erf 是误差函数 . 显然积分遇到了困难. 如果改用极坐标来计算, 也可用手工确定积分限. 输入Integrate(fx,y/.x-r*Cost,y-r*Sint)*r,t,0,2 Pi,r,0,1 则输出所求二重积分的计算结果e如果输入NIntegrate(fx,y/.x-r*Cost,y-r*Sint)*r,t,0,2 Pi,r,0,1 则输出积分的近似值1.98587例 2.3 (教材 例 2.3) 计
23、算 , 其中 由曲面 与dxyzx)(22yxz围成.2yxz先作出区域 的图形. 输入g1=ParametricPlot3DSqrt2*Sinfi*Costh, Sqrt2*Sinfi*Sinth, Sqrt2*Cosfi,fi,0,Pi/4,th,0,2Pig2=ParametricPlot3Dz*Cost,z*Sint,z,z,0,1,t,0,2PiShowg1,g2,ViewPoint-1.3,-2.4,1.0则分别输出三个图形(图 2.1(a), (b), (c)).-1-0.5 00.51 -1-0.500.5111.11.21.31.4(a)-1 -0.500.51 -1-0.5
24、00.5100.250.50.751(b)-1 -0.50 0.51-1 -0.50 0.5 100.51图 2.1考察上述图形, 可用手工确定积分限. 如果用直角坐标计算, 输入gx_,y_,z_=x2+y2+z;Integrategx,y,z,x,-1,1,y,-Sqrt1-x2, Sqrt1-x2,z,Sqrtx2+y2,Sqrt2-x2-y2执行后计算时间很长, 且未得到明确结果.现在改用柱面坐标和球面坐标来计算. 如果用柱坐标计算,输入Integrate(gx,y,z/.x-r*Coss,y-r*Sins)*r, r,0,1,s,0,2Pi,z,r,Sqrt2-r2则输出1528如果
25、用球面坐标计算,输入Integrate(gx,y,z/.x-r*Sinfi*Cost,y-r*Sinfi*Sint,z-r*Cosfi)*r2*Sinfi,s,0,2Pi,fi,0,Pi/4,r,0,Sqrt2则输出32165这与柱面坐标的结果相同.重积分的应用例 2.4 求由曲面 与 所围成的空间区域 的体yxf1,2,yxg积.输入Clearf,g;fx_,y_=1-x-y;gx_,y_=2-x2-y2;Plot3Dfx,y,x,-1,2,y,-1,2Plot3Dgx,y,x,-1,2,y,-1,2Show%,%一共输出三个图形, 最后一个图形是图 2.1.-1 01 2-1012-6-4
26、-202图 2.2首先观察到 的形状. 为了确定积分限, 要把两曲面的交线投影到 平面上输入Oxyjx=Solvefx,y=gx,y,y得到输出22451,4512 xyxy为了取出这两条曲线方程, 输入y1=jx1,1,2y2=jx2,1,2输出为24512x再输入tu1=Ploty1,x,-2,3,PlotStyle-Dashing0.02,DisplayFunction-Identity;tu2=Ploty2,x,-2,3,DisplayFunction-Identity;Showtu1,tu2,AspectRatio-1, DisplayFunction-$DisplayFunctio
27、n输出为图 2.2, 由此可见, 是下半圆(虚线), 是上半圆,因此投影区域是一个圆 .1y2y-0.5 0.5 1 1.5-0.50.511.5图 2.2设 的解为 与 ,则 为 的积分限. 输入21y1x21xxvals=Solvey1=y2,x输出为 6,6为了取出 , 输入21xx1=xvals1,1,2x2=xvals2,1,2输出为 612这时可以作最后的计算了. 输入Volume=Integrategx,y-fx,y,x,x1,x2,y,y1,y2/Simplify输出结果为 89例 2.5 (教材 例 2.4) 求旋转抛物面 在 平面上部的面积24yxzOx.S先调用软件包,
28、输入r*Cost,y-r*Sint;Integratez1*r,t,0,2 Pi,r,0,2/Simplify则输出所求曲面的面积176例 2.6 在 平面内有一个半径为 2 的圆, 它与 轴在原点 相切, 求它绕 轴旋Oxz zOz转一周所得旋转体体积.先作出这个旋转体的图形. 因为圆的方程是 ,42xz它绕 轴旋转所得的圆环面的方程为z, )(16)( 222yxzyx所以圆环面的球坐标方程是 输入.sin4rSphericalPlot3D4 Sint,t,0,Pi,s,0,2 Pi,PlotPoints-30,ViewPoint-4.0,0.54,2.0输出为图 2.4. -4-2024
29、-4 -2 0 2 4-2-1012图 2.4这是一个环面, 它的体积可以用三重积分计算(用球坐标). 输入Integrater2*Sint,s,0,2 Pi,t,0,Pi,r,0,4 Sint得到这个旋转体的体积为 216计算曲线积分例 2.7 (教材 例 2.5) 求 , 其中 积分路径为Ldszyxf),(,103,2yxzyxf: 32tt.0注意到,弧长微元 , 将曲线积分化为定积分,输入zdsttt22Clearx,y,z; luj=t,t2,3t2;Dluj,t则输出 对 的导数zyx,t6,21再输入ds=SqrtDluj,t.Dluj,t;Integrate(Sqrt1+30
30、 x2+10y/.x-t, y-t2,z-3t2)*ds,t,0,2则输出所求曲线积分的结果:326/3.例 2.8 (教材 例 2.6) 求 , 其中drFL. .20,sinco2)(,(356 tjttrjxyi输入vecf=x*y6,3x*(x*y5+2);vecr=2*Cost,Sint;Integrate(vecf.Dvecr,t)/.x-2Cost,y-Sint, t,0,2 Pi则输出所求积分的结果12例 2.9 求锥面 与柱面 的交线的长度.0,22zyx xy2先画出锥面和柱面的交线的图形. 输入g1=ParametricPlot3DSinu*Cosv, Sinu*Sinv
31、,Sinu, u,0,Pi,v,0,2Pi,DisplayFunction-Identity;g2=ParametricPlot3DCost2,Cost*Sint,z,t,0,2Pi,z,0,1.2, DisplayFunction-Identity;Showg1,g2,ViewPoint-1,-1,2,DisplayFunction-$DisplayFunction输出为图 2.5.-1 -0.500.51-1-0.500.5 100.51图 2.5输入直接作曲线的命令ParametricPlot3DCost2,Cost*Sint,Cost,t,-Pi/2,Pi/2, ViewPoint-1
32、,-1,2,Ticks-False输出为图 2.6.图 2.6为了用线积分计算曲线的弧长, 必须把曲线用参数方程表示出来. 因为空间曲线的投影曲线的方程为 , 它可以化成 , 再代入锥面方程xy2 tx2cos,sinty, 得22zyx.2/,costz因为空间曲线的弧长的计算公式是, 212t dtztytxs因此输入Clearx,y,z;x=Cost2;y=Cost*Sint;z=Cost;qx=x,y,z;IntegrateSqrtDqx,t. Dqx,t/Simplify,t,-Pi/2,Pi/2输出为 2Elliptice-1这是椭圆积分函数. 换算成近似值. 输入%/N输出为3.
33、8202计算曲面积分例 2.10 (教材 例 2.7) 计算曲面积分 , 其中 为锥面dSzxy)( 被柱面 所截得的有限部分.2yxzxy22注意到,面积微元 , 投影曲线 的极坐标方程为dyzdS1 xy22,costtr将曲面积分化作二重积分,并采用极坐标计算重积分.输入Clearf,g,r,t;fx_,y_,z_=x*y+y*z+z*x;gx_,y_=Sqrtx2+y2;mj=Sqrt1+Dgx,y,x2+Dgx,y,y2/Simplify;Integrate(fx,y,gx,y*mj/.x-r*Cost,y-r* Sint)*r,t,-Pi/2,Pi/2,r,0,2Cost则输出所求
34、曲面积分的计算结果15264例 2.11 计算曲面积分 其中 为球面,33dxyzydzx 的外侧. 22axy可以利用两类曲面积分的关系, 化作对曲面面积的曲面积分 . 这里ndsA. 因为球坐标的体积元素 注意到在球azyxnzyxA/,3 ,sin2drdv面 上 , 取 后得到面积元素的表示式:ar1dr.0,si2 把对面积的曲面即直接化作对 的二重积分. 输入,ClearA,fa,ds;A=x3,y3,z3;fa=x,y,z/a;ds=a2*Sinu;Integrate(A.fa/.x-a*Sinu*Cosv,y-a*Sinu*Sinv,z-a*Cosu)*ds/Simplify,
35、u,0,Pi,v,0,2Pi输出为512a如果用高斯公式计算, 则化为三重积分, 其中 为 . dvzyx22322azyx采用球坐标计算, 输入r*Sinu*Cosv,y-r*Sinu*Sinv,z-r*Cosu)*r2Sinu,v,0,2Pi,u,0,Pi,r,0,a输出结果相同.实验 3 最小二乘拟合(基础实验)实验目的 了解曲线拟合问题与最小二乘拟合原理. 学会观察给定数表的散点图, 选择恰当的曲线拟合该数表.最小二乘拟合原理给定平面上的一组点 ,2,1),(nkyx寻求一条曲线 使它较好的近似这组数据, 这就是曲线拟合. 最小二乘法是进行曲)(xfy线拟合的常用方法.最小二乘拟合的原
36、理是, 求 使),(fnkkyx12达到最小. 拟合时, 选取适当的拟合函数形式 ),()()()(10 xccxf m其中 称为拟合函数的基底函数.为使 取到极小值, 将 的表达式,)(,10xm )(xf代入, 对变量 求函数 的偏导数, 令其等于零, 就得到由 个方程组成的方程组, 从中ic 1m可解出 ).,210(mi曲线拟合例 3.1 (教材 例 3.1) 为研究某一化学反应过程中温度 对产品得率 的影响, )(Cx(%)y测得数据如下:x 100 110 120 130 140 150 160 170 180 190y 45 51 54 61 66 70 74 78 85 89试
37、求其拟合曲线.输入点的坐标, 作散点图, 即输入b2=100,45,110,51,120,54,130,61,140,66,150,70,160,74,170,78,180,85,190,89;fp=ListPlotb2则输出题设数据的散点图.120 140 160 18060708090通过观察发现散点基本位于一条直线附近, 可用直线拟合. 输入Fitb2,1,x,x (*用 Fit 作拟合, 这里是线性拟合*)则输出拟合直线-2.73939+0.48303x作图观察拟合效果. 输入gp=Plot%,x,100,190,PlotStyle-RGBColor1,0,0,DisplayFunct
38、ion-Identity; (*作拟合曲线的图形 *)Showfp,gp,DisplayFunction-$DisplayFunction(*显示数据点与拟合曲线*)则输出平面上的点与拟合抛物线的图形(图 3.1).120 140 160 18060708090图 3.1例 3.2 (教材 例 3.2) 给定平面上点的坐标如下表: 3627.105.94.708.43.6978.5.30.5124.yx试求其拟合曲线.输入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,
39、9.0253,0.9,10.3627;pd=ListPlotdata;则输出题设数据的散点图.0.4 0.6 0.8678910观察发现这些点位于一条抛物线附近. 用抛物线拟合, 即取基底函数 输入.,12xf=Fitdata,1,x,x2,x则输出5.30661-1.83196x+8.17149x2再输入fd=Plotf,x,0,1,DisplayFunction-Identity;Showpd,fd,DisplayFunction-$DisplayFunction则输出平面上的点与拟合抛物线的图形(图 3.2).0.2 0.4 0.6 0.8 167891011图 3.2下面的例子说明 F
40、it 的第二个参数中可以使用复杂的函数, 而不限于 等.2,x例 3.3 (教材 例 3.3) 使用初等函数的组合进行拟合的例子.先计算一个数表. 输入ft=TableN1+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/(,sin1xFitft,1,Sinx,Exp-x/3,Exp-x,x则输出拟合函数 102450489. 63/15 xSinx其中有些基函数的系数非常小, 可将它们删除. 输入Ch
41、op%则输出 3/.21xe实际上,我们正是用这个函数做的数表. 注:命令 Chop 的基本格式为Chopexpr, 其含义是去掉表达式 expr 的系数中绝对值小于 的项, 的默认值为 .10实验 4 水箱的流量问题(综合实验)实验目的 掌握应用最小二乘拟合原理分析和解决实际问题的思想和方法,能通过观察测试数据的散点图,建立恰当的数学模型,并用所学知识分析和解决所给问题.问题 (1991 年美国大学生数学建模竞赛的 A 题. 问题中使用的长度单位为 E(英尺, 1 E=30.24cm), 容积单位是 G(加仑, 1 G=3.785L).某些州的用水管理机构需估计公众的用水速度(单位:G/h)
42、和每天的总用水量. 许多供水单位由于没有测量流入或流出量的设备, 而只能测量水箱中的水位( 误差不超过 5%). 当水箱水位低于水位 L 时, 水泵开始工作将水灌入水箱, 直至水位达到最高水位 H 为止. 但是依然无法测量水泵灌水流量, 因此, 在水泵工作时无法立即将水箱中的水位和水量联系起来. 水泵一天灌水 12 次, 每次约 2h. 试估计在任一时刻(包括水泵灌水期间) t 流出水箱的流量并估计一天的总用水量.)(tf表 1 给出了某镇某一天的真实用水数据. 水箱是直径为 57E, 高为 40E 的正圆柱体. 当水位落到 27E 以下, 水泵自动启动把水灌入水箱; 当水位回升至 35.5E
43、 时, 水泵停止工作.表 1 时间/s 水位 E20时间/s 水位 E21003316663510619139371792121240252232854332284359323175311030542994294728922850279527522697泵水4663649953539365725460574645546853571854750217925482649335032603167308730122927284227672697泵水泵水393323943543318泵水35503445859688995393270347533973340模型假设(1) 影响水箱流量的唯一因素是该区公众
44、对水的普通需求. 所给数据反映该镇在通常情况下一天的用水量, 不包括任何非常情况, 如水泵故障、水管破裂、自然灾害等 . 并且认为水位高度、大气情况、温度变化等物理因素对水的流速均无直接影响;(2) 水泵的灌水速度为常数;(3) 从水箱中流出水的最大流速小于水泵的灌水速度. 为了满足公众的用水需求不让水箱中的水用尽, 这是显然的要求;(4) 因为公众对水的消耗量是以全天的活动(诸如洗澡、做饭、洗衣服等) 为基础的, 所以,可以认为每天的用水量分布都是相似的;(5) 水箱的水流量速度可用光滑曲线来近似.问题分析与模型建立为方便起见,记 V 表示水的容积; 表示时刻 (单位:h)水的容积; 表示流出水箱的iVit )(tf水的流速( 单位;G/h), 它是时间的函数;p 表示水泵的灌水速度(G/h).先将表 1 中数据作变换, 时间单位用小时(h), 水位高转换成水的体积( 单位:,2hrV). 输入G48.7E,03tt=0,3316,6635,10619,13937,17921,21240,25223, 28543,32284,35932,39332,39435,43