实验4常微分方程数值解

1,实验4常微分方程数值解,凯里学院理学院,主讲:潘东云,ExperimentsinMathematics,2,为什么要学习微分方程数值解,微分方程是研究函数变化规律的重要工具，有着广泛的应用。如:物体的运动,电路的电压,人口增长的预测许多微分方程没有解析解，数值解法是求解的重要手段，如,3,实验4的基本内容,3.实际问题用微分方程建模，并求解,2.龙格-库塔方法的MATLAB实现,*4.数值算法的收敛性、稳定性与刚性方程,两个最常用的数值算法:,欧拉（Euler）方法龙格-库塔（Runge-Kutta）方法,4,实例1海上缉私,海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船正以速度a向正北方向行驶，缉私艇立即以最大速度b(a)前往拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终指向走私船。,建立任意时刻缉私艇位置及航线的数学模型,并求解;求出缉私艇追上走私船的时间。,5,实例1海上缉私,建立坐标系如图:t=0艇在(0,0),船在(c,0);船速a,艇速b时刻t艇位于P(x,y),船到达Q(c,at),模型:,由方程无法得到x(t),y(t)的解析解需要用数值解法求解,6,“常微分方程初值问题数值解”的提法,不求解析解y=y(x)(无解析解或求解困难),而在一系列离散点,通常取等步长h,7,欧拉方法,基本思路,x取不同点,各种欧拉公式,向前欧拉公式,显式公式,8,欧拉方法,向后欧拉公式隐式公式,9,向前欧拉公式,向后欧拉公式,二者平均得到梯形公式,仍为隐式公式，需迭代求解,将梯形公式的迭代过程简化为两步,预测,校正,改进欧拉公式,10,龙格-库塔方法,向前，向后欧拉公式：用xn,xn+1内个点的导数代替f(x,y(x),梯形公式，改进欧拉公式：用xn,xn+1内个点导数的平均值代替f(x,y(x),龙格-库塔方法的基本思想,在xn,xn+1内多取几个点，将它们的导数加权平均代替f(x,y(x)，设法构造出精度更高的计算公式。,11,常用的(经典)龙格库塔公式,不足:收敛速度较慢,L级?阶,龙格-库塔方法的一般形式,使精度尽量高,4级4阶,12,微分方程组和高阶方程初值问题的数值解,欧拉方法和龙格-库塔方法可直接推广到微分方程组,高阶方程需要先降阶为一阶微分方程组,13,龙格库塔方法的MATLAB实现,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,计算步长按精度要求自动调整,14,微分方程的符号解,注意:在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分.任何D后所跟的字母为因变量，自变量可以指定或由系统规则选定为确省.例如，微分方程,应表达为：D2y=0.,结果：u=tan(t-c),15,解输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x),结果为:y=3e-2xsin（5x）,16,解输入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t)；x=simple(x)%将x化简y=simple(y)z=simple(z),结果为：x=(c1-c2+c3+c2e-3t-c3e-3t)e2ty=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t,17,解：用MATLAB编程如下：,functiondx=fun113(t,x)dx=t+x;,ts=0:0.1:1;x0=1;t,x=ode45(fun113,ts,x0);y=2*exp(t)-t-1;t,x,y,x-yplot(t,x,r-,t,y,b-.),grid,gtext(x(t),gtext(y(t),title(微分方程数值解例);xlabel(timet);ylabel(solutionx);legend(x,y);,微分方程的数值解,18,微分方程的数值解,解:令y1=x，y2=y1,1、建立m-文件fun114.m如下：functiondx=fun114(t,x)dx=zeros(2,1);dx(1)=x(2);dx(2)=-x(2)/t+(1-0.25/t2)*x(1);,2、取初值，输入命令：t,x=ode45(fun114,pi/22,pi/2-2/piplot(t,x(:,1),-),19,实例1海上缉私(续),模型的数值解,设:船速a=20(海里/小时)艇速b=40(海里/小时)距离c=15(海里),求:缉私艇的位置x(t),y(t)缉私艇的航线y(x),jisi.m,seajisi.m,20,实例1海上缉私(续),模型的数值解,a=20,b=40,c=15,走私船的位置x1(t)=c=15y1(t)=at=20t,t=0.5时缉私艇追上走私船,缉私艇的航线y(x),21,实例1海上缉私(续),模型的数值解,设b，c不变，a变大为30,35,接近40,观察解的变化:,a=35,b=40,c=15,t=?缉私艇追上走私船,累积误差较大,提高精度!,22,实例1海上缉私(续),模型的数值解,a=35,b=40,c=15,opt=odeset(RelTol,1e-6,AbsTol,1e-9);t,x=ode45(jisi,ts,x0,opt);,t=1.6时缉私艇追上走私船,缉私艇的航线y(x),判断“追上”的有效方法?,23,实例1海上缉私(续),模型的解析解,24,实例1海上缉私(续),模型的解析解,缉私艇的航线y(x)的解析解,x=c时,缉私艇追上走私船的y坐标,缉私艇追上走私船的时间:,a=20,b=40,c=15t1=0.5,a=35,b=40,c=15t1=1.6,25,微分方程数值解法的误差分析,数值解法:计算微分方程精确解y(xn)的近似值yn,按照步长h一步步计算,每步都有误差;,每一步的误差会逐步积累,称累积误差.,讨论计算一步出现的误差,假定yn=y(xn),估计yn+1的误差:y(xn+1)-yn+1,局部截断误差,26,误差分析,估计欧拉公式的局部截断误差,y(xn+1)在xn处作Taylor展开:,向前欧拉公式,yn=y(xn),局部截断误差主项为,27,误差分析,估计欧拉公式的局部截断误差,向后欧拉公式,局部截断误差主项为,梯形公式,向前、向后欧拉公式的平均,28,误差分析,算法精度的阶的定义,一个算法的局部截断误差为O(hp+1),局部截断误差,精度,向前欧拉公式,O(h2),1阶,向后欧拉公式,O(h2),1阶,梯形公式,O(h3),2阶,改进欧拉公式,O(h3),2阶,经典龙格-库塔公式,O(h5),4阶,29,实例2弱肉强食,问题,自然界中同一环境下两个种群之间的生存方式,相互竞争,相互依存,弱肉强食,弱肉强食,种群甲靠丰富的自然资源生存,食饵(Prey),种群乙靠捕食种群甲为生,捕食者(Predator),两个种群的数量如何演变?,历史背景一次世界大战期间地中海渔业的捕捞量下降(食用鱼和鲨鱼同时捕捞)，但是其中鲨鱼的比例却增加，为什么？,30,实例2弱肉强食,模型,食饵(甲)的密度x(t),捕食者(乙)的密度y(t),甲独立生存的增长率r,乙使甲的增长率减小,减小量与y成正比,乙独立生存的死亡率d,甲使乙的死亡率减小,减小量与x成正比,Volterra模型,x(t),y(t)无解析解,31,实例2弱肉强食,模型的数值解,猜测,x(t),y(t)是周期函数;,y(x)是封闭曲线,数值积分计算一个周期的平均值：,shier.m,shier1.m,32,实例2弱肉强食,模型的解析解,c由初始条件确定,相轨线是封闭曲线(c在一定范围内),求x(t),y(t)一周期的平均值:,可以证明,33,实例2弱肉强食,模型的解析解,x(t),y(t)一周期的平均值:,r食饵增长率,a捕食者对食饵的捕获能力,d捕食者死亡率,b食饵对捕食者的喂养能力,结果解释,与计算结果同,既相互制约又相互依存,34,x(t)的“相位”领先y(t),进一步分析,初值,相轨线的方向,35,一次大战期间地中海渔业的捕捞量下降，但是其中鲨鱼的比例却在增加，为什么？,rr-1,dd+1,捕捞,战时捕捞,rr-2,dd+2,20微分方程稳定,向前欧拉公式,向后欧拉公式,经典龙格-库塔公式,40,刚性现象与刚性方程,振动系统或电路系统的数学模型,现象,k=2000.5,r=1000,a=1,b=-1999.5,f(t)=1,瞬态解与稳态解,e-2000t快瞬态解,e-t/2慢瞬态解,计算到t=0.005时已衰减到4.510-5,计算到t=20时才衰减到4.510-5,精度达到10-4需算到t=20(由慢瞬态解=1/2决定),选取步长h由快瞬态解=2000决定,h2.785/2000=0.0014,龙格-库塔公式,t=20需14286步,快、慢瞬态解的特征根相差悬殊,求稳态解,41,刚性现象与刚性方程,刚性方程,振动、电路及化学反应中的线性常系数方程组,A的特征根k(k=1,2,n)的实部Re(k)10,刚性非线性方程组,线性常系数方程组,42,刚性现象与刚性方