常微分方程数值解

1、大学数学实验大学数学实验Mathematical Experiments 实验实验4 常微分方程数值解常微分方程数值解为什么要学习微分方程数值解为什么要学习微分方程数值解 微分方程是研究函数变化规律的重要工具，有着广泛微分方程是研究函数变化规律的重要工具，有着广泛的应用。如的应用。如: 物体的运动物体的运动, 电路的电压电路的电压, 人口增长的预测人口增长的预测 许多微分方程没有解析解，数值解法是求解的重要手许多微分方程没有解析解，数值解法是求解的重要手段，如段，如2,dyyxdx( )( )x taxyy taxyby 实验实验4 4的基本内容的基本内容3. 实际问题用微分方程建模，并求解实

2、际问题用微分方程建模，并求解2. 龙格龙格-库塔方法的库塔方法的MATLAB实现实现*4. 数值算法的收敛性、稳定性与刚性方程数值算法的收敛性、稳定性与刚性方程1. 两个最常用的数值算法两个最常用的数值算法: 欧拉（欧拉（Euler）方法）方法 龙格龙格-库塔（库塔（Runge-Kutta）方法）方法实例实例1 1 海上缉私海上缉私 海防某部缉私艇上的雷达发现正东方向海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船海里处有一艘走私船正以速度正以速度a向正北方向行驶，缉私艇立即以最大速度向正北方向行驶，缉私艇立即以最大速度b(a)前往前往拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始

