第六章数值计算命令与例题

1、 第六章 数值计算命令与例题数值计算命令与例题l北京交通大学6.1 求近似函数求近似函数l在生产和实验中在生产和实验中, 人们经常遇到需要通过某个未知的函数人们经常遇到需要通过某个未知的函数f(x)在有限个给定点的函数值在有限个给定点的函数值:xi, yi, i=1,2,., n, 这里这里 f(xi) = yi 去获得函数去获得函数f(x)的近似函数的近似函数 (x), 求近似函数求近似函数 (x)的方法的方法主要有主要有拟合拟合方法和方法和插值插值方法。方法。6.1.1 曲线拟合曲线拟合l 曲线拟合主要用来求一元近似函数, 它是根据最小二乘原理的意义下获得近似函数的, 此近似函数具有在数据

2、点处的误差平方和最小的特点。记函数集合:M=Span 0, 1, 2, m= (x)| (x)= a0 0(x)+a1 1(x)+am m(x), ai Rl称集合称集合M为函数为函数 0, 1, 2, m张成的空间，张成的空间，m+1个函数个函数 0(x), 1(x), 2(x), m(x)称为拟合基函数集合称为拟合基函数集合, 它们都是它们都是已知的函数。已知的函数。lMathematica曲线拟合的一般形式为曲线拟合的一般形式为: Fit数据点集合数据点集合, 拟合基函数集合拟合基函数集合, 自变量名自变量名具体的拟合命令有具体的拟合命令有: 命令形式命令形式1：Fitx1,y1,x2,

3、y2,.,xn,yn, 0, 1, 2, m ,x功能：根据数据点集功能：根据数据点集x1,y1,x2,y2,.,xn,yn求出具有拟合函数为求出具有拟合函数为 (x)= a 0 0(x)+a1 1(x)+a m m(x) 形式的近似函数形式的近似函数 (x) 命令形式命令形式2：Fity1,y2,.,yn, 0, 1, 2, m ,x功能：根据数据点集功能：根据数据点集1,y1,2,y2,.,n,yn求出具有拟合函数为求出具有拟合函数为 (x)= a 0 0(x)+a1 1(x)+a m m(x)形式的近似函数形式的近似函数 (x)命令形式命令形式3：Fitx1,y1,x2,y2,.,xn,

4、 yn, Tablexi,i,0,m ,x功能：根据数据点集功能：根据数据点集x1,y1,x2,y2,.,xn,yn求出拟合函数为求出拟合函数为m次多次多项式的近似函数项式的近似函数 (x) =a 0 +a1x+ a2x2 +a mx m多项式拟合算法多项式拟合算法 l输入n+1个拟合点： (xi, yi),i=0,1,nl根据散点图确定拟合多项式的次数ml计算相应正规线性方程组的系数和右端项l解正规正规线性方程组，得解：a0*,a1*,a m*l写出拟合多项式*(x)= a0*+ a1*x+ a2*x2+ + am*xm求求m次多项式拟合程序次多项式拟合程序lClearxi,xx,yi;lx

5、i=Inputxi=lyi=Inputyi=ln=Lengthxi;lh=ListPlotTablexii,yii,i,1,n,PlotStyle-PointSize0.04lm=Input多项式次数m=ls=TableSumxiki,k,1,n,i,0,2m;la=Tablesi+j-1,i,1,m+1,j,1,m+1;lPrinta=,MatrixForma;lb=TableSumxiki*yik,k,1,n,i,0,m;lPrintb=,MatrixFormb;lxx=Tablexi,i,1,m+1;lg=Solvea.xx=b,xx;lfa=Sumxi*t(i-1),i,1,m+1/.

6、g1;lp=fa/Nlp1=Plotp,t,xi1,xin,DisplayFunction-Identity;lShowp1,h,DisplayFunction-$DisplayFunction;l说明：说明：本程序用于求m次多项式拟合。程序执行后，按要求通过键盘输入拟合基点xi:x0 , x1, . , xn 、对应函数值yi: y0 , y1 , , yn 后，计算机给出散点图和请求输入拟合多项式次数的窗口，操作者可以根据散点图确定拟合多项式的次数通过键盘输入，程序即可给出对应的正规方程组系数矩阵a、常数项b、m次拟合多项式和由拟合函数图形和散点图画在一起的图形。l程序中变量说明程序中变量

