数学实验 常微分方程数值解

大学数学实验MathematicalExperiments实验4常微分方程数值解为什么要学习微分方程数值解•微分方程是研究函数变化规律的重要工具，有着广泛的应用。如:物体的运动,电路的电压,人口增长的预测•许多微分方程没有解析解，数值解法是求解的重要手段，如2,dyyxdx()()xtaxyytaxyby实验4的基本内容3.实际问题用微分方程建模，并求解2.龙格-库塔方法的MATLAB实现*4.数值算法的收敛性、稳定性与刚性方程1.两个最常用的数值算法:•欧拉（Euler）方法•龙格-库塔（Runge-Kutta）方法实例1海上缉私海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船正以速度a向正北方向行驶，缉私艇立即以最大速度b(a)前往拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终指向走私船。•建立任意时刻缉私艇位置及航线的数学模型,并求解;•求出缉私艇追上走私船的时间。a北bc艇船实例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,sindxdybbdtdt2222()()()()()()dxbcxdtcxatydybatydtcxaty由方程无法得到x(t),y(t)的解析解需要用数值解法求解“常微分方程初值问题数值解”的提法00'(,),()yfxyyxy设的解y=y(x)存在且唯一不求解析解y=y(x)(无解析解或求解困难)12nxxx而在一系列离散点()(1,2,)nyxnn求的近似值y通常取等步长h0nxxnh0x0y)(xyyyxy1y2yn)(nxy1x2xnxyP0x0x1x2x3xy=y(x)y0P1P2P3欧拉方法00'(,),()yfxyyxy基本思路1[,]nnxx在小区间上1(,)[,]nnfxyxx中的x取内的某一点1'[()()]/,nnyyxyxh11()()(,()),[,]nnnnyxyxhfxyxxxxx取不同点各种欧拉公式nxx取左端点1()()(,())nnnnyxyxhfxyx1(,),0,1,nnnnyyhfxyn向前欧拉公式显式公式11(),()nnnnyyxyyxP3欧拉方法11()()(,()),[,]nnnnyxyxhfxyxxxx11(),()nnnnxyyxyyx取右端点，111(,),0,1,nnnnyyhfxyn向后欧拉公式隐式公式yP0x0x1x2x3xy0y=y(x)P1P2n+1y右端未知，需迭代求解(0)1(,)nnnnyyhfxy初值(1)()111(,)0,1,2,,0,1,2,kknnnnyyhfxykn()11limknnkyy向前欧拉公式向后欧拉公式1(,)nnnnyyhfxy111(,)nnnnyyhfxy二者平均得到梯形公式111[(,)(,)],0,1,2nnnnnnhyyfxyfxyn仍为隐式公式，需迭代求解将梯形公式的迭代过程简化为两步1(,)nnnnyyhfxy预测111[(,)(,)],0,1,2nnnnnnhyyfxyfxyn校正改进欧拉公式龙格-库塔方法11()()(,()),[,]nnnnyxyxhfxyxxxx•向前，向后欧拉公式：用[xn,xn+1]内１个点的导数代替f(x,y(x))•梯形公式，改进欧拉公式：用[xn,xn+1]内２个点导数的平均值代替f(x,y(x))龙格-库塔方法的基本思想在[xn,xn+1]内多取几个点，将它们的导数加权平均代替f(x,y(x))，设法构造出精度更高的计算公式。1,10,1,,111ijijiLiiijiiacac满足常用的(经典)龙格—库塔公式不足:收敛速度较慢111222111(,)(,)(,),3,4,,LnniiinnnniininiijjjyyhkkfxykfxchychkkfxchychakiLL级?阶龙格-库塔方法的一般形式112341213243(22)6(,)(/2,/2)(/2,/2)(,)nnnnnnnnnnhyykkkkkfxykfxhyhkkfxhyhkkfxhyhk使精度尽量高4级4阶微分方程组和高阶方程初值问题的数值解欧拉方法和龙格-库塔方法可直接推广到微分方程组0000(,,)(,,)(),()yfxyzzgxyzyxyzxz向前欧拉公式),,(yyxfy1yy),,(21221yyxfyyy高阶方程需要先降阶为一阶微分方程组,2,1,0),,(),,(11nzyxhgzzzyxhfyynnnnnnnnnn龙格—库塔方法的MATLAB实现TnTnfffxxxxtxxtftx),(,),(,)(),,()(1100[t,x]=ode23(@f,ts,x0,opt)3级2阶龙格-库塔公式[t,x]=ode45(@f,ts,x0,opt)5级4阶龙格-库塔公式f是待解方程写成的函数m文件：functiondx=f(t,x)dx=[f1;f2;;fn];ts=[t0,t1,…,tf]输出指定时刻t0,t1,…,tf的函数值ts=t0:k:tf输出[t0,tf]内等分点处的函数值x0为函数初值(n维)输出t=ts,x为相应函数值(n维)opt为选项，缺省时精度为：相对误差10-3，绝对误差10-6,计算步长按精度要求自动调整实例1海上缉私(续)模型的数值解2222()()()()()()dxbcxdtcxatydybatydtcxaty(0)0,(0)0xy0yxc(x(t),y(t))ab设:船速a=20(海里/小时)艇速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=20tt=0.5时缉私艇追上走私船051015200246810xy缉私艇的航线y(x)tx(t)y(t)0000.051.99840.06980.103.98540.29240.155.94450.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=?缉私艇追上走私船tx(t)y(t)y1(t)00000.13.95610.50583.50.27.59282.13087.00.310.52404.828310.50.412.53848.275514.00.513.755112.083017.51.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('RelTol',1e-6,'AbsTol',1e-9);[t,x]=ode45(@jisi,ts,x0,opt);t=1.6时缉私艇追上走私船051015200204060xy缉私艇的航线y(x)判断“追上”的有效方法?tx(t)y(t)y1(t)00000.13.9561040.5058133.50.27.5928222.1306787.00.310.5219214.82930810.50.412.5394548.26984014.00.513.75397412.07534417.51.214.99961640.00000542.01.314.99996344.00000545.51.414.99999348.00000549.01.514.99999852.00000552.51.615.00002055.99993156.0实例1海上缉私(续)模型的解析解2222)()()(,)()()(yatxcyatbdtdyyatxcxcbdtdx)(xcyatdxdyatydxdyxc)(22()dydtcxadxdxdsbdt22()()dsdxdy222()1()dyadycxdxbdxdxdyp令bakpkdxdpxc,1)(2(0)0pdydsdsdtdxdtdxdtdtdydxdy实例1海上缉私(续)模型的解析解21)(pkdxdpxc(0)0p21()kcxppc21()kcxppc1[()()]2kkcxcxpcc/1kab(0)0y11211[()()]2111kkccxcxkcykckck缉私艇的航线y(x)的解析解x=c时缉私艇追上走私船的y坐标缉私艇追上走私船的时间:122bctbaa=20,b=40,c=15t1=0.5a=35,b=40,c=15t1=1.62221kcabcykba微分方程数值解法的误差分析数值解法:计算微分方程精确解y(xn)的近似值yn按照步长h一步步计算,每步都有误差;每一步的误差会逐步积累,称累积误差.讨论计算一步出现的误差假定yn=y(xn),估计yn+1的误差:y(xn+1)-yn+1局部截断误差误差分析估计欧拉公式的局部截断误差y(xn+1)在xn处作Taylor展开:)()(2)()()(321hOxyhxyhxyxynnnn向前欧拉公式1(,)nnnnyyhfxyyn=y(xn))()())(,()(1nnnnnnxyhxyxyxhfxyy232111()()()()2nnnnhTyxyyxOhOh局部截断误差主项为)(22nxyh误差分析估计欧拉公式的局部截断误差)()(2)()()(321hOxyhxyhxyxynnnn向后欧拉公式),(111nnnnyxhfyy232111()()()()2nnnnhTyxyyxOhOh局部截断误差主项为)(22nxyh梯形公式向前、向后欧拉公式的平均3431()()()12nnhTyxOhOhxnxn+1xyPnAQB误差分析算法精度的阶的定义一个算法的局部截断误差为O(hp+1)该算法具有p阶精度局部截断误差精度向前欧拉公式O(h2)1阶向后欧拉公式O(h2)1阶梯形公式O(h3)2阶改进欧拉公式O(h3)2阶经典龙格-库塔公式O(h5)4阶实例2弱肉强食问题自然界中同一环境下两个种群之间的生存方式相互竞争相互依存弱肉强食弱肉强食种群甲靠丰富的自然资源生存食饵(Prey)种群乙靠捕食种群甲为生捕食者(Predator)两个种群的数量如何演变?•历史背景——一次世界大战期间地中海渔业的捕捞量下降(食用鱼和鲨鱼同时捕捞)，但是其中鲨鱼的比例却增加，为什么？实例2弱肉强食模型食饵(甲)的密度x(t),捕食者(乙)的密度y(t)0,/rrxx甲独立生存的增长率r0,/aay