3、终拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终指向走私船。指向走私船。 建立任意时刻缉私艇位置及建立任意时刻缉私艇位置及 航线的数学模型航线的数学模型,并求解并求解; 求出缉私艇追上走私船的时间。求出缉私艇追上走私船的时间。 a北北bc艇艇船船实例实例1 1 海上缉私海上缉私 建立坐标系如图建立坐标系如图: t=0 艇在艇在(0, 0), 船在船在(c, 0); 船速船速a, 艇速艇速b 时刻时刻 t 艇位于艇位于P(x, y), 船到达船到达 Q(c, at)模型模型:0yxcR(c,y ) Q(c,at)P(x,y)bcos ,sindxdybbdtdt2222()()()()(

4、)()dxb cxdtcxatydyb atydtcxaty由方程无法得到由方程无法得到x(t), y(t)的解析解的解析解需要用数值解法求解需要用数值解法求解 “常微分方程初值问题数值解常微分方程初值问题数值解”的提法的提法00( , ), ()yf x yy xy设的解y=y(x)存在且唯一不求解析解不求解析解 y = y(x) (无解析解或求解困难无解析解或求解困难)12nxxx而在一系列离散点而在一系列离散点()(1, 2 ,)nyxn n求的 近 似 值 y通常取等步长通常取等步长h0nxxnh0 x0y)(xyy yxy1y2yn)(nxy1x2xnxyP0 x0 x1x2x3 x

5、y=y(x)y0P1P2P3欧 拉 方 法00( , ), ()yf x y y xy基本基本思路思路1,nnx x在小区间上1( , ),nnf x yx x中的x取内的某一点1 ()()/ ,nnyy xy xh11()( )( , ( ), ,nnnny xy xhf x y xxx xx取不同点取不同点各种欧拉公式各种欧拉公式nxx取左端点1()()( , ()nnnny xy xhf x y x1( ,),0,1,nnnnyyhf x yn向前欧拉公式向前欧拉公式显式公式显式公式11(),()nnnnyy xyy xP3欧 拉 方 法11()()( , ( ),nnnny xy xh

6、f x y xxx x11(),()nnnnxyy xyy x取右端点，111(,),0,1,nnnnyyhf xyn向后欧拉公式向后欧拉公式 隐式公式隐式公式yP0 x0 x1x2x3 xy0y=y(x)P1P2n+1y右端未知，需迭代求解(0)1(,)nnnnyyhfxy初值(1)()111(,)0,1,2,0,1,2,kknnnnyyhfxykn()11limknnkyy 向前欧拉公式向前欧拉公式向后欧拉公式向后欧拉公式1(,)nnnnyyhf xy111(,)nnnnyyhfxy二者平均得到二者平均得到梯形公式梯形公式111(,)(,),0,1,2nnnnnnhyyfxyfxyn仍为隐

7、式公式，需迭代求解仍为隐式公式，需迭代求解将梯形公式的迭代过程简化为两步将梯形公式的迭代过程简化为两步1(,)nnnnyyhf x y预测预测111(,)(,),0,1,2nnnnnnhyyfxyfxyn校正校正改进欧拉公式改进欧拉公式龙格龙格- -库塔方法库塔方法11()()( , ( ),nnnny xy xhf x y xxx x 向前，向后欧拉公式：向前，向后欧拉公式： 用用xn, xn+1内个点的导内个点的导 数代替数代替 f(x, y(x) 梯形公式，改进欧拉公式：梯形公式，改进欧拉公式：用用xn, xn+1内个点导数的平均值代替内个点导数的平均值代替 f(x, y(x)龙格龙格-

8、 -库塔方法的基本思想库塔方法的基本思想在在xn, xn+1内多取几个点，将它们的导数加权平内多取几个点，将它们的导数加权平均代替均代替 f(x, y(x)，设法构造出设法构造出精度更高精度更高的计算公式。的计算公式。 1, 10, 1,111ijijiLiiijiiacac满足常用的常用的(经典经典)龙格龙格库塔库塔公式公式不足不足:收敛速度较慢收敛速度较慢111222111(,)(,)(,),3,4,Lnniiinnnniininiijjjyyhkkf xykf xc h yc hkkf xc h yc ha kiLL级?阶龙格龙格-库塔库塔方法的方法的一般形式一般形式1123412132

9、43(22)6(,)(/2,/2)(/2,/2)(,)nnnnnnnnnnhyykkkkkf xykf xhyhkkf xhyhkkf xh yhk使精度尽量高使精度尽量高4级4阶微分方程组和高阶方程初值问题的数值解微分方程组和高阶方程初值问题的数值解 欧拉方法和龙格欧拉方法和龙格-库塔方法可直接推广到微分方程组库塔方法可直接推广到微分方程组0000( , , )( , , )(), ()yf x y zzg x y zy xyz xz向前欧拉公式向前欧拉公式 ),(yyxfy 1yy),(21221yyxfyyy高阶方程需要先降阶为一阶微分方程组高阶方程需要先降阶为一阶微分方程组 , 2 ,

10、 1 , 0),(),(11nzyxhgzzzyxhfyynnnnnnnnnn龙格龙格库塔方法的库塔方法的 MATLAB 实现实现 TnTnfffxxxxtxxtftx),(,),(,)(),()(1100t,x=ode23(f,ts,x0,opt) 3级级2阶龙格阶龙格-库塔公式库塔公式 t,x=ode45(f,ts,x0,opt) 5级级4阶龙格阶龙格-库塔公式库塔公式 f是待解方程写成是待解方程写成的函数的函数m文件：文件： function dx=f(t,x)dx=f1; f2; fn;ts = t0,t1, ,tf输出指定时刻输出指定时刻 t0,t1, ,tf 的函数值的函数值ts

11、= t0:k:tf输出输出 t0,tf 内等分点处的函数值内等分点处的函数值x0为函数初值为函数初值(n维维) 输出输出t=ts, x为相应函数值为相应函数值(n维维) optopt为选项，缺省时精度为：相对误差为选项，缺省时精度为：相对误差1010-3-3，绝对误差绝对误差1010-6-6, , 计算步长按精度要求自动调整计算步长按精度要求自动调整实例实例1海上缉私海上缉私( (续续) ) 模型的数值解模型的数值解2222()()()()()()dxb cxdtcxatydyb atydtcxaty(0)0,(0)0 xy0yxc(x(t),y(t)ab设设: 船速船速a=20 (海里海里/

12、小时小时) 艇速艇速b=40 (海里海里/小时小时) 距离距离c=15 (海里海里) 求求: 缉私艇的位置缉私艇的位置x(t),y(t) 缉私艇的航线缉私艇的航线y(x) jisi.m, seajisi.m 实例实例1海上缉私海上缉私( (续续) ) 模型的数值解模型的数值解 a=20, b=40, c=15走私船的位置走私船的位置x1(t)= c=15y1(t)=at=20t t=0.5时缉私艇追上走私船时缉私艇追上走私船 051015200246810 xy缉私艇的航线缉私艇的航线y(x)tx(t)y(t)0000.051.99840.06980.103.98540.29240.155.9

13、4450.69060.207.85151.28990.259.67052.11780.3011.34963.20050.3512.81704.55520.4013.98066.17730.4514.74518.02730.5015.00469.9979y1(t)01.02.03.04.05.06.07.08.09.010.0实例实例1 海上缉私海上缉私( (续续) ) 模型的数值解模型的数值解 设设b，c不变，不变，a变大变大为为30, 35, 接近接近40, 观察解的变化观察解的变化 :a=35, b=40, c=15t=? 缉私艇追上走私船缉私艇追上走私船 t x(t) y(t) y1(t

14、)0 0 0 00.13.95610.50583.50.27.59282.13087.00.310.52404.828310.50.412.53848.275514.00.513.755112.083017.51.214.998640.016442.01.314.999644.016545.51.415.011748.018349.01.515.002352.014652.51.614.986655.948656.0累积误差较大累积误差较大提高精度提高精度! ! 实例实例1海上缉私海上缉私( (续续) ) 模型的数值解模型的数值解 a=35, b=40, c=15opt=odeset(RelT

15、ol,1e-6, AbsTol,1e-9);t,x=ode45(jisi,ts,x0,opt);t=1.6时缉私艇追上走私船时缉私艇追上走私船 051015200204060 xy缉私艇的航线缉私艇的航线y(x)判断判断“追上追上”的有效方的有效方法法?t x(t) y(t)y1(t)0 0 0 00.13.9561040.5058133.50.27.5928222.1306787.00.310.5219214.82930810.50.412.5394548.26984014.00.513.75397412.07534417.51.214.99961640.00000542.01.314.99

16、996344.00000545.51.414.99999348.00000549.01.514.99999852.00000552.51.615.00002055.99993156.0实例实例1 海上缉私海上缉私( (续续) ) 模型的解析解模型的解析解 2222)()()(,)()()(yatxcyatbdtdyyatxcxcbdtdx)(xcyatdxdyatydxdyxc )(22()d ydtcxadxdxdsbdt22()()dsdxdy222()1 ()d yadycxdxbdxdxdyp 令bakpkdxdpxc,1)(2(0)0pdydsdsdtdxdtdxdtdtdydxdy

17、实例实例1海上缉私海上缉私( (续续) ) 模型的解析解模型的解析解 21)(pkdxdpxc(0)0p21()kcxppc21()kcxppc1()() 2kkcxcxpcc/1ka b(0)0y11211()()2 111kkccxcxkcykckck缉私艇的航线缉私艇的航线y y( (x x) )的解析解的解析解x=c时时 缉私艇追上走私船的缉私艇追上走私船的y y坐标坐标 缉私艇追上走私船的时间缉私艇追上走私船的时间: 122bctbaa=20, b=40, c=15 t1=0.5 a=35, b=40, c=15 t1=1.6 2221kcabcykba微分方程数值解法的误差分微分方

18、程数值解法的误差分析析数值解法数值解法: 计算微分方程精确解计算微分方程精确解y(xn)的近似值的近似值yn按照步长按照步长h一步步计算一步步计算, 每步都有误差每步都有误差;每一步的误差会逐步积累每一步的误差会逐步积累, , 称称累积误差累积误差. .讨论计算一步出现的误差讨论计算一步出现的误差假定假定yn= y(xn) , 估计估计yn+1的误差的误差: y(xn+1)- yn+1局部截断误差局部截断误差 误差分析误差分析估计欧拉公式的局部截断误差估计欧拉公式的局部截断误差 y(xn+1)在在xn处作处作Taylor展开展开: )()(2)()()(321hOxyhxyhxyxynnnn

19、向前欧拉公式向前欧拉公式1(,)nnnnyyhf xyyn= y(xn)()()(,()(1nnnnnnxyhxyxyxhfxyy232111()()()()2nnnnhTy xyyxO hO h局部截断误差主项为局部截断误差主项为误差分析误差分析估计欧拉公式的局部截断误差估计欧拉公式的局部截断误差 )()(2)()()(321hOxyhxyhxyxynnnn 向后欧拉公式向后欧拉公式),(111nnnnyxhfyy232111()()()()2nnnnhTy xyyxO hO h 局部截断误差主项为局部截断误差主项为)(22nxyh 梯形梯形公式公式向前、向后欧拉公式的平均向前、向后欧拉公式

20、的平均3431()()()12nnhTyxO hO h xn xn+1 x y Pn A Q B误误差差分分析析算法算法精度的阶精度的阶的定义的定义 一个算法的局部截断误差为一个算法的局部截断误差为O(hp+1) 该算法具有该算法具有p阶精度阶精度 局部截断误差局部截断误差精度精度向前欧拉公式向前欧拉公式O(h2)1阶阶向后欧拉公式向后欧拉公式O(h2)1阶阶梯形公式梯形公式O(h3)2阶阶改进欧拉公式改进欧拉公式O(h3)2阶阶经典龙格经典龙格- -库塔公式库塔公式O(h5)4阶阶实实 例例 2弱弱 肉肉 强强 食食 问问题题自然界中同一环境下两个种群之间的生存方式自然界中同一环境下两个种群

21、之间的生存方式 相互竞争相互竞争 相互依存相互依存 弱肉强食弱肉强食 弱肉弱肉强食强食种群甲靠丰富的自然资源生存种群甲靠丰富的自然资源生存 食饵食饵(Prey) (Prey) 种群乙靠捕食种群甲为生种群乙靠捕食种群甲为生 捕食者捕食者(Predator) (Predator) 两个种群的数量如何演变两个种群的数量如何演变? ? 历史背景历史背景一次世界大战期间地中海渔业一次世界大战期间地中海渔业的捕捞量下降的捕捞量下降(食用鱼和鲨鱼同时捕捞食用鱼和鲨鱼同时捕捞)，但，但是其中是其中鲨鱼的比例却增加，为什么？鲨鱼的比例却增加，为什么？实实 例例 2弱弱 肉肉 强强 食食 模型模型食饵食饵(甲甲)

22、 的密度的密度x(t), 捕食者捕食者(乙乙)的密度的密度y(t)0,/rrxx 甲独立生存的增长率甲独立生存的增长率r r0,/aayrxx 乙使甲的增长率减小乙使甲的增长率减小, ,减小量与减小量与 y y 成正比成正比0,/ddyy 乙独立生存的死亡率乙独立生存的死亡率d0),(/bbxdyy 甲使乙的死亡率减小甲使乙的死亡率减小, ,减小量与减小量与 x成正比成正比00()()(0),(0)xray xrxaxyydbx ydybxyxxyyVolterra模型模型 x(t), y(t)无解析解无解析解实例实例2弱肉强食弱肉强食 模型的数值解模型的数值解 00)0(,)0(yyxxbx

23、ydyyaxyrxx001,0.50.1,0.0225,2rdabxy051015020406080100 x(t)y(t)020406080100051015202530 xy猜测猜测 x(t), y(t)是周期函数是周期函数; y(x)是封闭曲线是封闭曲线 数值积分计算一个周期的平均值：数值积分计算一个周期的平均值： 25,10 xyshier.m, shier1.mybxdxayrdydx)()(实例实例2弱肉强食弱肉强食 模型的解析解模型的解析解 ybxdyxayrx)()(ceyexayrbxd)(相轨线dyyayrdxxbxdc由初始条件确定由初始条件确定 相轨线是封闭曲线相轨线是

24、封闭曲线(c在一定范围内在一定范围内)求求x(t),y(t)一周期的平均值一周期的平均值:yx ,可以可以证明证明ybxdty)()(bdyytx/ )/()(TdttxTx0)(1bdx x(t), y(t)是周期函数是周期函数(周期记作周期记作T) )0(ln)(ln(1bdTbyTyT实例实例2弱肉强食弱肉强食 模型的解析解模型的解析解 ryax(t),y(t)一周期的平均值一周期的平均值:bdx 02. 0, 1 . 0, 5 . 0, 1badr10,25yxr 食饵增长率食饵增长率a 捕食者对食饵的捕获能力捕食者对食饵的捕获能力 d 捕食者死亡率捕食者死亡率b 食饵对捕食者的喂养能

25、力食饵对捕食者的喂养能力 结果解释结果解释ybxdyxayrx)()(与计算结果同与计算结果同yar ,xbd ,既相互制约既相互制约又相互依存又相互依存02040608010012005101520253000yx，00yx，00yx00yx0PT2T3T4T1P024681012020406080100120 x(t)y(t)T1 T2 T3 T4x(t) 的的“相位相位”领先领先 y(t)xayrtx)()(ybxdty)()()()(:1tytxT)()(:2tytxT)()(:3tytxT)()(:4tytxT进一步分析进一步分析)/,/(arbdP),(000yxP初值初值相轨线的

26、方向相轨线的方向一次大战期间地中海渔业的捕捞量下降，一次大战期间地中海渔业的捕捞量下降，但是其中但是其中鲨鱼的比例却在增加，为什么？鲨鱼的比例却在增加，为什么？rr- 1, dd+ 1捕捞捕捞战时战时捕捞捕捞rr- 2, dd+ 2 , 2 0微分方程稳定微分方程稳定 ),(1nnnnyxhfyy向前欧拉公式向前欧拉公式向后欧拉公式向后欧拉公式11nnnyyh y经典经典龙格龙格- -库塔公式库塔公式2.785/h算法稳定算法稳定)(,()(,(),(*yyyxfxxyxfyxfyyx0,yy(特征根特征根 - )nnh)1 (1nkn11h2/h111nnhh任意任意xycenyh )1 ( 刚性现象与刚性方刚性现象与刚性方程程bxaxrktfrxxkx) 0 (,) 0 (),0,()( 振动系统或电路系统的数学模型振动系统或电路系统的数学模型 现象现象 k=2000.5, r=1000, a=1, b=-1999.5, f(t)=1 瞬态解与稳态解瞬态解与稳态解 e-2000t 快瞬态解快瞬态解 e