常微分方程数值解与matlab.ppt

实验4-常微分方程数值解,1.求解常微分方程数值方法介绍（1）一阶微分方程求方程(1)的数值解,就是计算(精确)解在一系列离散点的近似值.通常取相等的步长h,于是xn=x0+nh(n=1,2,).(a)欧拉方法基本思想是在小区间xn,xn+1上用差商代替方程(1)左端的导数而方程右端函数f(x,y(x)中的x取xn,xn+1上得某一点,公式为（2）,实验4-常微分方程数值解,(b)Runge-Kutta方法基本思想是用小区间xn,xn+1上的若干个点的导数的线性组合代替方程(2)右端的,一般形式为(3)满足并使(3)的局部截断误差-L级p阶Runge-Kutta公式,实验4-常微分方程数值解,(2)常微分方程组和高阶方程的数值方法欧拉方法和Runge-Kutta方法可直接推广到求常微分方程组,如对欧拉公式为Runge-Kutta公式有类似的形式.对高阶方程(5)需先降阶化为一阶常微分方程组,降阶方法不唯一.简单、常用的方法是令y1=y,将(5)化为,实验4-常微分方程数值解,2.Runge-Kutta方法的MatLab实现对微分方程(组)的初值问题Runge-Kutta方法用MatLab命令实现:t,x=ode23(f,ts,x0,options)%用3级2阶Runge-Kutta公式t,x=ode45(f,ts,x0,options)%用5级4阶Runge-Kutta公式命令的输入f是待解方程写成的函数M文件:functiondx=f(t,x)Dx=f1;f2;fn;,实验4-常微分方程数值解,2.Runge-Kutta方法的MatLab实现举例：仿真模拟著名的Lorenz系统混沌图其中，先建立一个函数M文件functionxdot=lorenz(t,x)sigma=10;r=28;row=8/3;xdot=-sigma*x(1)+sigma*x(2);(r-x(3)*x(1)-x(2);x(1)*x(2)-row*x(3);,实验4-常微分方程数值解,2.Runge-Kutta方法的MatLab实现画出Lorenz系统图clearall;clf;options=odeset(RelTol,1e-5,AbsTol,1e-5);tspan=0,100;x0=1,2,3;t,x=ode45(lorenz,tspan,x0,options);l=length(x(:,1);a=1;b=l;figure(1)plot3(x(a:b,3),x(a:b,1),x(a:b,2),b);gridon;%画出三维相图xlabel(z);ylabel(x);zlabel(y);figure(2)subplot(311);plot(t,x(a:b,1);%画三分量演化图subplot(312);plot(t,x(a:b,2)subplot(313);plot(t,x(a:b,3),实验4-常微分方程数值解,2.Runge-Kutta方法的MatLab实现作业报告：著名的Duffing系统（描述弹簧系统性质）其中类似的，分别画出F=1，2，3，4，6等时的相图翻阅一些参考书，你能得到一些什么结论？,实验4-常微分方程数值解,3.实例问题缉私艇追击走私船海上边防缉私艇发现距d公里处有一走私船正以匀速a沿直线行驶,缉私艇立即以最大匀速度v追赶,在雷达的引导下,缉私艇的方向始终指向走私船.问缉私艇何时追赶上走私船?并求出缉私艇追赶的路线.,S,（1）建立模型,走私船初始位置在点(S0,0),行驶方向为x轴正方向,缉私艇的初始位置在点(0,M0),在时刻t：走私船的位置到达点:(S0+at,0)缉私艇到达点M(x,y),S,（2）模型求解（a）求解析解,令：,令：,（2）模型求解,(a)求解析解,当y=0时:,走私船a=0.4千米/秒，分别取v=0.6,0.8,1.0千米/秒时，缉私艇追赶路线的图形。,clearall;clf;a=0.4;v=0.60.81.0;%取不同的速度r=0.4./v;t=20*r./(a*(1-r.2)%追上的时间fori=1:3y=20:-0.01:0;x(:,i)=-0.5*(-40*r(i)+20(-r(i)*(r(i)-1)*y.(1+r(i)+20r(i)*(r(i)+1)*y.(1-r(i)/(1-r(i)2);plot(x(:,i),y);axis(030020);holdonend,追赶时间分别为：T=60.0000,33.3333,23.8095（秒）,2),当,时，,缉私艇不可能追赶上走私船。,3）,，,，,当,时，,，,缉私艇不可能追赶上走私船。,（b）用MATLAB软件求解析解,MATLAB软件5.3以上版本提供的解常微分方程解析解的指令是Dsolve,完整的调用格式是:dsolve(eqn1,eqn2,.)其中eqn1,eqn2,.是输入宗量,包括三部分:微分方程、初始条件、指定变量,若不指定变量,则默认小写字母t为独立变量.书P-69,微分方程的书写格式规定:当y是因变量时,用“Dny”表示y的n阶导数.,例求微分方程,的通解。,dsolve(Dy=x+x*y,x),Ans=-1+exp(1/2*x2)*C1,dsolve(Dx=1/2*(y/20)r-(20/y)r),x(20)=0,y),Ans=1/2*20(-r)*y(r+1)/(r+1)+1/2*20r/(r-1)*y*(1/y)r-20*r/(r2-1),（c）用MATLAB软件防真,当建立动态系统的微分方程模型很困难时,我们可以用计算机仿真法对系统进行分析研究.所谓计算机仿真就是利用计算机对实际动态系统的结构和行为进行编程、模拟和计算,以此来预测系统的行为效果.,追赶方向可用方向余弦表示为：%两点形成的向量的方向余弦,时间步长为,，,则在时刻,时：,仿真算法：,第一步：设置时间步长,速度a,v及初始距离d，,第二步：,计算动点缉私艇D在时刻,时的坐标,，,计算走私船R在时刻,时的坐标,，,第三步：计算缉私艇与走私船这两个动点之间的距离：,根据事先给定的距离，判断缉私艇是否已经追上了走私船，从而判断退出循环还是让时间产生一个步长，返回到第二步继续进入下一次循环；,第四步：当从上述循环退出后，由点列,和,可分别绘,制成两条曲线即为缉私艇和走私船走过的轨迹曲线。,缉私艇初始位置,，,走私船初始位置,追击问题的数值模拟(P-66)clear;clf;d=120;v=90;a=80;s0=8;%给出初始条件T=10;dt=0.001;%选取时间区间T(可以偏大一点),时间微元dtt=0:dt:T;%离散时间表tn=length(t);%离散时间表t长度x(1)=0;y(1)=d;s(1)=s0;%初始位置、初始距离fori=1:nx(i+1)=x(i)+v*dt*(s0+a*t(i)-x(i)/sqrt(s0+a*t(i)-x(i)2+y(i)2);y(i+1)=y(i)+v*dt*(-y(i)/sqrt(s0+a*t(i)-x(i)2+y(i)2);%递推算式、d=sqrt(s0+a*t(i)-x(i)2+y(i)2);%t(i)时刻的距离ifd10的方程便可认为是刚性方程,实际问题中可出现s达的情况.刚性是问题本身的性质,与解法无关.但正是由于这种性质,用数值方法求解时需要计算到最慢瞬态解衰减成可忽略的小量为止,使得积分区间很长,而为保证计算的稳定性,当最快瞬态解的很大时,又必须使步长充分小,这就出现了在大区间上用小步长计算的困难情况.,实验4-常微分方程数值解,4.刚性现象与刚性方程MatLab求解Matlab中求解常微分方程的命令ode23,ode45.由于其步长是按稳定性要求和指定的精度加以调整的,所以用它们解刚性微分方程时步长会自动变小,对于大的区间会导致计算时间很长.Matlab中有专门求解刚性方程的命令ode23s,ode15s,用法与ode23,ode45相同.例解方程特征根,刚性比.解析解为,实验4-常微分方程数值解,4.刚性现象与刚性方程MatLab求解functiondx=stiff1(t,x)dx=x(1)+2*x(2);-(106+1)*x(1)-(106+2)*x(2);t=0:0.1:1;%t=0:0.1:10;x1=(106/4+1)*exp(-t)-exp(-106*t);x2=-(106/4+1)*exp(-t)+(106+1)/2*exp(-106*t);A=t;x1;x2x0=106/4,106/4-1/2;t,x=ode23s(stiff1,t,x0);%很快出结果B=t,x%t,y=ode23(stiff1,t,x0);%几十秒出结果%C=t,y%要计算0,10才能保证精度,ode23薛要很长很长时间.,谢谢同学们！,