7、说明lxi：存放拟合基点x0 , x1, . , xn lyi: 存放对应函数值y0 , y1 , , ynlm: 存放拟合多项式次数la: 存放正规方程组系数矩阵lb: 存放正规方程组常数项lp: 存放m次拟合多项式lh: 存放散点图lp1：存放拟合函数图形lxx:定义正规方程组变量，存放m次拟合多项式的系数l注：注：语句s=TableSumxiki,k,1,n,i,0,2m、a=Tablesi+j-1,i,1,m+1,j,1,m+1、lb=TableSumxiki*yik,k,1,n,i,0,m是用简化的正规方程组编程的。例例1已知一组实验数据x 1 3 4 5 6 7 8 9 10f(x

8、) 10 5 4 2 1 1 2 3 4用多项式拟合求其拟合曲线。l解：解：执行m次多项式拟合程序后，在输入的两个窗口中按提示分别输入1,3,4,5,6,7,8,9,10，10,5,4,2,1,1,2,3,4l每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上画出散点图。l由于该散点图具有2次多项式形状，因此在确定选择多项式次数窗口输入2，按OK”按扭后得如下输出结果。la= 9 53 381 53 381 3017 381 3017 25317lb=32147102513.4597 - 3.60531 t + 0.267571 t2l于是得求出的二次多项式拟合函数为(t)=13.4597

9、 - 3.60531 t + 0.267571 t2而且从图形上看拟合效果很好。l例例 2已知一组实验数据lx 6 8 10 12 14 16 18 20 22 24lf(x) 4.6 4.8 4.6 4.9 5.0 5.4 5.1 5.5 5.6 6.0l修改多项式拟合程序使其具有可以选择不同多项式拟合函数功能，并用此程序确定本题多项式拟合曲线。l解：解： 在拟合多项式程序中加入评价拟合效果的判定人机交互语句 tu=InputSatisfatory?0(No)or1Yes 和While语句来调控何时终止实验，调控变量用tu取值确定，取值0表示不满意，1表示满意。此外，去掉正规方程组系数及拟合

10、多项式的输出，代之以在图形上表出拟合多项式的次数提示，修改后的程序如下。lxi=Tablei,i,6,24,2;lyi=4.6,4.8,4.6,4.9,5,5.4,5.1,5.5,5.6,6;ln=Lengthxi;lh=ListPlotTablexii,yii,i,1,n,PlotStyle-PointSize0.04ltu=0;lWhiletu=0,lm=Input多项式次数m;ls=TableSumxiki,k,1,n,i,0,2m;la=Tablesi+j-1,i,1,m+1,j,1,m+1;l(*Printa=,MatrixForma;*)lb=TableSumxiki*yik,k,

11、1,n,i,0,m;l(*Printb=,MatrixFormb;*)lxx=Tablexi,i,1,m+1;lg=Solvea.xx=b,xx;lfa=Sumxi*t(i-1),i,1,m+1/.g1;lp=fa/N;lp1=Plotp,t,xi1,xin,PlotLabel-m拟合多项式,lDisplayFunction-Identity;lShowp1,h,DisplayFunction-$DisplayFunction;ltu=InputSatisfatory?0(No)or1Yesl执行该程序后，屏幕上出现拟合多项式次数输入窗口和散点图。l由于该散点图不好确定拟合次数，先用3次拟合多

12、项式计算，因此输入“3”并用鼠标点击窗口的“OK”按扭，得如下输出图形。l屏幕出现提示是否满意的输入窗口，因为不太满意，输入：0，单击“OK”按扭并在屏幕上出现拟合多项式次数输入窗口中输入：6，单击OK，屏幕出现下一个组合图形。l继续做实验，得到如下若干图形。l由于总共有10个数据点，所以拟合多项式最高次数只能到9次，因此实验到9次拟合多项式后，在屏幕出现提示是否满意的输入窗口中，输入：1，单击“OK”按扭退出实验。线性模型拟合线性模型拟合 l线性模型拟合算法线性模型拟合算法l l1.输入n+1个拟合点： (xi, yi),i=0,1,nl2.根据散点图确定拟合基函数组l3.计算相应正规线性方

