常微分方程数值解课件

*,1,常微分方程数值解,1常微分方程数值解,2,实例,1,海上缉私,海防某部缉私艇上的雷达发现正东方向,c,海里处有一艘走私船正以速度,a,向正北方向行驶，缉私艇立即以最大速度,b(a),前往拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终指向走私船。,建立任意时刻缉私艇位置及,航线的数学模型,并求解,;,求出缉私艇追上走私船的时间。,a,北,b,c,艇,船,2实例1 海上缉私 海防某部缉私艇上的雷达发现正东方向c海,3,实例,1,海上缉私,建立坐标系如图,:,t,=0,艇在,(0,0),船在,(c,0);,船速,a,艇速,b,时刻,t,艇位于,P,(,x,y,),船到达,Q,(,c,at,),模型,:,0,y,x,c,R,(,c,y,),Q,(,c,at,),P,(,x,y,),b,由方程无法得到,x,(,t,),y,(,t,),的解析解,需要用数值解法求解,3实例1 海上缉私 建立坐标系如图:t=0 艇在(0,4,龙格,库塔方法的,MATLAB,实现,t,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=t0:k:tf,输出,t0,tf,内等分点处的函数值,x0,为函数初值,(n,维,),输出,t=ts,x,为相应函数值,(n,维,),opt,为选项，缺省时精度为：相对误差,10,-3,，绝对误差,10,-6,计算步长按精度要求自动调整,4龙格库塔方法的 MATLAB 实现 t,x=ode2,5,实例,1,海上缉私,(,续,),模型的数值解,0,y,x,c,(,x,(,t,),y,(,t,),),a,b,设,:,船速,a,=20(,海里,/,小时,),艇速,b,=40(,海里,/,小时,),距离,c,=15(,海里,),求,:,缉私艇的位置,x,(,t,),y,(,t,),缉私艇的航线,y,(,x,),jisi.m,seajisi.m,5实例1海上缉私(续)模型的数值解0yxc(x(t),6,%Creat the function for jisi,%Let x(1)=x,x(2)=y,function dx=jisi(t,x,a,b,c),s=sqrt(c-x(1)2+(a*t-x(2)2);,%dx=b*(c-x(1)/s;b*(a*t-x(2)/s;,dx=1;1;,dx(1)=b*(c-x(1)/s;,dx(2)=b*(a*t-x(2)/s;,6%Creat the function for jisi,7,clear,clf,shg,%Set the definied time,%ts=0:0.05:0.5;,%ts=0:0.1:1.6;,n=length(ts);,x0=0 0;,a=35;b=40;c=15;,opt1=odeset(RelTol,1e-6,AbsTol,1e-9);,t,x=ode45(jisi,ts,x0,opt1,a,b,c);,%a=35;b=40;c=15;,%opt1=odeset(RelTol,1e-6,AbsTol,1e-9);,%t,x=ode45(jisi,ts,x0,opt1,a,b,c);,%exact solution x1=c,y1=a*t;,%output t,x(t),y(t)and draw x(t),y(t),t,x,y1,plot(t,x),grid,gtext(x(t),FontSize,16),gtext(y(t),FontSize,16),pause,%draw y(x):the position of tatch jisi,plot(x(:,1),x(:,2),r*),grid,xlabel(x,FontSize,16),ylabel(y,FontSize,16),7clear,clf,shg,8,实例,1,海上缉私,(,续,),模型的数值解,a,=20,b,=40,c,=15,走私船的位置,x,1,(,t,)=,c,=15,y,1,(,t,)=,at,=20,t,t,=0.5,时缉私艇追上走私船,缉私艇的航线,y,(,x,),t,x,(,t,),y,(,t,),0,0,0,0.05,1.9984,0.0698,0.10,3.9854,0.2924,0.15,5.9445,0.6906,0.20,7.8515,1.2899,0.25,9.6705,2.1178,0.30,11.3496,3.2005,0.35,12.8170,4.5552,0.40,13.9806,6.1773,0.45,14.7451,8.0273,0.50,15.0046,9.9979,y,1,(,t,),0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,8 实例1海上缉私(续)模型的数值解 a=20,b=,9,实例,1,海上缉私,(,续,),模型的数值解,设,b,，,c,不变，,a,变大为,30,35,接近,40,观察解的变化,:,a,=35,b,=40,c,=15,t,=?,缉私艇追上走私船,t,x,(,t,),y,(,t,),y,1,(,t,),0,0,0,0,0.1,3.9561,0.5058,3.5,0.2,7.5928,2.1308,7.0,0.3,10.5240,4.8283,10.5,0.4,12.5384,8.2755,14.0,0.5,13.7551,12.0830,17.5,1.2,14.9986,40.0164,42.0,1.3,14.9996,44.0165,45.5,1.4,15.0117,48.0183,49.0,1.5,15.0023,52.0146,52.5,1.6,14.9866,55.9486,56.0,累积误差较大,提高精度,!,9实例1 模型的数值解 设b，c不变，a变大为30,35,10,实例,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,),判断“追上”的有效方法,?,t,x,(,t,),y,(,t,),y,1,(,t,),0,0,0,0,0.1,3.956104,0.505813,3.5,0.2,7.592822,2.130678,7.0,0.3,10.521921,4.829308,10.5,0.4,12.539454,8.269840,14.0,0.5,13.753974,12.075344,17.5,1.2,14.999616,40.000005,42.0,1.3,14.999963,44.000005,45.5,1.4,14.999993,48.000005,49.0,1.5,14.999998,52.000005,52.5,1.6,15.000020,55.999931,56.0,10实例1模型的数值解 a=35,b=40,c=1,11,实例,1,海上缉私,(,续,),模型的解析解,11实例1 海上缉私(续)模型的解析解,12,实例,1,海上缉私,(,续,),模型的解析解,缉私艇的航线,y,(,x,),的解析解,x=c,时,缉私艇追上走私船的,y,坐标,缉私艇追上走私船的时间,:,a,=20,b,=40,c,=15,t,1,=0.5,a,=35,b,=40,c,=15,t,1,=1.6,12实例1海上缉私(续)模型的解析解 缉私艇的航线y(,13,实 例,2,弱 肉 强 食,问题,自然界中同一环境下两个种群之间的生存方式,相互竞争,相互依存,弱肉强食,弱肉强食,种群甲靠丰富的自然资源生存,食饵,(Prey),种群乙靠捕食种群甲为生,捕食者,(Predator),两个种群的数量如何演变,?,历史背景,一次世界大战期间地中海渔业的捕捞量下降,(,食用鱼和鲨鱼同时捕捞,),，但是其中,鲨鱼的比例却增加，为什么？,13实 例 2弱 肉 强 食 问题自然界中同一环境,14,实 例,2,弱 肉 强 食,模型,食饵,(,甲,),的密度,x,(,t,),捕食者,(,乙,),的密度,y,(,t,),甲独立生存的增长率,r,乙使甲的增长率减小,减小量与,y,成正比,乙独立生存的死亡率,d,甲使乙的死亡率减小,减小量与,x,成正比,Volterra,模型,x,(,t,),y,(,t,),无解析解,14实 例 2弱 肉 强 食 模型食饵(甲)的密,15,实例,2,弱肉强食,模型的数值解,猜测,x,(,t,),y,(,t,),是周期函数,;,y,(,x,),是封闭曲线,数值积分计算一个周期的平均值：,shier.m,shier1.m,15实例2弱肉强食 模型的数值解 猜测 x(t),y(t,16,function y=shier(t,x),r=1;d=0.5;a=0.1;b=0.02;,%y=diag(r-a*x(2),-d+b*x(1)*x;,y=(r-a*x(2)*x(1),(-d+b*x(1)*x(2);,16function y=shier(t,x),17,ts=0:0.1:15;,x0=25,2;,t,x=ode45(shier,ts,x0);,t,x,plot(t,x),grid,gtext(x(t),FontSize,16),gtext(y(t),FontSize,16),pause,plot(x(:,1),x(:,2),grid,xlabel(x,FontSize,16),ylabel(y,FontSize,16),%averge value of x:T is about 10.7,%trapz(t(1:108),x(1:108,:)/10.7,17ts=0:0.1:15;,18,实例,2,弱肉强食,模型的解析解,c,由初始条件确定,相轨线是封闭曲线,(,c,在一定范围内,),求,x,(,t,),y,(,t,),一周期的平均值,:,可以证明,x,(,t,),y,(,t,),是周期函数,(,周期记作,T,),18实例2弱肉强食 模型的解析解 c由初始条件确定 相,19,实例,2,弱肉强食,模型的解析解,x,(,t,),y,(,t,),一周期的平均值,:,r,食饵增长率,a,捕食者对食饵的捕获能力,d,捕食者死亡率,b,食饵对捕食者的喂养能力,结果解释,与计算结果同,既相互制约,又相互依存,19实例2弱肉强食 模型的解析解 x(t),y(t)一周,20,T,2,T,3,T,4,T,1,P,T,1,T,2,T,3,T,4,x,(,t,),的“相位”领先,y,(,t,),进一步分析,初值,相轨线的方向,20T2T3T4T1PT1 T2 T3,21,一次大战期间地中海渔业的捕捞量下降，但是其中,鲨鱼的比例却在增加，为什么？,r,r-,1,d,d+,1,捕捞,战时捕捞,r,r-,2,d,d+,2,2,1,x,y,食饵,(,鱼,),减少，,捕食者,(,鲨鱼,),增加,自然环境,还表明：对,害虫,(,食饵,),益虫,(,捕食者,),系统，使用灭两种,虫的,杀虫剂,会使害虫增加，益虫减少。,21一次大战期间地中海渔业的捕捞量下降，但是其中鲨鱼的比例却,22,食饵,-,捕食者模型,(Volterra),的缺点与改进,Volterra,模型,改写,多数,食饵,捕食者系统观察不到周期震荡,而是趋向某个平衡状态,即存在稳定平衡点,加,Logistic,项,有,稳定平衡点,22食饵-捕食者模型(Volterra)的缺点与改进Volt,23,相轨线是封闭曲线，结构不稳定,一旦离开某一条闭轨线，就进入另一条闭轨线，不恢复原状。,自然界存在的周期性平衡生态系统是结构稳定的，即偏离周期轨道后，内部制约使系统恢