9_基于Mathematica的数值计算.ppt

第九讲数值计算,6.1求近似函数在生产和实验中,人们经常遇到需要通过某个未知的函数f(x)在有限个给定点的函数值:xi,yi,i=1,2,.,n,这里f(xi)=yi去获得函数f(x)的近似函数(x),求近似函数(x)的方法主要有拟合方法和插值方法。6.1.1曲线拟合曲线拟合主要用来求一元近似函数,它是根据最小二乘原理的意义下获得近似函数的,此近似函数具有在数据点处的误差平方和最小的特点。记函数集合:M=Span0,1,2,m=(x)|(x)=a00(x)+a11(x)+amm(x),aiR称集合M为函数0,1,2,m张成的空间，m+1个函数0(x),1(x),2(x),m(x)称为拟合基函数集合,它们都是已知的函数。,Mathematica曲线拟合的一般形式为:Fit数据点集合,拟合基函数集合,自变量名具体的拟合命令有:命令形式1：Fitx1,y1,x2,y2,.,xn,yn,0,1,2,m,x功能：根据数据点集x1,y1,x2,y2,.,xn,yn求出具有拟合函数为(x)=a00(x)+a11(x)+amm(x)形式的近似函数(x)命令形式2：Fity1,y2,.,yn,0,1,2,m,x功能：根据数据点集1,y1,2,y2,.,n,yn求出具有拟合函数为(x)=a00(x)+a11(x)+amm(x)形式的近似函数(x)命令形式3：Fitx1,y1,x2,y2,.,xn,yn,Tablexi,i,0,m,x功能：根据数据点集x1,y1,x2,y2,.,xn,yn求出拟合函数为m次多项式的近似函数(x)=a0+a1x+a2x2+amxm,多项式拟合算法,输入n+1个拟合点：(xi,yi),i=0,1,n根据散点图确定拟合多项式的次数m计算相应正规线性方程组的系数和右端项解正规正规线性方程组，得解：a0*,a1*,am*写出拟合多项式*(x)=a0*+a1*x+a2*x2+am*xm,求m次多项式拟合程序,Clearxi,xx,yi;xi=Inputxi=yi=Inputyi=n=Lengthxi;h=ListPlotTablexii,yii,i,1,n,PlotStyle-PointSize0.04m=Input多项式次数m=s=TableSumxiki,k,1,n,i,0,2m;a=Tablesi+j-1,i,1,m+1,j,1,m+1;Printa=,MatrixForma;b=TableSumxiki*yik,k,1,n,i,0,m;Printb=,MatrixFormb;xx=Tablexi,i,1,m+1;g=Solvea.xx=b,xx;fa=Sumxi*t(i-1),i,1,m+1/.g1;p=fa/Np1=Plotp,t,xi1,xin,DisplayFunction-Identity;Showp1,h,DisplayFunction-$DisplayFunction;,说明：本程序用于求m次多项式拟合。程序执行后，按要求通过键盘输入拟合基点xi:x0,x1,.,xn、对应函数值yi:y0,y1,yn后，计算机给出散点图和请求输入拟合多项式次数的窗口，操作者可以根据散点图确定拟合多项式的次数通过键盘输入，程序即可给出对应的正规方程组系数矩阵a、常数项b、m次拟合多项式和由拟合函数图形和散点图画在一起的图形。,程序中变量说明xi：存放拟合基点x0,x1,.,xnyi:存放对应函数值y0,y1,ynm:存放拟合多项式次数a:存放正规方程组系数矩阵b:存放正规方程组常数项p:存放m次拟合多项式h:存放散点图p1：存放拟合函数图形xx:定义正规方程组变量，存放m次拟合多项式的系数注：语句s=TableSumxiki,k,1,n,i,0,2m、a=Tablesi+j-1,i,1,m+1,j,1,m+1、b=TableSumxiki*yik,k,1,n,i,0,m是用简化的正规方程组编程的。,例1已知一组实验数据x1345678910f(x)1054211234用多项式拟合求其拟合曲线。,解：执行m次多项式拟合程序后，在输入的两个窗口中按提示分别输入1,3,4,5,6,7,8,9,10，10,5,4,2,1,1,2,3,4每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上画出散点图。,由于该散点图具有2次多项式形状，因此在确定选择多项式次数窗口输入2，按OK”按扭后得如下输出结果。a=953381533813017381301725317b=32147102513.4597-3.60531t+0.267571t2,于是得求出的二次多项式拟合函数为(t)=13.4597-3.60531t+0.267571t2而且从图形上看拟合效果很好。,直接利用Fit命令：g=1,10,3,5,4,4,5,2,6,1,7,1,8,2,9,3,10,4;Hh=Fitg,Tablexk,k,0,6,x;Tu1=PlotHh,x,1,10;Tu2=ListPlotg;ShowTu1,Tu2,例2已知一组实验数据x681012141618202224f(x)4.64.84.64.95.05.45.15.55.66.0修改多项式拟合程序使其具有可以选择不同多项式拟合函数功能，并用此程序确定本题多项式拟合曲线。解：在拟合多项式程序中加入评价拟合效果的判定人机交互语句tu=InputSatisfatory?0(No)or1Yes和While语句来调控何时终止实验，调控变量用tu取值确定，取值0表示不满意，1表示满意。此外，去掉正规方程组系数及拟合多项式的输出，代之以在图形上表出拟合多项式的次数提示，修改后的程序如下。,xi=Tablei,i,6,24,2;yi=4.6,4.8,4.6,4.9,5,5.4,5.1,5.5,5.6,6;n=Lengthxi;h=ListPlotTablexii,yii,i,1,n,PlotStyle-PointSize0.04tu=0;Whiletu=0,m=Input多项式次数m;s=TableSumxiki,k,1,n,i,0,2m;a=Tablesi+j-1,i,1,m+1,j,1,m+1;(*Printa=,MatrixForma;*)b=TableSumxiki*yik,k,1,n,i,0,m;(*Printb=,MatrixFormb;*)xx=Tablexi,i,1,m+1;g=Solvea.xx=b,xx;fa=Sumxi*t(i-1),i,1,m+1/.g1;p=fa/N;p1=Plotp,t,xi1,xin,PlotLabel-m拟合多项式,DisplayFunction-Identity;Showp1,h,DisplayFunction-$DisplayFunction;tu=InputSatisfatory?0(No)or1Yes执行该程序后，屏幕上出现拟合多项式次数输入窗口和散点图。,由于该散点图不好确定拟合次数，先用3次拟合多项式计算，因此输入“3”并用鼠标点击窗口的“OK”按扭，得如下输出图形。,屏幕出现提示是否满意的输入窗口，因为不太满意，输入：0，单击“OK”按扭并在屏幕上出现拟合多项式次数输入窗口中输入：6，单击OK，屏幕出现下一个组合图形。,继续做实验，得到如下若干图形。由于总共有10个数据点，所以拟合多项式最高次数只能到9次，因此实验到9次拟合多项式后，在屏幕出现提示是否满意的输入窗口中，输入：1，单击“OK”按扭退出实验。,直接利用Fit命令：g=6,4.6,8,4.8,10,4.6,12,4.9,14,5,16,5.4,18,5.1,20,5.5,22,5.6,24,6;Hh=Fitg,Tablexk,k,0,9,x;Tu1=PlotHh,x,6,24;Tu2=ListPlotg;ShowTu1,Tu2,线性模型拟合,线性模型拟合算法1.输入n+1个拟合点：(xi,yi),i=0,1,n2.根据散点图确定拟合基函数组3.计算相应正规线性方程组的系数和右端项4.解正规正规线性方程组，得解：a0*,a1*,am*5.写出线性拟合模型*(x)=a0*0(x)+a1*1(x)+a2*2(x)+am*m(x),求线性模型拟合程序,Clearx,xi,xx,yi;xi=Inputxi=yi=Inputyi=n=Lengthxi;h=ListPlotTablexii,yii,i,1,n,PlotStyle-PointSize0.04m1=Input输入拟合基函数组m=Lengthm1p=Tablem1/.x-xik,k,1,na=TableSumpk,i*pk,j,k,1,n,i,1,m,j,1,m/N(*Printa=,MatrixForma;*)b=TableSumpk,i*yik,k,1,n,i,1,m/N(*Printb=,MatrixFormb;*)xx=Tablexti,i,1,mg=Solvea.xx=b,xxfa=Sumxti*m1i,i,1,m/.g1pp=fa/N;p1=Plotpp,x,xi1,xin,DisplayFunction-Identity;Showp1,h,DisplayFunction-$DisplayFunction,说明：本程序用于求线性模型拟合。程序执行后，按要求通过键盘输入插值基点xi:x0,x1,.,xn、对应函数值yi:y0,y1,yn后，计算机给出散点图和请求输入拟合拟合基函数组0(x)，1(x)，2(x)、m(x)的窗口，操作者可以根据散点图确定拟合基函数组通过键盘输入，程序即可给出对应的正规方程组系数矩阵a、常数项b、线性模型拟合函数和由拟合函数图形与散点图画在一起的图形。,程序中变量说明：xi：存放拟合基点x0,x1,.,xnyi:存放对应函数值y0,y1,ynm1:存放拟合基函数组0(x)，1(x)，2(x)、m(x)m:存放拟合基函数组个数a:存放正规方程组系数矩阵b:存放正规方程组常数项p:存放拟合基函数组在拟合基点x0,x1,.,xn的函数值pp:存放求出的线性模型拟合函数h:存放散点图p1：存放拟合函数图形xx:定义正规方程组变量，存放线性模型拟合系数注：（1）语句m1=Input输入拟合基函数组产生的输入应用函数表，即用“0x，1x，mx”的形式输入。（2）Mathematica数学软件中也有一个求线性模型拟合的命令，命令形式为Fitx1,y1,x2,y2,.,xn,yn,0,1,2,m,x它可以根据数据点集x1,y1,x2,y2,.,xn,yn求出具有拟合函数为(x)=a00(x)+a11(x)+amm(x)形式的近似函数(x)。,例3已知函数在自变量x=1,2,10上数据为2.89229,2.86323,0.473147,-2.25209,-2.87003,-0.835768,1.97187,2.96841,1.23648,-1.63202，试用合适的函数进行拟合。解：执行线性模型拟合程序后，在输入的两个窗口中按提示分别输入1,2,3,4,5,6,7,8,9,10、2.89229,2.86323,0.473147,-2.25209,-2.87003,-0.835768,1.97187,2.96841,1.23648,-1.63202每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上画出散点图。,由于该散点图具有类似正弦曲线的形状，因此在确定选择拟合基函数窗口输入“1,Sinx”，按“OK”按扭后得如下输出结果。a=10.1.411191.411195.00143b=4.8155215.42390.0482789+3.07027Sinx,于是，我们得到了很好的拟合函数(x)=0.0482789+3.07027sin(x)。对于本题，如果看到散点图具有一个弯曲而选用三次多项式拟合，则输入“1,x,x2,x3”后会得到如下输出。,a=10.55.385.3025.55.385.3025.25333.385.3025.25333.220825.3025.25333.220825.1.97841106b=4.8155214.0236104.284820.71110.4765-7.37686x+1.41989x2-0.0796299x3,显然这个拟合很不满意。,直接利用Fit命令：g=2.89229,2.86323,0.473147,-2.25209,-2.87003,-0.835768,1.97187,2.96841,1.23648,-1.63202;Hh=Fitg,1,Sinx,x;Tu1=PlotHh,x,1,10;Tu2=ListPlotg;ShowTu1,Tu2,6.1.2函数插值多项式插值多项式插值是在给定n个数据点:xi,yi,i=1,2,.,n后,求出一个次数不超过n-1的多项式(x)作为函数y=f(x)的近似表达式，它就是计算方法中常说的拉格朗日(Lagrange)插值或Newton插值多项式。Mathematica多项式插值的一般形式为:InterpolatingPolynomial数据点集合,自变量名具体的多项式插值命令有:命令形式1：InterpolatingPolynomialx1,y1,x2,y2,.,xn,yn,x功能：根据数据点集x1,y1,x2,y2,.,xn,yn求出n-1次插值多项式(x)命令形式2：InterpolatingPolynomialy1,y2,.,yn,x功能：根据数据点集1,y1,2,y2,.,n,yn求出n-1次插值多项式(x),例给定数据表x0123y=f(x)13512用Lagrange插值法求三次插值多项式,并给出函数f(x)在x=1.4的近似值。,Mathematica命令：Hx=InterpolatingPolynomial0,1,1,3,2,5,3,12,x;tu1=ListPlot0,1,1,3,2,5,3,12;tu2=PlotHx,x,0,3;Showtu1,tu2Hx/.x1.43.52,例2给定数据表x4.00024.01044.02334.0294y0.60208170.60318770.60458240.6052404(1)用插值法求三次插值多项式(2)求f(4.011)的近似值,Mathematica命令：Hx=InterpolatingPolynomial4.0002,0.602082,4.0104,0.603188,4.0233,0.604582,4.0294,0.60524,x;tu1=ListPlot4.0002,0.602082,4.0104,0.603188,4.0233,0.604582,4.0294,0.60524;tu2=PlotHx,x,4,4.03;Showtu1,tu2Hx/.x4.0110.603253,例3多项式插值的误差估计式中可以看到，当插值节点越多时误差会越小，这个结论正确吗？通过实验说明该结论的正确性。解:考虑函数f(x)=1/(1+x2)在区间-4,4内选取不同个数的等距插值节点做观察，这里分别选-4,4内的9个和11个的等距节点来做实验，将对应的插值函数图与被插函数f(x)=1/(1+x2)画在一起做观察，为简单起见，这里用Mathematica命令做实验，对应命令为,In1:=u=Tablex,(1+x2)(-1),x,-4,4;（*采取f(x)在-4,4内的9个插值点In2:=g=ListPlotu,PlotStyle-PointSize0.04（*将散点图图形文件存放在变量g中In3:=s=InterpolatingPolynomialu,x;（*将插值函数存放在变量s中In4:=t=Plots,(1+x2)(-1),x,-4,4,PlotStyle-Thickness0.005,Thickness0.006（*将插值函数s与f(x)画在一起的图形文件存放在变量t中In5:=Showt,g（*将散点图,插值函数s,f(x)画在一个坐标系中,在-4,4中选9个等距节点的插值函数与被插函数图，粗线为被插函数图,In6:=u=Tablex,(1+x2)(-1),x,-4,4,0.8;（*采取f(x)在-4,4内的11插值点In7:=g=ListPlotu,PlotStyle-PointSize0.04In8:=s=InterpolatingPolynomialu,x;In9:=t=Plots,(1+x2)(-1),x,-4,4,PlotStyle-Thickness0.005,Thickness0.006In10:=Showt,g,在-4,4中选11个等距节点的插值函数与被插函数图,分段多项式插值,Mathematica分段多项式插值的一般形式为:Interpolation数据点集合具体的分段多项式插值命令有:命令形式1：Interpolationx1,y1,x2,y2,.,xn,yn功能：根据数据点集x1,y1,x2,y2,.,xn,yn求出分段插值多项式(x)命令形式2：Interpolationy1,y2,.,yn功能：根据数据点集1,y1,2,y2,.,n,yn求出分段插值多项式(x),例5:用分段多项式插值重做例4的插值问题,并验证,随着插值点的增多,分段插值函数的误差会越来越小。(见图)解:In27:=r=InterpolationTablex,(1+x2)(-1),x,-4,4Out27=InterpolatingFunction-4,4,In28:=(1+x2)-1/.x-3.7Out28=0.0680735In29:=r3.7Out29=0.0734In30:=d=Tablex,(1+x2)(-1),x,-4,4,0.5;In31:=r1=InterpolationdOut31=InterpolatingFunction-4,4.,In32:=r13.7Out32=0.0681761In33:=Plotr1x,(1+x2)(-1),x,-4,4,下一页,返回,分段线性插值,Interpolationx1,y1,x2,y2,.,xn,yn，InterpolationOrder1用如上Mathematica命令求出的分段线性插值函数没有给出具体的分段函数表达式,而是用“InterpolatingFunctionx1,xn,”作为所求的分段插值函数,通常可以用变量=Interpolation数据点集合把所得的分段插值函数存放在变量中,如果要计算插值函数在某一点的如p点的值,只要输入“变量p”即可，如果要表示这个插值函数应该用“变量x”。此外，命令Interpolationx1,y1,x2,y2,.,xn,yn，InterpolationOrder3可以给出更好的分段插值函数。,例4给定数据表xi1234yi0-5-63（1）用分段线性插值函数求函数f(x)在x=1.4的近似值（2）用Mathematica命令画出分段线性插值图形,解：d=Interpolation1,0,2,-5,3,-6,4,3,InterpolationOrder-1Plotdx,x,1,4,得输出如下。InterpolatingFunction1,4,6.2解常微分方程,自变量、未知函数以及未知函数的导数（或微分）之间的一个（或一组）关系式称为微分方程（或微分方程组）。求出常微分方程（或微分方程组）的未知函数（或未知函数组），称为解常微分方程。含有任意常数的解称为通解，否则称为特解。n阶常微分方程的一般形式为：利用Mathematica可以象高等数学一样求出常微分方程的解析解（准确解），如果求不出解析解Mathematica可以求出微分方程的数值解（即近似解。Mathematica解常微分方程的命令有求常微分方程（组）的解析解命令和求常微分方程（组）的数值解命令。,6.2.1求常微分方程（组）的解析解,命令形式1：DSolveeqn,yx,x功能：求出常微分方程eqn的未知函数yx的解析通解。命令形式2：DSolveeqn1,eqn2,.,y1x,y2x,.,x功能：求出常微分方程组eqn1,eqn2,.的所有未知函数y1x,y2x,.的解析通解。注意：每个常微分方程的中的等号输入时应该用两个等号代替未知函数及未知函数的导数要用“自变量名”指示未知函数的自变量常微分方程组和未知函数组应该用大括号括起来命令形式2可以用来求常微分方程的特解eqn表示常微分方程,x是自变量名,yx表示未知函数,自变量名和未知函数可以是其他的变量名,例6:求常微分方程y=2ax的通解，a为常数解:Mathematica命令:In34:=DSolveyx=2ax,yx,xOut34=yx-ax2+C1即得本问题的通解例7:求常微分方程y+y=0的通解解:Mathematica命令:In35:=DSolveyx+yx=0,yx,xOut35=yx-C2Cosx-C1Sinx即得本问题的通解:C1，C2是任意常数。,例8:求常微分方程的特解。解:Mathematica命令:In36:=DSolveyx=x/yx+yx/x,y1=2,yx,xOut36=yx-Sqrtx2(4+2Logx)本问题的特解为:下面的操作可以检验所求特解的正确性:In37:=yx_:=Sqrtx2*(4+2Logx)In38:=y1Out38=2In39:=Dyx,x-x/yx-yx/x;In40:=Simplify%Out40=0,例9:求常微分方程组y(t)=x(t),x(t)=y(t)的通解。解:Mathematica命令:DSolveyt=xt,xt=yt,xt,yt,t6.2.2求常微分方程（组）的数值解命令命令形式1：NDSolveeqns,y,x,xmin,xmax功能：求出自变量范围为xmin,xmax且满足给定常微分方程及初值条件eqns的未知函数y的数值解命令形式2：NDSolveeqns,y1,y2,.,x,xmin,xmax功能：求出自变量范围为xmin,xmax且满足给定常微分方程及初值条件eqns的未知函数y1,y2,.的数值解。,例10:求微分方程初值问题y=x+y,y(0)=1在区间0,0.5内的数值解解:Mathematica命令:In42:=NDSolveyx=x+yx,y0=1,y,x,0,0.5Out42=y-InterpolatingFunction0.,0.5,上面显示的是所求数值解的替换形式，为得到本问题数值解，再键入:In43:=f=y/.%1Out43=InterpolatingFunction0.,0.5,In44:=Plotfx,x,0,0.5In45:=Tablefx,x,0,0.5,0.1Out45=1.,1.11034,1.24281,1.39972,1.58365,1.79744,例11:已知常微分方程组求函数x(t)和y(t)在0,10上的数值解。解:Mathematica命令:In46:=q=NDSolvext=-xt2-yt,yt=2xt-yt,x0=y0=1,x,y,t,0,10Out46=x-InterpolatingFunction0.,10.,y-InterpolatingFunction0.,10.,In47:=f=x/.Firstq;g=y/.Lastq;In48=f0.4Out48=0.375078In49=g4Out49=-0.0912007In50=ParametricPlotft,gt,t,0,10,Euler方法,Euler方法算法输入函数f(x,y)、初值y0、自变量区间端点a,b步长h计算节点数n和节点xk用Euler公式yn+1=yn+hf(xn,yn)求数值解,Euler方法程序Clearx,y,ffx_,y_=Input函数f(x,y)=y0=Input初值y0=a=Input左端点a=b=Input右端点b=h=Input步长h=n=(b-a)/h;Fori=0,iInterpolatingFunctionrange,的形式给出的,其中的InterpolatingFunctionrange,是所求的插值函数表示的数值解,range就是所求数值解的自变量范围。,用Euler方法求初值问题的数值解。取步长h=0.1,并在一个坐标系中画出数值解与准确解的图形。解：执行Euler方法程序后，在输入的窗口中按提示分别输入y-2x/y、1、0、1、0.1，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果。,y(0.1)=1.1y(0.2)=1.19182y(0.3)=1.27744y(0.4)=1.35821y(0.5)=1.43513y(0.6)=1.50897y(0.7)=1.58034y(0.8)=1.64978y(0.9)=1.71778y(1.)=1.78477,在一个坐标系中画出数值解与准确解的图形如下。图中点图是数值解，曲线为准确解。从图中可以知道数值解与准确解比较接近，但还是有误差的。,改进的Euler方法算法,1.输入函数f(x,y)、初值y0、自变量区间端点