sxjm53微分方程数值解ppt课件

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

2、11)(),()(),()4(yxyyyxfyyxyyyxfy本节主要研究问题本节主要研究问题(1)的数值解法，对的数值解法，对(2)(4)只作简单介绍。只作简单介绍。 00)(),()1(yxyyxfy考虑一阶常微分方程初值问题考虑一阶常微分方程初值问题其中，其中，y = y(x) 是未知函数，是未知函数，y(x0) = y0 是初是初值条件，而值条件，而f (x, y) 是给定的二元函数是给定的二元函数. 二二. 数值解的提法数值解的提法（其中（其中L为为Lipschitz常数则初值问常数则初值问题题1存在唯一的连续解。存在唯一的连续解。 2121),(),(yyLyxfyxf 若若f(x

3、)在在xa,b连续且连续且 f 满足对满足对 y 的的Lipschitz条件：条件：求问题求问题1的数值解，就是要寻找解函数的数值解，就是要寻找解函数在一系列离散节点在一系列离散节点x1 x2 xn xn+1 上的近似值上的近似值y1, y 2,yn 。为了计算方便，可取为了计算方便，可取xn=x0+nh，(n=0,1,2,)， h称为步长。称为步长。l常微分方程的数值解法有单步法和多常微分方程的数值解法有单步法和多步法之分：步法之分：l单步法：在计算单步法：在计算yn1 时只用到前一时只用到前一点点yn 的值的值 ；l多步法：计算多步法：计算yn1 时不仅利用时不仅利用yn，还要利用还要利用

4、yn-1, yn-2,.一般一般k步法要用到步法要用到 yn, yn-1, yn-2,., yn-k+1。1、欧拉、欧拉Euler方法方法在在x= x0 处，用差商代替导数：处，用差商代替导数：hxyxyxxxyxyxy)()()()()(0101010 由由00000)(),()(yxyyxfxy 10001),()(yyxhfyxy 得得三三. 简单的数值方法简单的数值方法同理，在同理，在x= xn 处，用差商代替导数：处，用差商代替导数：hxyxyxxxyxyxynnnnnnn)()()()()(111 ),()()(1nnnnyxhfxyxy ),()(nnnyxfxy 由由得得若记若

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

6、nnnyxyhyy 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. 计算结果如下：计算结果如下：xn y(xn) yn yn-y(xn)1.60.44940.49720.0478

7、1.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由表中数据可以看到，微分方程初值问题的数值解由表中数据可以看到，微分方程初值问题的数值解和解析解的误差一般在小数点后第二位或第三位小和解析解的误差一般在小数点后第二位或第三位小数上，这说明数上，这说明Euler方法的精度是比较差的。方法的精度是比较差的。 : 数值解数值解 : 准确解准确解 h=0.2;y

8、(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利用数值积分中的左矩形公式：利用数值积分中的左矩形公式：)(,()(,(1nnxxxyxfhdxxyxfnn ),(1nnnnyxhfyy 此即为此即为Euler公式。公式。设设y(xn)

9、= yn，则得，则得若用右矩形公式：若用右矩形公式：)(,()(,(111 nnxxxyxfhdxxyxfnn),(111 nnnnyxhfyy得得),()0(1nnnnyxhfyy ),()(11)1(1knnnknyxhfyy ),(),(11)(111)1(1 nnknnnknyxfyxfhyy1)(1 nknyyhL2、梯形方法、梯形方法dxxyxfxyxynnxxnn 1)(,()()(1由由 )(,()(,(2)(,(111 nnnnxxxyxfxyxfhdxxyxfnn ),(),(2111 nnnnnnyxfyxfhyy得得),()0(1nnnnyxhfyy ),(),(211

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

11、测校正校正上面两式统称预测校正法，上面两式统称预测校正法，又称改进的又称改进的Euler方法。方法。单步法的一般形式为单步法的一般形式为 ： （与与f 有关）有关）),(11hyyxhyynnnnn 显式单步法形式为：显式单步法形式为：),(1hyxhyynnnn 整体截断误差：从整体截断误差：从x0开始，考虑每一步产生的开始，考虑每一步产生的误差，直到误差，直到xn，则有误差，则有误差nnnyxye)(称为数值方法在节点称为数值方法在节点xn处的整体截断误差。处的整体截断误差。但但en不易分析和计算，故只考虑从不易分析和计算，故只考虑从xn到到xn+1的局的局部情况。部情况。四四. 单步法的

12、局部截断误差和精度单步法的局部截断误差和精度定义：设定义：设y(x)是初值问题是初值问题(1)的精确解，则称的精确解，则称),(,()()(11hxyxhxyxyTnnnnn 为显式单步法在节点为显式单步法在节点xn+1处的局部截断误差处的局部截断误差 。若存在最大整数若存在最大整数p使局部截断误差满足使局部截断误差满足)(),()()(11 pnhOhyxhxyhxyT 则称显式单步法具有则称显式单步法具有p阶精度或称阶精度或称p阶方法。阶方法。注：将注：将Tn+1表达式各项在表达式各项在xn处作处作Taylor展开，展开，可得具体表达式。可得具体表达式。Euler方法的局部截断误差：方法的

13、局部截断误差：)(,()()(11nnnnnxyxhfxyxyT )()()(nnnxyhxyhxy )()()(2)()(2nnnnnxyhxyyhxyhxy )()(2)(2322hOxyhyhnn ),(1 nnnxx )()2/(2nxyh 其中其中（设（设yn=y(xn)）)()(2)()(111 nnnnnxyxyhxyxyT)(! 3)(2)(32nnnxyhxyhxyh )()(2)()()(242hOxyhxyhxyxyhnnnn )()(1243hOxyhn 故故Tn+1= O(h3)，p=2，)()12/(3nxyh 1、Runge-Kutta方法的基本思想方法的基本思想

14、)(!)(2)()()()()(211nppnnnnnnxyphxyhxyhxyyhxyxy Tn+1= O(hp+1)，若提高，若提高p，可提高精度。，可提高精度。),(yxfy ),(),(),()(yxfyxfyxfxyyx 五五. RungeKutta方法方法 ),(111nnnnyxfKhKyy ),(),()22(121211hKyhxfKyxfKKKhyynnnnnnRunge-Kutta方法的基本思想方法的基本思想： ), 3 , 2(),(),(11111riKhyhxfKyxfKKchyyijjijnininnriiinn 上面第一个式子的右端在上面第一个式子的右端在(xn

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

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

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

18、Oyhyhy 即即 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= c2= ， 2=21 =1，即，即为改进为改进Eu

19、ler公式。公式。若取若取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+h*k1)/x(1)-2*(y0+h*k1)2;y

20、(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*)3、三阶与四阶、三阶与四阶Runge-Kutta方法方法 ),(),(),()(232131331212213322111hKhKyhx

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

22、yynnnnnnnnnn一阶常微分方程组的数值解法：一阶常微分方程组的数值解法： 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(lot1,t0,t1,y0);%求解求解plot(t,y)

25、%绘图绘图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)调用调用ode23 0000)()(),(yxyyxyyyxfy引进新的变量引进新的变量,令令z=y，即可将上述二阶方，即可将上述二阶方程化为如下的一阶方程组的初值问题程化为如下的一阶方程组的初值问题: 0000)(,)(),(yxyzyyxzzyxfz高阶常微分方程的数值解法：高阶常微分方程的数值解法：03 yyyy1)0(, 1)0(, 0)0( yyy)20( xyy

26、yy 3yy 1yy 2yy 3 1)0(, 1)0(, 0)0(332112333221yyyyyyyyyyy考虑方程：考虑方程：bxayyxfy ),()1.7(;)(,)( byay)2 . 7(结合下述三种边界条件之一：结合下述三种边界条件之一：)3 . 7()4 . 7(;)(,)( byay;)()(,)()(1010 bybyayay边界问题的解法：打靶法、有限差分法边界问题的解法：打靶法、有限差分法(7.4)式中式中0, 0, 00000 它们分别称为第一、第二、第三边界问题。它们分别称为第一、第二、第三边界问题。七七.二阶常微分方程边值问题数值方法二阶常微分方程边值问题数值方

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

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

29、为考虑初值问题(7.5)，取，取 0010yy 以以 为待定参数。为待定参数。0y，以，以 为待定参数。为待定参数。0y2 2、有限差分法、有限差分法则离散化成差分方程则离散化成差分方程nihyyyxfhyyyiiiiiii, 2 , 1),2,(211211 ,1 nabh将区间将区间a,b进行等分：进行等分：nihyyyiii, 2 , 1,211 1, 1 , 0, nixxi设在设在, 1, 1 , 0, niihaxi处的数值解为处的数值解为 。iy用中心差分近似微分，即用中心差分近似微分，即 nihyyyyiiii, 2 , 1,2211 二阶中心差二阶中心差商商)7 .7()8 .7(对应的边界条件也离散成对应的边界条件也离散成 10,nyy 第一边界问题：第一边界问题： hyyhyynn 101,第二边界问题：第二边界问题：第三边界问题：第三边界问题：,)1 (1001hyhy )()()()()(),(xrxyxqxyxpyyxf 其中其中 为已知函数，为已知函数，)(),(),(xrxqxp假设假设 是是 的线性函数时，的线性函数时，f 可写成可写成),(yyxf yy ,hyyhnn110)1 (