13、程组的系数和右端项l4.解正规正规线性方程组，得解：a0*,a1*,a m*l5.写出线性拟合模型*(x)= a0*0(x)+ a1*1(x)+ a2*2(x)+ + am*m(x)求线性模型拟合程序求线性模型拟合程序lClearx,xi,xx,yi;lxi=Inputxi=lyi=Inputyi=ln=Lengthxi;lh=ListPlotTablexii,yii,i,1,n,PlotStyle-PointSize0.04lm1=Input输入拟合基函数组lm=Lengthm1lp=Tablem1/.x-xik,k,1,nla=TableSumpk,i*pk,j,k,1,n,i,1,m,j

14、,1,m/Nl(*Printa=,MatrixForma;*)lb=TableSumpk,i*yik,k,1,n,i,1,m/Nl(*Printb=,MatrixFormb;*)lxx=Tablexti,i,1,mlg=Solvea.xx=b,xxlfa=Sumxti*m1i,i,1,m/.g1lpp=fa/N;lp1=Plotpp,x,xi1,xin,DisplayFunction-Identity;lShowp1,h,DisplayFunction-$DisplayFunctionl说明：说明：本程序用于求线性模型拟合。程序执行后，按要求通过键盘输入插值基点xi:x0 , x1, . ,

15、xn 、对应函数值yi: y0 , y1 , , yn 后，计算机给出散点图和请求输入拟合拟合基函数组0(x)，1(x)，2(x)、m(x)的窗口，操作者可以根据散点图确定拟合基函数组通过键盘输入，程序即可给出对应的正规方程组系数矩阵a、常数项b、线性模型拟合函数和由拟合函数图形与散点图画在一起的图形。l程序中变量说明：程序中变量说明：lxi：存放拟合基点x0 , x1, . , xn lyi: 存放对应函数值y0 , y1 , , ynlm1: 存放拟合基函数组0(x)，1(x)，2(x)、m(x)lm: 存放拟合基函数组个数la: 存放正规方程组系数矩阵lb: 存放正规方程组常数项lp:

16、存放拟合基函数组在拟合基点x0 , x1, . , xn 的函数值lpp: 存放求出的线性模型拟合函数lh: 存放散点图lp1：存放拟合函数图形lxx:定义正规方程组变量，存放线性模型拟合系数l注：（注：（1）语句m1=Input输入拟合基函数组产生的输入应用函数表，即用l“0 x，1x，mx”的形式输入。l（2）Mathematica数学软件中也有一个求线性模型拟合的命令，命令形式为l Fitx1,y1,x2,y2,.,xn,yn, 0, 1, 2, m ,xl它可以根据数据点集x1,y1,x2,y2,.,xn,yn求出具有拟合函数为l (x)= a 0 0(x)+a11(x)+a mm(x

17、)l形式的近似函数(x)。l例例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，试用合适的函数进行拟合。l解：解：执行线性模型拟合程序后，在输入的两个窗口中按提示分别输入l1,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.63202l每次输

18、入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上画出散点图。l由于该散点图具有类似正弦曲线的形状，因此在确定选择拟合基函数窗口输入“1,Sinx”，按“OK”按扭后得如下输出结果。la=10. 1.41119l 1.41119 5.00143lb=4.81552l 15.4239l0.0482789 + 3.07027 Sinxl于是，我们得到了很好的拟合函数(x)=0.0482789 + 3.07027 sin(x) 。l对于本题，如果看到散点图具有一个弯曲而选用三次多项式拟合，则输入“1,x,x2,x3”后会得到如下输出。la=l 10. 55. 385. 3025.l 55. 385.

19、3025. 25333.l 385. 3025. 25333. 220825.l 3025. 25333. 220825. 1.97841 106lb=4.81552l 14.0236l 104.284l 820.711l10.4765 - 7.37686 x + 1.41989 x2 - 0.0796299 x3l显然这个拟合很不满意。6.1.2 函数插值函数插值l多项式插值多项式插值 多项式插值是在给定多项式插值是在给定n个数据点个数据点:xi, yi, i=1,2,., n 后后, 求出一个次求出一个次数不超过数不超过n-1的多项式的多项式 (x)作为函数作为函数 y=f(x) 的近似表

20、达式，它就是的近似表达式，它就是计算方法中常说的计算方法中常说的拉格朗日拉格朗日(Lagrange)插值或插值或Newton插值多项式。插值多项式。lMathematica多项式插值的一般形式为多项式插值的一般形式为: InterpolatingPolynomial 数据点集合数据点集合, 自变量名自变量名l具体的多项式插值命令有具体的多项式插值命令有:命令形式命令形式1：InterpolatingPolynomial x1,y1,x2,y2,.,xn,yn,x功能：功能：根据数据点集x1,y1,x2,y2,.,xn,yn求出n-1次插值多项式 (x)命令形式命令形式2：Interpolati

