实验名称微分方程数值解剖析

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

2、解y(x)在a,b上的有限个值y(Xk)的近似值yk称为常微分方程初值问题数值解，其中x k= xk-1 + h k-1 , k=1,2,3,N。Xk称为节点，hk称为步长。通常，步长h不变，取为等距步长h k=h= ( b-a) /N , N为等分区间a,b分割数。3. 常微分方程初值问题数值解法：求常微分方程初值问题数值解yk的方法称为微分方程初值问题数值解法。4单步法：在计算yk+1之时只用到yk的方法，其计算公式有：显式单步法计算公式：y k+1=yk+h(x k,y k；h )隐式单步法计算公式：y k+1 = y k +h(x k ,y k,y k+1,h )式中函数是连续函数，称

3、为增量函数。5. 多步法：在计算yk+1之时不仅用到yk，还要用yk-1 , yk-2,一般m步法要用到yk,y k-1 ,y k-2 ,yk-m+1, 多步法也有显式方法和隐式方法之分。6. 数值解法的局部截断误差：假设某常微分方程初值问题数值解法在x= x k没有误差，即y(x k)= y k,称Tk+1 =y(x k+1)- y k+1为该数值方法的局部截断误差。显式单步法有其局部截断误差为：Tk+1 =y(x k+1)- y(x k)- h (x k,y(x k),h )7. 数值解法的阶：如果某常微分方程初值问题数值解法的局部截断误差为：Tk+1 =0(h p+1)则称该数值方法的阶

4、为P,或称该数值方法为P阶方法。数值方法的阶越高，方法越好。、程序中 Mathematica 语句解释:1. k+,+k,k,k k+ 表示赋值关系 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

5、 表示 x = x * k3. Forstat,test, in cr,body以stat为初值，重复计算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展开法。

