实验名称微分方程数值解

1、探索实验8 常微分方程初值问题数值解一、 实验目的 了解常微分方程初值问题的数值解概念，掌握解常微分方程初值问题的Euler方法及改进的Euler方法和Runge-Kutta方法解常微分方程初值问题的算法构造和计算。能用程序实现解常微分方程初值问题的Euler方法、改进的Euler方法和经典的Runge-Kutta方法，学习用计算机求常微分方程初值问题数值解的一些科学计算方法和简单的编程技术。二、概念与结论1. 常微分方程初值问题： 常微分方程特解问题称为初值问题，通常其形式为 2常微分方程初值问题数值解： 常微分方程初值问题的解y(x)在a,b上的有限个值y(xk)的近似值yk称为常微分方程

2、初值问题数值解，其中 xk= xk-1 + hk-1 ，k=1,2,3,N。xk称为节点， hk 称为步长。 通常，步长h不变，取为等距步长 hk=h=（b-a）/N，N为等分区间a,b分割数。3.常微分方程初值问题数值解法： 求常微分方程初值问题数值解yk的方法称为微分方程初值问题数值解法。4.单步法： 在计算yk+1之时只用到yk 的方法，其计算公式有： 显式单步法计算公式: yk+1=yk+h(xk,yk;h） 隐式单步法计算公式: yk+1 = yk +h(xk ,yk,yk+1,h） 式中函数是连续函数，称为增量函数。5.多步法： 在计算yk+1之时不仅用到yk ，还要用yk-1，y

3、k-2，,一般m 步法要用到yk,yk-1,yk-2， yk-m+1， 多步法也有显式方法和隐式方法之分。6.数值解法的局部截断误差： 假设某常微分方程初值问题数值解法在x= xk没有误差，即y(xk)= yk，称Tk+1 =y(xk+1)- yk+1为该数值方法的局部截断误差。 显式单步法有其局部截断误差为：Tk+1 =y(xk+1)- y(xk)- h(xk,y(xk),h）7.数值解法的阶： 如果某常微分方程初值问题数值解法的局部截断误差为：Tk+1 =O(h p+1)则称该数值方法的阶为P,或称该数值方法为P阶方法。 数值方法的阶越高，方法越好。三、程序中Mathematica语句解释

4、： 1.k+,+k,k-,-kk+ 表示赋值关系 k = k+1 , 如: k=1;Table+k,5获得表2,3,4,5,6+k 表示先处理k的值,再做赋值 k=k+1, 如: k=1;Tablek+,5获得表1,2,3,4,5k- 表示赋值关系 k = k-1, 如: k=1;Tablek-,5获得表1, 0, -1, -2, -3-k 表示先处理k的值,再做赋值 k=k-1,如:k=1;Table-k,4获得表0,-1,-2,-32. x+=k, x*=kx+=k 表示 x = x + kx*=k 表示 x = x * k3.Forstat,test,incr,body 以stat为初值

5、,重复计算incr和body直到test为False终止 。这里start为初始值,test为条件,incr为循环变量修正式,body为循环体,通常由incr项控制test的变化。 注意: 上述命令形式中的start可以是由复合表达式提供的多个初值,如果循环体生成 Break 语句,则退出For循环; 如果循环体生成Continue 语句,则由incr的增量进入For语句的下一次循环。 四、方法与程序在实际问题中，常常会遇到一些常微分方程初值问题。对这些问题如果只求助于高等数学中介绍的用求通解加确定边界条件的方法去求解往往是无能为力的。因为一方面通解可能求不出来，另一方面所求的解可能是不能用初

6、等函数表达的形式。人们处理这类问题主要采用常微分方程初值问题数值解的方法来求解。这类方法有很多，这里主要介绍解常微分方程初值问题的Euler方法、改进的Euler方法和Runge-Kutta方法的构造过程和算法程序。1. Euler方法Euler方法是最简单的求微分方程数值解的方法。这个方法由于精度不高，实用中较少使用。人们常用Euler方法来说明求微分方程数值解涉及到的一些问题。1.1 Euler方法的构造过程： Euler方法涉及的计算公式有很多方法，这里介绍在求微分方程数值解用的较多的Taylor展开法。 设 f(x,y) 充分光滑, xn+1=xn+h,将y(xn+1)在x n点作Ta

7、ylor展开，得:y (x n+1)=y(xn)+hy¢(xn)+(h2/2!)y²(xn)+O(h3)取其关于h 的线性部分，有：y (x n+1)»y(xn)+hy¢(xn)注意到y¢(xn)= f(xn,y(xn),用yk代替y(xk)并将“»”换为等号“=”,得到Euler公式：y n+1= yn+h f(xn,yn) 于是，由初始条件y0=y(a),借助Euler公式就可以依次计算出微分方程初值问题的数值解：y1 ,y2 ,yk称由Euler公式求数值解的方法为Euler方法。显然Euler方法是单步显式方法。 因为对Eul

8、er公式有局部截断误差为Tn+1= O(h2)因此Euler方法是一阶方法。1.2 Euler方法算法： 1 输入函数f(x,y)、初值y0、自变量区间端点a,b步长h2 计算节点数n和节点x k 3 用Euler公式y n+1= yn+h f(xn,yn) 求数值解1.3 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;F

9、ori=0,i<n,i+,xk=a+i*h;y1=y0+h*fxk,y0;Print"y(",xk+h/N,")=",y1/N;y0=y1说明：本程序用Euler公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微方程初值问题数值解。程序中变量说明：fx,y: 存放函数f(x,y)y0: 存放初值y0及数值解a：存放自变量区间左端点b: 存放自变量区间右端点n: 存放节点个数h: 存放节点步长xk：存放节点xiy1: 存放数值解注：1）语句Print"

10、;y(",xk+h/N,")=",y1/N是将数值解用6位有效数字显示出来，如果要显示n位数有效数字，可以将语句改为Print"y(",xk+h/N,")=",Ny1,n2）Mathematica中有求微分方程初值问题数值解的命令，形式为:NDSolve yx=fx,yx,ya=y0, y, x, a, b 由NDSolve命令得到的解是以未知函数名->InterpolatingFunctionrange, <>的形式给出的,其中的InterpolatingFunctionrange, <>是所

11、求的插值函数表示的数值解, range就是所求数值解的自变量范围。1.4 例题与实验例1. 用Euler方法求初值问题的数值解。取步长h=0.1,并在一个坐标系中画出数值解与准确解的图形。解：执行Euler方法程序后，在输入的窗口中按提示分别输入x-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

12、(1.)=1.78477本题的准确解为y(x)= (1 + 2 x)1/2，在一个坐标系中画出数值解与准确解的图形如下：图中点图是数值解，曲线为准确解。从图中可以知道数值解与准确解比较接近，但还是有误差的。2. 改进的Euler方法 改进的Euler方法比Euler方法精度高，其中把微分方程初值问题数值隐式解法计算公式变为预测、校正公式的做法很有代表性。2.1改进的Euler方法构造过程： 将微分方程y¢=f(x,y) 在xk,xk+1积分，得对右端的定积分用数值积分的梯形公式，有用yk代替y(xk)并将“»”换为等号“=”可得求微分方程初值问题的数值解的梯形公式：这个公式

13、是单步隐式公式。可以证明它是二阶方法，比Euler方法阶高。不过，在给定初始条件y0=y(a)后,要求出数值解，每一次计算yk的值都要进行非线性方程求根的迭代解法来完成，因此计算量很大。为减少计算量，通常采用迭代一两次就转入下一步计算的方法，特别，如果先用Euler公式计算一次以对要计算的下一步值解进行预测，然后再用梯形公式对其进行校正的方法得到下一步的值，它的计算格式为：这个计算格式称为改进的Euler公式，而称由如上计算格式求数值解的方法为改进的Euler方法。2.2 改进的Euler方法算法：1.输入函数f(x,y)、初值y0、自变量区间端点a,b步长h2.计算节点数n和节点x k 3.

14、用改进的Euler公式求数值解2.3 改进的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,i<n,i+,xk=a+i*h;y1=y0+h*fxk,y0;y1=y0+h/2*(fxk,y0+fxk+h,y1);Print"y(",xk+h/N,")=",y1/N;y

15、0=y1注：语句h=Input"步长h="的步长输入应该用带小数点的数输入，这样可以加快计算速度，否则，由于Mathematica软件本身的原因，可能会出现计算时间过长等计算不出结果的情况。说明：本程序用改进的Euler公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微方程初值问题数值解。程序中变量说明：fx,y: 存放函数f(x,y)y0: 存放初值y0及数值解a：存放自变量区间左端点b: 存放自变量区间右端点n: 存放节点个数h: 存放节点步长xk：存放节点xiy1: 存放数值解

16、2.5例题与实验例2. 用改进的Euler方法求初值问题的数值解。取步长h=0.1和0.01计算，并与准确解比较。解：执行改进的Euler方法程序后，在输入的窗口中按提示分别输入-2x*y2、1、0、1、0.1，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果：y(0.1)=0.99 y(0.2)=0.961366 y(0.3)=0.917246 y(0.4)=0.861954y(0.5)=0.800034 y(0.6)=0.735527 y(0.7)=0.671587 y(0.8)=0.610399y(0.9)=0.553289 y(1.)=0.500919再执行一次

17、改进的Euler方法程序后，在输入的窗口中按提示分别输入-2x*y2、1、0、1、0.05，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果：y(0.05)=0.9975 y(0.1)=0.990087 y(0.15)=0.977978 y(0.2)=0.961519y(0.25)=0.941158 y(0.3)=0.917417 y(0.35)=0.890863 y(0.4)=0.862076y(0.45)=0.831624 y(0.5)=0.800043 y(0.55)=0.767819 y(0.6)=0.735383y(0.65)=0.7031 y(0.7)=0.

18、671277 y(0.75)=0.640158 y(0.8)=0.609935y(0.85)=0.580748 y(0.9)=0.552699 y(0.95)=0.52585 y(1.)=0.500236本题的准确解为y(x)=(1+x2)-1它在0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0的值为y(0.1)=0.990099, y(0.2)=0.961538, y(0.3)=0.917431, y(0.4)=0.862069,y(0.5)=0.8, y(0.6)=0.735294, y(0.7)=0.671141, y(0.8)=0.609756, y(0.

19、9)= 0.552486, y(1.0)= 0.5,比较上面的计算结果，显然，步长为0.05的结果比步长为0.1的结果要好。3.Runge-Kutta方法Runge-Kutta方法是显式单步高阶的求微分方程数值解的方法，其中四阶Runge-Kutta方法最为常用。四阶Runge-Kutta方法又称为经典的Runge-Kutta方法，它具有精度高、程序简单，容易编程、计算过程稳定饿易于调节步长的特点。3.1 Runge-Kutta方法的构造过程： 将将微分方程y¢=f(x,y) 在xk,xk+1积分并利用积分中值定理，有：为提供增加求解的精度，把f(x,y(x)写成一个线性组合的形式并

20、用yk代替y(xk)，就得到Runge-Kutta方法的一般形式：上式中选择不同的m,ci,x i值就得到不同的Runge-Kutta计算公式。通常，为方便获得Runge-Kutta计算公式，常把Runge-Kutta方法的一般形式写为：于是，利用二元Taylor展开将公式中的f在(xk,yk)展开并适当选择参数就可以得到具体的Runge-Kutta计算公式了。称由Runge-Kutta计算公式式求数值解的方法为Runge-Kutta方法。下面以m=2的二阶Runge-Kutta计算公式为例来说明获得Runge-Kutta计算公式的方法。 m=2的二阶Runge-Kutta计算公式为它的增量函

21、数为j (x,y(x),h)=c1f(x,y(x)+ c2f(x+a2h,y(x)+hb21f(x,y(x)，考虑它的局部截断误差Tk+1 =y(xk+1)- y(xk)- hj (xk,y(xk),h）。利用f(x+a2h,y(x)+hb21f(x,y(x) 在(x,y)做二元Taylor展开使其具阶为2，则有：这是有四个参数三个方程的方程组，它有无穷多组解。例如，取c1=0,可以得到c2=1，a2=0.5, b21=0.5,于是得到一个m=2的二阶Runge-Kutta计算公式：它称为中点公式。类似地，可以得到很多其他Runge-Kutta计算公式。经典的Runge-Kutta计算公式是四

22、阶的，形式为：于是，由初始条件y0=y(a),借助经典的Runge-Kutta计算公式就可以依次计算出微分方程初值问题的数值解：y1 ,y2 ,yn。3.2 Runge-Kutta方法算法： 4 输入函数f(x,y)、初值y0、自变量区间端点a,b步长h5 计算节点数n和节点x k 6 用Runge-Kutta r公式求数值解3.3 Runge-Kutta方法程序：Clearx,y,ffx_,y_= Input"函数f(x,y)="y0=Input"初值y0 ="a=Input"左端点a="b=Input"右端点b=&quo

23、t;h=Input"步长h="n=(b-a)/h;Fori=0,i<10,i+,xk=a+i*h;k1=fxk,y0;k2=fxk+h/2,y0+k1*h/2;k3=fxk+h/2,y0+k2*h/2;k4=fxk+h,y0+k3*h;y1=y0+h/6*(k1+2k2+2k3+k4);Print"y(",xk+h/N,")=",y1/N;y0=y1说明：本程序用经典Runge-Kutta公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微

24、方程初值问题数值解。程序中变量说明：fx,y: 存放函数f(x,y)y0: 存放初值y0及数值解a：存放自变量区间左端点b: 存放自变量区间右端点n: 存放节点个数h: 存放节点步长xk：存放节点xiy1: 存放数值解注：语句h=Input"步长h="的步长输入应该用带小数点的数输入，这样可以加快计算速度，否则，由于Mathematica软件本身的原因，可能会出现计算时间过长等计算不出结果的情况。1.5 例题与实验例3. 用经典的Runge-Kutta方法求初值问题的数值解。取步长h=0.1，并与准确解比较。解：执行经典的Runge-Kutta方法程序后，在输入的窗口中按提

25、示分别输入-2x*y2、1、0、1、0.1，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果：y(0.1)=0.990099y(0.2)=0.961538y(0.3)=0.917431y(0.4)=0.862068y(0.5)=0.799999y(0.6)=0.735294y(0.7)=0.671141y(0.8)=0.609756y(0.9)=0.552487y(1.)=0.500001本题的准确解为y(x)=(1+x2)-1，它在0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0的值为y(0.1)=0.990099, y(0.2)=0.9

26、61538, y(0.3)=0.917431, y(0.4)=0.862069,y(0.5)=0.8, y(0.6)=0.735294, y(0.7)=0.671141, y(0.8)=0.609756, y(0.9)= 0.552486, y(1.0)= 0.5,比较上面的计算结果，可见，计算结果比改进的Euler法要好得多。例4. 用经典的Runge-Kutta方法求初值问题的数值解。分别取步长h=0.1，0.025,0.01计算并与准确解在xi=0.1i,i=0,1,2,10各点比较。解：执行经典的Runge-Kutta方法程序后，在输入的窗口中按提示分别输入-50y+50x2+2x、1

27、/3、0、1、0.1，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出步长h=0.1的如下数值解结果：y(0.1)=4.60549 y(0.2)=63.0625 y(0.3)=864.049 y(0.4)=11843.6 y(0.5)=162355.y(0.6)=2.22561´106 y(0.7)=3.05094´107 y(0.8)=4.18232´108y(0.9)=5.73327´109 y(1.)=7.85936´1010 再执行一次经典的Runge-Kutta方法程序，在输入的窗口中按提示分别输入-50y+50x2+2x、1/3、0、1、0.025，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出步长h=0.025的部分数值解结果：y(0.1)=0.0130149y(0.2)=0.0400633y(0.3)=0.090037y(0.4)=0.160037y(0.5)=0.250037y(0.6)=0.360037y(0.7)=0.490037y(0.8)=0.640037y(0.9)=0.810037y(1.)=1.00004又执行一次经