21、ngPolynomial y1,y2,.,yn,x功能：功能：根据数据点集1,y1,2,y2,.,n,yn求出n-1次插值多项式 (x)Lagrange插值插值l求求n次次Lagrange插值多项式算法插值多项式算法l l输入n+1个插值点： (xi, yi),i=0,1,nl计算插值基函数l0n(x), l1n(x) , lnn(x)l给出n次Lagrange插值多项式：Ln(x)= y0 l0n(x)+ y1 l1n(x) +yn lnn(x)l求求Lagrange插值多项式程序插值多项式程序lClearlag,xi,x,yi;lxi=Inputxi=lyi=Inputyi=ln=Leng

22、thxi-1;lp=Sumyii*(Product(x-xij)/(xii-xij),j,1,i-1l*Product(x-xij)/(xii-xij),j,i+1,n+1),i,1,n+1;llagx_=Simplifypl说明：说明：本程序用于求n次Lagrange插值多项式。程序执行后，按要求通过键盘输入插值基点xi:x0 , x1, . , xn 和对应函数值yi: y0 , y1 , , yn 后，程序即可给出对应的n次Lagrange插值多项式lagx。l程序中变量说明程序中变量说明lxi：存放插值基点x0 , x1, . , xn lyi: 存放对应函数值y0 , y1 , ,

23、ynllagx: 存放求出的n次Lagrange插值多项式Ln(x)l注：注：语句lagx_=Simplifyp用简化形式给出对应的n次Lagrange插值多项式。l例例给定数据表 x 0 1 2 3 y=f(x) 1 3 5 12 用Lagrange插值法求三次插值多项式,并给出函数f(x)在x =1.4的近似值。l解：解： 执行Lagrange插值程序后，在输入的两个窗口中按提示分别输入0, 1, 2, 3、1, 3, 5, 12，每次输入后用鼠标点击窗口的“OK”按扭，得如下插值函数。 6 + 22 x - 15 x2 + 5 x3 - 6l所以得到三次插值多项式L3(x)=1+11 x

24、/3-5 x2/2+5 x3/6 接着键入“lag1.4”，则输出3.52，因此f(x)在x =1.4的近似值为3.52，即f(1.4)3.52.Newton插值插值l 求求Newton插值多项式算法插值多项式算法l l1. 输入n+1个插值点： (xi, yi),i=0,1,nl2. 计算差商表l3.给出n次Newton插值多项式。求求Newton插值多项式程序插值多项式程序lClearnewt,s,x;lxi=Inputxi=lyi=Inputyi=ln=Lengthxi;l(*计算差商表*)lf=Table0,n,n;lDofi,1=yii,i,1,nlDofi,j+1=(fi,j-fi

25、+1,j)/(xii-xii+j), j,1,n-1,i,1,n-jlPrint差商表lDoPrintxii, ,fi,i,1,nl(*求Newton插值多项式*)lfa=1;ls=f1,1;lDofa=(x-xik)*fa;s=s+fa*f1,k+1,k,1,n-1lnewtx_=slSimplify%l说明：说明：本程序用于求n次Newton插值多项式。程序执行后，按要求通过键盘输入插值基点xi:x0 , x1, . , xn 和对应函数值yi: y0 , y1 , , yn 后，程序依次给出输入的数据表、计算出的差商表、Newton插值多项式、Newton插值多项式的简化形式。l程序中变

26、量说明：程序中变量说明：lxi：存放插值基点x0 , x1, . , xn lyi: 存放对应函数值y0 , y1 , , ynlf:存放函数值y0 , y1 , , yn及所有差商lnewtx: 存放求出的n次newton插值多项式Nn(x)l注：注：l（1)语句f=Table0,n,n用于产生一个nn的矩阵变元用于存放函数值y0 , y1 , , yn及所有差商。l（2)在Mathematica中有一个求n次插值多项式的命令，命令形式 InterpolatingPolynomialx01,y0，x1,y1,x2,y2,xn,yn, xl它可以求过n+1个插值点x01,y0，x1,y1,x2