7、设f(x,y)充分光滑，Xn+i=Xn+h,将y(Xn+i)在x n点作Taylor展开，得:23y (x n+i)=y(x n)+hy (xn)+(h /2!)y (xn)+)(h )取其关于h的线性部分，有：y (x n+1 )y(xn)+hy (xn)注意到y(xn)= f(x n,y(xn),用yk代替y(xk)并将“ 换为等号=”,得到Euler公式：y n+1 = yn+h f(x n,yn)于是，由初始条件yo=y(a)，借助Euler公式就可以依次计算出微分方程初值问题的数值解:y1 ,y2,y称由Euler公式求数值解的方法为Euler方法。显然Euler方法是单步显式方法。

8、因为对Euler公式有局部截断误差为2Tn+1= O(h )因此Euler方法是一阶方法。1.2 Euler方法算法:1. 输入函数f(x,y)、初值yo、自变量区间端点a,b步长h2. 计算节点数n和节点x k3. 用 Euler 公式 y n+i= 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;Fori=0,i lnterpolatingFunctionrange, 的形 式给出

9、的，其中的InterpolatingFunctionrange, 是所求的插值函数表示的数值解 ，range就 是所求数值解的自变量范围。1.4例题与实验例1. 用Euler方法求初值问题0 x ： 1, 2x y =y_ y.y(0) =1的数值解。取步长 h=0.1,并在一个坐标系中画出数值解与准确解的图形。解：执行Euler方法程序后，在输入的窗口中按提示分别输入x-2x/y、1、0、1、0.1，每次输入后用鼠标点击窗口的“0K”按扭，计算机在屏幕上给出如下数值解结果：y(0.1)=1.1y(0.2)=1.19182y(0.3)=1.27744y(0.4)=1.35821y(0.5)=1

10、.43513y(0.6)=1.50897y(0.7)=1.58034y(0.8)=1.64978y(0.9)=1.71778y(1.)=1.78477本题的准确解为y(x)= (1 + 2 x) 1/2，在一个坐标系中画出数值解与准确解的图形如下：图中点图是数值解，曲线为准确解。从图中可以知道数值解与准确解比较接近, 差的。但还是有误2.改进的Euler方法改进的Euler方法比Euler方法精度高，其中把微分方程初值问题数值隐式解法计算公式 变为预测、校正公式的做法很有代表性。2.1改进的Euler方法构造过程：将微分方程y =f(x,y)在xk,xk+i积分，得.xk+y(Xk 1) -

11、y(xj = x f (x, y(x)dxxk对右端的定积分用数值积分的梯形公式，有f (x, y(x)dx ：号(f(Xk, y(Xk) f(Xk 1, y(Xk 1)用yk代替y(Xk)并将“：”换为等号“=”可得求微分方程初值问题的数值解的梯形公式:y(Xk 1) - y(Xk)=Xk 1Xkyk .1 = yk- (f (Xk, yk) f (Xk 1, yk 1)这个公式是单步隐式公式。可以证明它是二阶方法，比Euler方法阶高。不过，在给定初始条件yo=y(a)后,要求出数值解，每一次计算yk的值都要进行非线性方程求根的迭代解法来完成，因此计算量很大。为减少计算量，通常采用迭代一两

12、次就转入下一步计算的方法，特别,如果先用Euler公式计算一次以对要计算的下一步值解进行预测，然后再用梯形公式对其进 行校正的方法得到下一步的值，它的计算格式为：?k 1 二 ykhf(Xk，yk)h yk 1 訥 -(f(xk ,yk) f (xk 1, ?k 1)这个计算格式称为改进的Euler公式，而称由如上计算格式求数值解的方法为改进的Euler方法。2.2改进的Euler方法算法：1输入函数f(x,y)、初值y0、自变量区间端点a,b步长h2计算节点数n和节点x k3.用改进的Euler公式求数值解2.3 改进的Euler方法程序：Clearx,y,ffx_,y_= Input函数

13、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);Prin ty(”,xk+h/N,)=,y1/N;y0=y1注：语句h=Input”步长h=的步长输入应该用带小数点的数输入，这样可以加快计算速度， 否则，由于Mathematica软件本身的原因，可能会出现计算时间过长等计算不出结果的情况。 说明：本程序用改进的Euler公式求常微方程初值问题数值解。程序执行后

14、，按要求通过键盘依次输入输入函数f(x,y)、初值yo、自变量区间端点 a,b步长h后，计算机则给出常微方程初值问题数值解。程序中变量说明：fx,y:存放函数 f(x,y)y0:存放初值yo及数值解 a：存放自变量区间左端点 b:存放自变量区间右端点 n:存放节点个数 h:存放节点步长xk :存放节点 xiy1:存放数值解2.5例题与实验例2. 用改进的Euler方法求初值问题2y=2xyOcxcl、y(O)=1的数值解。取步长 h=O.1和O.O1计算，并与准确解比较。解：执行改进的Euler方法程序后，在输入的窗口中按提示分别输入-2x*yA2、1、O、1、O.1，每次输入后用鼠标点击窗口

15、的“OK ”按扭，计算机在屏幕上给出如下数值解结果：y(O.1)=O.99y(O.5)=O.8OOO34y(O.9)=O.553289y(O.2)=O.961366y(O.6)=O.735527y(1.)=O.5OO919y(O.3)=O.917246y(O.7)=O.671587y(O.4)=O.861954y(O.8)=O.61O399再执行一次改进的Euler方法程序后，在输入的窗口中按提示分别输入-2x*yA2、1、0、1、O.O5,每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果：y(0.05)=0.9975 y(0.1)=0.990087y(0.15)=0.

16、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.7031y(0.7)=0.671277y(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本题的准确解为2 -1y(x)=(1+x )它在

17、 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.9)= 0.552486,y(1.0)= 0.5,比较上面的计算结果，显然，步长为0.05的结果比步长为0.1的结果要好。3.Runge-Kutta 方法Runge-Kutta方法是显式单步高阶的求微分方程数值解的方法，其中四阶Runge-

18、Kutta方法最为常用。四阶 Runge-Kutta方法又称为经典的 Runge-Kutta方法，它具有精度高、程序 简单，容易编程、计算过程稳定饿易于调节步长的特点。3.1 Runge-Kutta方法的构造过程：将将微分方程y=f(x,y)在xk,xk+1积分并利用积分中值定理，有：y(xk 十)一y( Xk) = f + f (x, y(x)dx 知(, y(專)Xkyk代替y(Xk)，就得到为提供增加求解的精度，把f( ,y()写成一个线性组合的形式并用Run ge-Kutta方法的一般形式：myk 1 *k h Ci f ( i, y( i )i =上式中选择不同的m,G, i值就得到

19、不同的Runge-Kutta计算公式。通常，为方便获得Runge-Kutta计算公式，常把 Runge-Kutta方法的一般形式写为:yk41 =myk +h瓦 qKji丄K1 =fgyk)r VKr =f(Xk +ah, yk +h瓦 bjKj) (r=1,2,，m) j丄于是，利用二元 Taylor展开将公式中的f在(Xk,yQ展开并适当选择参数就可以得到具体的Run ge-Kutta计算公式了。称由 Run ge-Kutta计算公式式求数值解的方法为Run ge-Kutta方法。下面以m=2的二阶Runge-Kutta计算公式为例来说明获得Runge-Kutta计算公式的方法。m=2的二

20、阶Runge-Kutta计算公式为yk 仃yk h(ciKi C2K2)2=f(x,yk)K2= f(Xk+a2hyk+hlK1)它的增量函数为 申(x,y(x),h)=c if(x,y(x)+ c 2f(x+a2h,y(x)+hb 2if(x,y(x)，考虑它的局部截断误差 Tk+1 =y(x k+i)- y(x k)- h 9 (x k,y(x k),h)。利用 f(x+a 2h,y(x)+hb 2if(x,y(x) 在(x,y)做二元 Taylor 展开使其具阶为2，则有：Ci + C2 = 11* 2 -c2b21 = 01 c pa 0j 222这是有四个参数三个方程的方程组，它有无

21、穷多组解。例如，取C1=0,可以得到C2=1，a2=0.5,b21=0.5,于是得到一个 m=2的二阶Runge-Kutta计算公式：丄h丄hyk 1 二yk h(f(xk -, yk 2 f(xk,yk)它称为中点公式。类似地，可以得到很多其他Run ge-Kutta计算公式。经典的Runge-Kutta计算公式是四阶的，形式为：y1K1 K2K3K4Yk罷2勺2K3心)6fdk)f(xk 皿 k1)22f(xk 皿 k2)22f(Xk h, yk hK3)于是，由初始条件yo=y(a)，借助经典的Runge-Kutta计算公式就可以依次计算出微分方程初值 冋题的数值解：yi ,y2 ,y。

22、3.2 Runge-Kutta 方法算法：4. 输入函数f(x,y)、初值yo、自变量区间端点a,b步长h5. 计算节点数n和节点x k6. 用 Runge-Kutta r 公式求数值解3.3 Runge-Kutta 方法程序：Clearx,y,ffx_,y_= Input 函数 f(x,y)= yo=Input 初值 yo = a=Input 左端点 a= b=Input 右端点 b= h=Input 步长 h= n=(b-a)/h;Fori=o,i1o,i+, xk=a+i*h;k1=fxk,yo; k2=fxk+h/2,yo+k1*h/2; k3=fxk+h/2,yo+k2*h/2; k

23、4=fxk+h,yo+k3*h; y1=yo+h/6*(k1+2k2+2k3+k4); Printy(,xk+h/N,)=,y1/N;yo=y1说明： 本程序用经典 Runge-Kutta 公式求常微方程初值冋题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值yo、自变量区间端点a,b步长h后，计算机则给出常微方程初值冋题数值解。程序中变量说明：fx,y: 存放函数 f(x,y)y0:存放初值yo及数值解 a：存放自变量区间左端点 b: 存放自变量区间右端点 n: 存放节点个数 h: 存放节点步长xk :存放节点 xiy1: 存放数值解 注： 语句 h=Input 步长 h

24、= 的步长输入应该用带小数点的数输入，这样可以加快计算速度， 否则，由于 Mathematica 软件本身的原因， 可能会出现计算时间过长等计算不出结果的情况。1.5 例题与实验例3. 用经典的Runge-Kutta方法求初值问题r 2y=_2xyOcxcly(O) =1的数值解。取步长 h=0.1，并与准确解比较。解：执行经典的Runge-Kutta方法程序后，在输入的窗口中按提示分别输入-2x*yA2、1、0、1、0.1，每次输入后用鼠标点击窗口的“0K ”按扭，计算机在屏幕上给出如下数值解结果：y(0.1)=0.990099y(0.2)=0.961538y(0.3)=0.917431y(

25、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+x )-，它在 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,

26、 y(0.9)= 0.552486,y(1.0)= 0.5,比较上面的计算结果，可见，计算结果比改进的Euler法要好得多。例4. 用经典的Runge-Kutta方法求初值问题y - -50y 50x2 2x 0 _x_11/3、0、1、0.025，每次输入后用鼠标点击窗口的“0K ”按扭，计算机在屏幕上给出步长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又执行一次经典的 Runge-Kutta方法程序，在输入的窗口中按提示分别输入-50y+50xA2+2x、1/3、0、1、0.01，每次输入后用鼠标点击窗口的 “OK ”按扭，计算机在屏幕上给出步长 h=0.01 的部分数值解结果