常微分方程数值解

1、常微分方程的数值解法常微分方程的数值解法常微分方程的数值解常微分方程的数值解n引言引言n简单的数值方法简单的数值方法nRunge-Kutta方法方法n一阶常微分方程组一阶常微分方程组n高阶方程高阶方程n微分方程边值问题数值解微分方程边值问题数值解在高等数学中我们见过以下常微分方程：在高等数学中我们见过以下常微分方程：1、引言、引言 0)(),()1(yaybxayxfy )(,)(),()2(0ayyaybxayyxfy nybyyaybxayyxfy)(,)(),()3(0(1)，(2)式称为初值问题，式称为初值问题，(3)式称为边值问题。式称为边值问题。在实际应用中还经常需要求解常微分方程

2、组：在实际应用中还经常需要求解常微分方程组： 2002212210012111)(),()(),()4(yxyyyxfyyxyyyxfy本章主要研究本章主要研究问题问题(1)的数值解法，对的数值解法，对(2)(4)只只作简单介绍。作简单介绍。 00)(),()1(yxyyxfy（其中（其中L为为Lipschitz常数）则初值问题（常数）则初值问题（1）存）存在唯一的连续解。在唯一的连续解。 2121),(),(yyLyxfyxf 考虑一阶常微分方程初值问题考虑一阶常微分方程初值问题其中，其中，y = y(x) 是未知函数，是未知函数，y(x0) = y0 是初值是初值条件，而条件，而f (x,

3、 y) 是给定的二元函数是给定的二元函数. 由常微分方程理论知，若由常微分方程理论知，若f(x)在在x a,b连续且连续且 f 满足对满足对 y 的的Lipschitz条件：条件：n常微分方程的数值解法有常微分方程的数值解法有单步法单步法和和多步法多步法之分：之分：q单步法：在计算单步法：在计算yn1 时只用到前一点时只用到前一点yn 的值的值 ；q多步法：计算多步法：计算yn1 时不仅利用时不仅利用yn，还要利用，还要利用yn-1, yn-2,.一般一般k步法要用到步法要用到 yn, yn-1, yn-2,., yn-k+1。求问题（求问题（1）的数值解，就是）的数值解，就是要寻找解函数在一

4、要寻找解函数在一系列离散节点系列离散节点x1 x2 xn xn+1 上的近似值上的近似值y1, y 2,yn 。为了计算方便，可取为了计算方便，可取xn=x0+nh，(n=0,1,2,)， h称为步长。称为步长。 2、 简单的数值方法简单的数值方法一、欧拉（一、欧拉（Euler）方法）方法在在x= x0 处，用差商代替导数：处，用差商代替导数：hxyxyxxxyxyxy)()()()()(0101010 由由00000)(),()(yxyyxfxy 10001),()(yyxhfyxy 得得同理，在同理，在x= xn 处，用处，用差商差商代替代替导数导数：hxyxyxxxyxyxynnnnnn

5、n)()()()()(111 ),()()(1nnnnyxhfxyxy ),()(nnnyxfxy 由由得得若记若记11)(,)( nnnnyxyyxy则上式可记为则上式可记为),(1nnnnyxhfyy 此即为求解初值问题的此即为求解初值问题的Euler方法方法，又称，又称显式显式Euler方法方法。Pn+1yOxx0 x1x2xnP0P1P2Pny=f(x)xn+1Euler方法的几何意义：方法的几何意义：（Euler折线法）折线法）例例: 用用Euler方法求解常微分方程初值问题方法求解常微分方程初值问题 .0)0()30(22yxyxyy并将数值解和该问题的解析解比较。并将数值解和该问

6、题的解析解比较。21)(xxxy 解解析析解解：解：解：Euler方法的具体格式：方法的具体格式：)2(21nnnnnyxyhyy xn y(xn) yn yn-y(xn)0.00000.20.19230.20000.00770.40.34480.38400.03920.60.44120.51700.07580.80.48780.58240.09461.00.50000.59240.09241.20.49180.57050.07871.40.47300.53540.0624取取h=0.2, xn=nh,(n=0,1,2,15), f(x,y)=y/x 2y2 计算中取计算中取f(0,0)=1.

7、 计算结果如下：计算结果如下：xn y(xn) yn yn-y(xn)1.60.44940.49720.04781.80.42450.46050.03592.00.40000.42680.02682.20.37670.39660.01992.40.35500.36980.01472.60.33510.34590.01082.80.31670.32460.00793.00.30000.30570.0057由表中数据可以看到，微分方程初值问题的数值解由表中数据可以看到，微分方程初值问题的数值解和解析解的误差一般在小数点后第二位或第三位小和解析解的误差一般在小数点后第二位或第三位小数上，这说明数上，

8、这说明Euler方法的精度是比较差的。方法的精度是比较差的。 : 数值解数值解 : 准确解准确解 h=0.2;y(1)=0.2;x=h:h:3;for n=1:14 xn=x(n);yn=y(n); y(n+1)=yn+h*(yn/xn-2*yn*yn);endx0=0:h:3;y0=x0./(1+x0.2);plot(x0,y0,x,y,x,y,o)00.511.522.5300.20.40.60.8若直接对若直接对y=f(x,y)在在xn, xn+1积分，积分，dxxyxfxyxynnxxnn 1)(,()()(1利用数值积分中的利用数值积分中的左矩形公式左矩形公式：)(,()(,(1nn

9、xxxyxfhdxxyxfnn ),(1nnnnyxhfyy 此即为此即为Euler公式公式。设设y(xn)= yn，则得则得若用若用右矩形公式右矩形公式：)(,()(,(111 nnxxxyxfhdxxyxfnn),(111 nnnnyxhfyy得得),()0(1nnnnyxhfyy ),()(11)1(1knnnknyxhfyy ),(),(11)(111)1(1 nnknnnknyxfyxfhyy1)(1 nknyyhL二、梯形方法二、梯形方法dxxyxfxyxynnxxnn 1)(,()()(1由由 )(,()(,(2)(,(111 nnnnxxxyxfxyxfhdxxyxfnn ),

10、(),(2111 nnnnnnyxfyxfhyy得得),()0(1nnnnyxhfyy ),(),(211)(111)1(1 nnknnnknyxfyxfhyy1)(12 nknyyhL ),(),(2)(11)1(1knnnnnknyxfyxfhyy 由以上分析可以看出，隐式方法的计算比显由以上分析可以看出，隐式方法的计算比显式方法复杂，需要用迭代法求解非线性方程式方法复杂，需要用迭代法求解非线性方程才能得出计算结果。才能得出计算结果。可采用将可采用将显式显式Euler格式格式与与梯形格式梯形格式结合使用结合使用的方法来避免求解非线性方程。的方法来避免求解非线性方程。记记),(1nnnnyx

11、hfyy 再用梯形格式计算：再用梯形格式计算： ),(),(2111 nnnnnnyxfyxfhyy预测预测校正校正上面两式统称上面两式统称预测校正法预测校正法，又称又称改进的改进的Euler方法方法。三、单步法的局部截断误差和精度三、单步法的局部截断误差和精度单步法的一般形式为单步法的一般形式为 ： （ 与与f 有关）有关）),(11hyyxhyynnnnn 显式单步法形式为：显式单步法形式为：),(1hyxhyynnnn整体截断误差整体截断误差：从：从x0开始，考虑每一步产生的开始，考虑每一步产生的误差，直到误差，直到xn，则有误差则有误差nnnyxye)(称为数值方法在节点称为数值方法在

12、节点xn处的整体截断误差。处的整体截断误差。但但en不易分析和计算，故只考虑从不易分析和计算，故只考虑从xn到到xn+1的局部的局部情况。情况。定义定义：设：设y(x)是初值问题是初值问题(1)的精确解，则称的精确解，则称),(,()()(11hxyxhxyxyTnnnnn 为显式单步法在节点为显式单步法在节点xn+1处的处的局部截断误差局部截断误差 。若存在最大整数若存在最大整数p使局部截断误差满足使局部截断误差满足)(),()()(11 pnhOhyxhxyhxyT 则称显式单步法具有则称显式单步法具有p阶精度阶精度或称或称p阶方法阶方法。注：将注：将Tn+1表达式各项在表达式各项在xn处

13、作处作Taylor展开，可展开，可得具体表达式。得具体表达式。Euler方法的局部截断误差方法的局部截断误差：)(,()()(11nnnnnxyxhfxyxyT )()()(nnnxyhxyhxy )()()(2)()(2nnnnnxyhxyyhxyhxy )()(2)(2322hOxyhyhnn ),(1 nnnxx )()2/(2nxyh Tn+1O(h2)ynxn其中其中（设（设yn=y(xn)）)()(2)()(111 nnnnnxyxyhxyxyT)(! 3)(2)(32nnnxyhxyhxyh )()(2)()()(242hOxyhxyhxyxyhnnnn )()(1243hOxy

14、hn 故故Tn+1O(h3)()12/(3nxyh 3、 RungeKutta方法方法一、一、Runge-Kutta方法的基本思想方法的基本思想)(!)(2)()()()()(211nppnnnnnnxyphxyhxyhxyyhxyxy Tn+1O(hp+1),(yxfy ),(),(),()(yxfyxfyxfxyyx ),(111nnnnyxfKhKyy ),(),()22(121211hKyhxfKyxfKKKhyynnnnnnRunge-Kutta方法的基本思想方法的基本思想： ), 3 , 2(),(),(11111riKhyhxfKyxfKKchyyijjijnininnriiin

15、n 上面第一个式子的右端在上面第一个式子的右端在(xn,yn)作泰勒展开后作泰勒展开后,按按h的幂次作升序排列的幂次作升序排列 ： 332211! 31! 21hhhyynn )()()()(1nnnnxyhxyhxyxy )()(!)(! 21)(2 pnppnhOxyphxyh例如，要求例如，要求)1(321, pnpnnnffff 就得到就得到p个方程个方程,从而定出参数从而定出参数ci ， i， ij ，再，再代入代入K1，K2， Kr的表达式，就可得到计算的表达式，就可得到计算微分方程初值问题的数值计算公式微分方程初值问题的数值计算公式 ： riiinnKchyy11上式称为上式称为

16、r 级级RungeKutta方法的计算公式。方法的计算公式。当当r=1时，就是时，就是Euler方法。方法。要使要使Runge-Kutta公式具有更高的阶公式具有更高的阶p，就要增就要增加加r 的值。下面我们只就的值。下面我们只就r =2推导推导R-K方法。方法。二、二阶二、二阶Runge-Kutta方法方法 ),(),()(12122122111hKyhxfKyxfKKcKchyynnnnnn 上式的局部截断误差为上式的局部截断误差为：),(),()(2122111nnnnnnnnhfyhxfcyxfchyxyT )(! 2)(321hOyhyhyxynnnn nnnynnxnnnnnnnf

17、yxfyxfxyxfdxdyfyxfy),(),()(,(),()(),(),(),(2212212hOhfyxfhyxffhfyhxfnnnynnxnnnn 利用二元函数的利用二元函数的Taylor展开，得展开，得代入代入Tn+1的表达式，得的表达式，得)(),(),()(),(23212212hOhfyxfhyxffcfchfyxfyxfhhfnnnynnxnnnnnynnxn )(),()21(),()21()1(3221222221hOhfyxfchyxfchfccnnnynnxn ),(),()(2122111nnnnnnnnhfyhxfcyxfchyxyT ),(),()(! 22

18、122132nnnnnnnnnhfyhxfcyxfchyhOyhyhy 即即 021021012122221 cccc 212112122221 ccccc1 = 1-a ， 2 = 21 =1/(2a)要使上式要使上式p=2阶，则需阶，则需方程组解不唯一，可令方程组解不唯一，可令c2=a 0 ，则则 满足上述条件的公式都为满足上述条件的公式都为2阶阶R-K公式公式。 )2,2(),(12121KhyhxfKyxfKhKyynnnnnn称称中点公式中点公式，相当于数值积分的中矩形公式：，相当于数值积分的中矩形公式：),(2,2(1nnnnnnyxfhyhxhfyy 如取如取a= ，则则c1=

19、c2= ， 2= 21 =1，即为改，即为改进进Euler公式。公式。若取若取a= 1，则则c1= 0，c2= 1， 2= 21 = ，得，得例例:蛇形曲线的初值问题蛇形曲线的初值问题令令f(x,y)=y/x 2y2, 取取f(0,0)=1, h=0.2, xn=nh , ( n = 1,2,15)2阶龙格阶龙格-库塔公式计算格式：库塔公式计算格式： k1=yn/xn 2yn2, k2 = (yn+hk1)/(xn+h) 2(yn+hk1)2 yn+1=yn + 0.5h k1 + k2 . 0)0()30(22yxyxyyx0=0;y0=0;h=.2;x=.2:h:3;k1=1;k2=(y0

20、+h*k1)/x(1)-2*(y0+h*k1)2;y(1)=y0+.5*h*(k1+k2); for n=1:14 k1=y(n)/x(n)-2*y(n)2; k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)2; y(n+1)=y(n)+0.5*h*(k1+k2);endy1=x./(1+x.2);subplot(221)plot(x,y1,x,y1,b*)subplot(222)plot(x,y,x,y,o)subplot(223)plot(x,y,x,y,o,x,y1,x,y1,b*)三、三阶与四阶三、三阶与四阶Runge-Kutta方法方法 ),(),(),()(2

21、32131331212213322111hKhKyhxfKhKyhxfKyxfKKcKcKchyynnnnnnnn 32313212321, ccc)(33221111KcKcKchyxyTnnn 类似二阶的推导过程，将类似二阶的推导过程，将K2, K3按二元函数展按二元函数展开，使开，使Tn+1= O(h4)，得得 61312113223233222332232313212321 cccccccc )2,()2,2(),()4(62131213211hKhKyhxfKKhyhxfKyxfKKKKhyynnnnnnnn ),()2,2()2,2(),()22(6342312143211hKyh

22、xfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn4、一阶常微分方程组和高阶微分方程的、一阶常微分方程组和高阶微分方程的 数值解法简介数值解法简介一阶常微分方程组的数值解法：一阶常微分方程组的数值解法： 00212002212210012111)(),()(),()(),(nnnnnnnyxyyyyxfyyxyyyyxfyyxyyyyxfy nyyyy21 ),(),(),(),(21yxfyxfyxfyxfn 00201002010)()()(nnyyyxyxyxyy 00)(),(yxyyxfy它在形式上跟单个微分方程的初值问题形式完它在形式上跟单个微分方程的初值

23、问题形式完全相同，只是函数变成了向量函数。故前面介全相同，只是函数变成了向量函数。故前面介绍的一切数值方法都适用，只要把函数换成向绍的一切数值方法都适用，只要把函数换成向量函数即可。量函数即可。 k1=f(xn,yn), k2=f(xn+0.5h,yn+0.5hk1)k3=f(xn+0.5h,yn+0.5hk2), k4=f(xn+h,yn+hk3) yn+1= yn+ hk1+2k2+2k3+k4/6(n = 0, 1, 2, ,N)四阶龙格四阶龙格-库塔公式库塔公式数值求解数值求解: 000)(),(yxyxxyxfy狐狸和野兔问题狐狸和野兔问题常微分方程组常微分方程组在一个生物圈中在一个

24、生物圈中,只有野兔和狐狸两种动物只有野兔和狐狸两种动物,记记y1 为野兔数量为野兔数量, y2 为狐狸数量为狐狸数量. 设有足够的食物供设有足够的食物供野兔生存野兔生存,而狐狸只以野兔为食物而狐狸只以野兔为食物. 模型如下模型如下 2122211102. 001. 0yyydxdyyyydxdy 自变量自变量x (0, 15),初值条件为初值条件为: y1(0)=20, y2(0)=20 21221102. 001. 0),(yyyyyyyxf 21yyy0510150100200300400t0=0;t1=15;%设置区域设置区域y0=20 20;%给定初值给定初值t,y=ode23(lot

25、1,t0,t1,y0);%求解求解plot(t,y)%绘图绘图function z=lot1(x,y)z(1,:)=y(1)-0.01*y(1).*y(2);z(2,:)=-y(2)+0.02*y(1).*y(2);求解常微分方程求解常微分方程:(1)定义函数定义函数;(2)调用调用ode235、高阶常微分方程的数值解法：、高阶常微分方程的数值解法： 0000)()(),(yxyyxyyyxfy引进新的变量引进新的变量,令令z=y，即可将上述二阶方程即可将上述二阶方程化为如下的一阶方程组的初值问题化为如下的一阶方程组的初值问题: 0000)(,)(),(yxyzyyxzzyxfz03 yyyy

26、1)0(, 1)0(,0)0( yyy)20( xyyyy 3yy 1yy2yy 31)0(, 1)0(, 0)0(332112333221yyyyyyyyyyy6 6、二阶常微分方程边值问题数值方法、二阶常微分方程边值问题数值方法考虑方程：考虑方程：bxayyxfy ),().(16;)(,)( byay).(26结合下述三种边界条件之一：结合下述三种边界条件之一：).(36).(46;)(,)( byay;)()(,)()(1010 bybyayay边界问题的解法：边界问题的解法：打靶法打靶法、有限差分法有限差分法(6.4)式中式中0, 0, 00000 它们分别称为第一、第二、第三边界问

27、题。它们分别称为第一、第二、第三边界问题。打靶法打靶法 基本思路：将基本思路：将边值边值问题转化为问题转化为初值初值问题考虑。或者问题考虑。或者说适当选择初始值使初值问题的解满足边值条件。然说适当选择初始值使初值问题的解满足边值条件。然后用求解后用求解初值初值问题的任一种有效的数值方法求解。问题的任一种有效的数值方法求解。以第一边界条件为例以第一边界条件为例考虑边值问题：考虑边值问题： )(,)(),(byaybxayyxfy取取y0= , 考虑初值问题考虑初值问题 00)(,)(),(yayyaybxayyxfy).( 56).(66)(),(000ygyyyynn 待定，由数值解法求解待定

28、，由数值解法求解(5.5)得到在得到在 处的处的解解 ，bxxn 0y ),(00yyyynn ,)(0时时 yg使使 ，这里，这里 为给定允许误差界，就为给定允许误差界，就停止迭代改进停止迭代改进, 输出作为数值解。输出作为数值解。 )(0yg0 0)(0 yg求解非线性方程求解非线性方程若若 , 则得所求数值解。当则得所求数值解。当 )(0yg该方程可用二分法、正割法或该方程可用二分法、正割法或Newton法等来求解。法等来求解。若求得若求得,0y 对第二类边界问题，也可转化为考虑初值问题对第二类边界问题，也可转化为考虑初值问题(6.5)，取取,0 y对第三类边界问题，仍可转化为考虑初值问

29、题对第三类边界问题，仍可转化为考虑初值问题(6.5)，取，取 0010yy 以以 为待定参数。为待定参数。0y，以，以 为待定参数。为待定参数。0y有限差分法有限差分法则离散化成差分方程则离散化成差分方程nihyyyxfhyyyiiiiiii, 2 , 1),2,(211211 (6.8),1 nabh将区间将区间a,b进行等分：进行等分：nihyyyiii, 2 , 1,211 (6.7)1, 1 , 0, nixxi设在设在, 1, 1 , 0, niihaxi处的数值解为处的数值解为 。iy用中心差分近似微分，即用中心差分近似微分，即 nihyyyyiiii, 2 , 1,2211 二阶中心差商二阶中心差商 对应的边界条件也离散成对应的边界条件也离散成 10,nyy(6.9) 第一边界问题：第一边界问题： hyyhyynn 101,(6.10) 第二边界问题：第二边界问题： 第三边界问题：第三边界问题：,)1 (1001hyhy (6.11)()()()()(),(xrxyxqxyxpyyxf 其中其中 为已知函数，则由常微分方程的理论知，通过为已知函数，则由常微分方程的理论知，通过)(),(),(xrxqxp )(,)()()()()(byayxrxyxqxy