27、,y2,xn,yn的n次插值多项式Pn(x)。l例例2给定数据表l x 4.0002 4.0104 4.0233 4.0294l y 0.6020817 0.6031877 0.6045824 0.6052404l(1)计算差商表l(2)用Newton插值法求三次插值多项式Nn(x)l(3)求f(4.011)的近似值l解：解： 执行Newton插值程序后，在输入的两个窗口中按提示分别输入4.0002, 4.0104, 4.0233, 4.0294，0.6020817, 0.6031877, 0.6045824, 0.6052404,每次输入后用鼠标点击窗口的“OK”按扭，得如下插值函数。4.0

28、002, 4.0104, 4.0233, 4.02940.6020817, 0.6031877, 0.6045824, 0.6052404 l差商表4.0002 0.6020817, 0.108431, -0.0136404, 0.02116294.0104 0.6031877, 0.108116, -0.0130225, 04.0233 0.6045824, 0.107869, 0, 04.0294 0.6052404, 0, 0, 00.6020817 + 0.108431 (-4.0002 + x) - 0.0136404 (-4.0104 + x) (-4.0002 + x) + 0.

29、0211629 (-4.0233 + x) (-4.0104 + x) (-4.0002 + x)-1.41642 + 1.23926 x - 0.268313 x2 + 0.0211629 x3l于是我们得到本题的差商表为xi yi f, f , , f , , ,4.0002 0.6020817, 0.108431, -0.0136404, 0.02116294.0104 0.6031877, 0.108116, -0.01302254.0233 0.6045824, 0.107869 4.0294 0.6052404 和三次插值多项式N3(x)= -1.41642 + 1.23926 x

30、 - 0.268313 x2 + 0.0211629 x3接着键入“newt4.011”，则输出0.603253，因此f(4.011) 0.603253。l例例3多项式插值的误差估计式中可以看到，当插值节点越多时误差会越小，这个结论正确吗？通过实验说明该结论的正确性。l解解: 考虑函数f(x) = (1+x2)-1 在区间-4,4内选取不同个数的等距插值节点做观察，这里分别选-4,4 内的9个和11个的等距节点来做实验，将对应的插值函数图与被插函数f(x) = (1+x2)-1画在一起做观察，为简单起见，这里用Mathematica 命令做实验，对应命令为lIn1:= u=Tablex,(1+

31、x2)-1,x,-4,4 ; （*采取f(x) 在-4,4 内的9个插值点 lIn2:= g=ListPlotu, PlotStyle-PointSize0.04 （*将散点图图形文件存放在变量g中lIn3:= s=InterpolatingPolynomialu , x ; （*将插值函数存放在变量s中lIn4:= t= Plots, (1+x2)-1, x,-4,4, PlotStyle-Thickness0.005,Thickness0.006 （*将插值函数s与f(x)画在一起的图形文件存放在变量t中lIn5:= Showt, g（*将散点图, 插值函数s, f(x)画在一个坐标系中l

32、在-4,4中选9个等距节点的插值函数与被插函数图，粗线为被插函数图 lIn6:= u=Tablex,(1+x2)-1,x,-4,4，0.8 ; （*采取f(x) 在-4,4 内的11插值点 lIn7:= g=ListPlotu, PlotStyle-PointSize0.04 lIn8:= s=InterpolatingPolynomialu , x ; lIn9:= t= Plots, (1+x2)-1, x,-4,4, PlotStyle-Thickness0.005,Thickness0.006 lIn10:= Showt, g l在-4,4中选11个等距节点的插值函数与被插函数图 分段

33、多项式插值分段多项式插值Mathematica分段多项式插值的一般形式为分段多项式插值的一般形式为: Interpolation 数据点集合数据点集合具体的分段多项式插值命令有具体的分段多项式插值命令有:l命令形式命令形式1：Interpolation x1,y1,x2,y2,.,xn,yn功能：功能：根据数据点集x1,y1,x2,y2,.,xn,yn求出分段插值多项式(x)l命令形式命令形式2：Interpolation y1,y2,.,yn功能：功能：根据数据点集1,y1,2,y2,.,n,yn求出分段插值多项式(x)例5: 用分段多项式插值重做例4的插值问题, 并验证, 随着插值点的增多

34、, 分段插值函数的误差会越来越小。(见图) 解解:In27:= r=InterpolationTablex,(1+x2)-1, x,-4,4 Out27= InterpolatingFunction-4, 4, In28:= (1+x2)-1/.x-3.7Out28= 0.0680735 In29:= r3.7Out29=0.0734 In30:= d= Tablex,(1+x2)-1, x,-4,4,0.5; In31 := r1=InterpolationdOut31= InterpolatingFunction-4, 4., In32:= r13.7Out32= 0.0681761In3

35、3:= Plotr1x, (1+x2)-1, x,-4,4 下一页下一页返回返回分段线性插值分段线性插值l用分段线性插值函数做近似计算的算法用分段线性插值函数做近似计算的算法l1. 输入n个插值点： (xi, yi),i=1,nl2. 输入要近似计算的自变量点xal3.寻找包含xa的小区间xi , xi+1l4.用 xi , xi+1 上的线性插值函数在xa处的值作为f(xa)的函数值l用分段线性插值函数做近似计算的程序用分段线性插值函数做近似计算的程序lClearx,a,b;llia_,b_,x_:=(b2-a2)/(b1-a1)*(x-a1)+a2lxi=Inputxi=lyi=Input

36、yi=lxa=Inputxa=ln=Lengthxi;lIfxaxin,Print超限;Break,l DoIfxa=xik,m=k,k,1,n-1;l Printm=,m,x=,xim;l lixim,yim,xim+1,yim+1,xal说明：说明：本程序用分段线性插值做近似计算。程序执行后，按要求通过键盘输入插值基点xi: x1, . , xn 和对应函数值yi: y1 , , yn 和要做近似计算的点xa后，程序将给出包含xa小区间xi , xi+1的下标i和xi ,和函数f(xa)的近似值。如果xa不在x1, xn 内程序给出超限提示。l程序中变量说明程序中变量说明lxi：存放插值基

37、点 x1, . , xn lyi: 存放对应函数值 y1 , , ynlxa: 存放要做近似计算的自变量值lm:存放包含xa小区间xi , xi+1的下标il注：在Mathematica中有一个求分段线性插值函数的命令，命令形式为l Interpolation x1,y1,x2,y2,.,xn,yn，InterpolationOrder - 1l用如上Mathematica命令求出的分段线性插值函数没有给出具体的分段函数表达式, 而是用l“InterpolatingFunctionx1, xn, ”作为所求的分段插值函数, 通常可以用l 变量=Interpolation 数据点集合l把所得的分

38、段插值函数存放在变量中, 如果要计算插值函数在某一点的如p点的值,只要输入 “变量p” 即可，如果要表示这个插值函数应该用 “变量x” 。l此外，命令l Interpolation x1,y1,x2,y2,.,xn,yn，InterpolationOrder - 3l可以给出更好的分段插值函数。l例例4给定数据表 xi 1 2 3 4 yi 0 -5 -6 3l（1）用分段线性插值函数求函数f(x)在x =1.4的近似值l（2）用Mathematica命令画出分段线性插值图形l解：解：执行求分段线性插值函数程序后，在输入的两个窗口中按提示分别输入l1，2，3，4，0，-5，-6，3，1.4l每

39、次输入后用鼠标点击窗口的“OK”按扭，得如下插值函数。lm=1 x=1l-2l此结果说明f(x)在x =1.4的近似值为2，它对应的小区间为x1 , x2。为完成第二个问题，键入命令ld=Interpolation 1,0,2,-5,3,-6,4,3,InterpolationOrder - 1lPlotdx,x,1,4l得输出如下。lInterpolatingFunction1, 4, 6.2 解常微分方程l自变量、未知函数以及未知函数的导数（或微分）之间的一个（或一组）关系式称为微分方程（或微分方程组）。求出常微分方程（或微分方程组）的未知函数（或未知函数组），称为解常微分方程。含有任意常

40、数的解称为通解，否则称为特解。n阶常微分方程的一般形式为：l利用Mathematica可以象高等数学一样求出常微分方程的解析解（准确解），如果求不出解析解Mathematica可以求出微分方程的数值解（即近似解。Mathematica解常微分方程的命令有求常微分方程（组）的解析解命令和求常微分方程（组）的数值解命令。 0),()(nyyyxF6.2.1求常微分方程（组）的解析解求常微分方程（组）的解析解命令形式命令形式1：DSolveeqn, yx, x 功能：功能：求出常微分方程eqn的未知函数yx的解析通解。命令形式命令形式2：DSolveeqn1,eqn2, . , y1x,y2x, .

41、 , x 功能：功能：求出常微分方程组eqn1,eqn2, . 的所有未知函数y1x,y2x, . 的解析通解。注意注意：l每个常微分方程的中的等号输入时应该用两个等号代替l未知函数及未知函数的导数要用“自变量名”指示未知函数的自变量l常微分方程组和未知函数组应该用大括号括起来l命令形式2可以用来求常微分方程的特解leqn表示常微分方程, x是自变量名, yx 表示未知函数, 自变量名和未知函数可以是其他的变量名例例6: 求常微分方程y =2ax的通解，a为常数解:Mathematica 命令: In34:= DSolveyx = 2 a x, yx, xOut34= yx - a x 2 +

42、 C1即得本问题的通解 例7: 求常微分方程y +y =0的通解 解:Mathematica 命令: In35:= DSolveyx+yx=0, yx, xOut35= yx - C2 Cosx - C1 Sinx即得本问题的通解:C1，C2是任意常数。 xcxcysincos1221xacy例8: 求常微分方程的特解。 解:Mathematica 命令: In36:= DSolveyx=x/yx+yx/x, y1=2, yx, x Out36= yx - Sqrtx2 (4 + 2 Logx) 本问题的特解为:l下面的操作可以检验所求特解的正确性:In37:= yx_:= Sqrtx2* (

43、4 + 2 Logx) In38:= y1 Out38=2In39:=Dyx,x-x/yx-yx/x; In40:= Simplify% Out40= 0 2|,1xyxyyxy)ln24(2xxy例例9: 求常微分方程组求常微分方程组y (t) =x(t), x (t) =y(t) 的通解。的通解。 解解:Mathematica 命令命令: In41:= DSolveyt=xt,xt=yt, xt,yt,t C1+ E 2 t C1-C2+ E 2 t C2l Out41= xt - -, 2 E t -C1+ E 2 t C1+C2+ E 2 t C2 yt - - 2 E t6.2.2

44、求常微分方程（组）的数值解命令求常微分方程（组）的数值解命令l命令形式命令形式1：NDSolve eqns, y, x, xmin, xmax 功能：功能：求出自变量范围为xmin, xmax且满足给定常微分方程及初值条件eqns的未知函数y的数值解 l命令形式命令形式2：NDSolve eqns, y1,y2, . , x, xmin, xmax 功能：功能：求出自变量范围为xmin, xmax且满足给定常微分方程及初值条件eqns的未知函数y1,y2, . 的数值解。例例10:求微分方程初值问题求微分方程初值问题y =x+y,y(0)=1在区间在区间0,0.5内的数值解内的数值解解解:Ma

45、thematica 命令命令: In42:= NDSolve yx=x+yx, y0=1, y, x,0,0.5 Out42= y - InterpolatingFunction0., 0.5, l上面显示的是所求数值解的替换形式，为得到本问题数值解，再键入上面显示的是所求数值解的替换形式，为得到本问题数值解，再键入:In43:= f= y/.% 1 Out43= InterpolatingFunction0., 0.5, In44:= Plotfx,x,0,0.5 In45:= Table fx, x, 0, 0.5, 0.1 Out45= 1., 1.11034, 1.24281, 1.3

46、9972, 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,10 Out46= x - InterpolatingFunction0., 10., , y - InterpolatingFunction0., 10., In47:= f=x/. q1; g=y/. q 1; In48= f0.4 Out48= 0.375078 In49= g4Out49

47、= -0.0912007In50= ParametricPlotft,gt,t,0,101)0()0(22yxyxyxyxEuler方法方法lEuler方法算法方法算法l l输入函数f(x,y)、初值y0、自变量区间端点a,b步长hl计算节点数n和节点x k l用Euler公式y n+1= yn+h f(xn,yn) 求数值解lEuler方法程序方法程序lClearx,y,flfx_,y_= Input函数f(x,y)=ly0=Input初值y0 =la=Input左端点a=lb=Input右端点b=lh=Input步长h=ln=(b-a)/h;lFori=0,iInterpolatingFunctionrange, 的形式给出的,其中的InterpolatingFunctionrange, 是所求的插值函数表示的数值解, range就是所求数值解的自变量范围。l用Euler方法求初值问题 的数值解。取步长h=0.1,并在一个坐标系中画出数值解与准确解的图形。l解：解：执行Euler方法程序后，在输入的窗口中按提示分别输入y-2x/y、1、0、1、0.1，每次输入后用鼠标点