1、实验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公式
2、,实验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,option
3、s) %用5级4阶Runge-Kutta公式 命令的输入f是待解方程写成的函数M文件: function dx=f(t,x) Dx=f1;f2;fn;,实验4-常微分方程数值解,2. Runge-Kutta方法的MatLab实现 举例：仿真模拟著名的Lorenz系统混沌图其中，先建立一个函数M文件function xdot=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实现
4、 画出Lorenz系统图 clear all;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);grid on;%画出三维相图 xlabel(z);ylabel(x);zlabel(y);figure(2) subplot(311);plot(t,x(a:b,1) ;%画三分量演化图 subplo
5、t(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）建立模型
6、,走私船初始位置在点(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千米/秒时， 缉私艇追赶路线的图形。,clear all;clf; a=0.4; v=0.6 0.8 1.0;%取不同的速度 r=0.4./v; t=20*r./(a*(1-r.2)%追上的时间for i=1:3 y=20:-0.01:0; x(:,i)=-0.5
7、*(-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(0 30 0 20); hold on end,追赶时间分别为： T=60.0000,33.3333, 23.8095（秒）,2),当,时，,缉私艇不可能追赶上走私船。,3）,，,，,当,时，,，,缉私艇不可能追赶上走私船。,（b）用MATLAB软件求解析解,MATLAB软件5.3以上版本提供的解常微分方程解析解的指令是 Dsolve, 完整的调用格式是: dsolve(eqn1,eqn2, .) 其中e
8、qn1,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软件防真,当建立动态系统的微分方程模型很困难时, 我们可
9、以用计算机仿真法对系统进行分析研究. 所谓计算机仿真就是利用计算机对实际动态系 统的结构和行为进行编程、模拟和计算, 以此 来预测系统的行为效果.,追赶方向可用方向余弦表示为：%两点形成的向量的方向余弦,时间步长为,，,则在时刻,时：,仿真算法：,第一步：设置时间步长, 速度a, v及初始距离d，,第二步：,计算动点缉私艇D在时刻,时的坐标,，,计算走私船R在时刻,时的坐标,，,第三步：计算缉私艇与走私船这两个动点之间的距离：,根据事先给定的距离，判断缉私艇是否已经追上了走私船，从而判断 退出循环还是让时间产生一个步长，返回到第二步继续进入下一次循环；,第四步：当从上述循环退出后，由点列,和,
10、可分别绘,制成两条曲线即为缉私艇和走私船走过的轨迹曲线。,缉私艇初始位置,，,走私船初始位置,追击问题的数值模拟(P-66) clear;clf; d=120;v=90;a=80;s0=8;%给出初始条件 T=10;dt=0.001;%选取时间区间T(可以偏大一点),时间微元dt t=0:dt:T;%离散时间表t n=length(t); %离散时间表t长度 x(1)=0;y(1)=d;s(1)=s0;%初始位置、初始距离 for i=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
11、*(-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)时刻的距离if d0.1i*dtbreakend%判断是否已追上，并显示追上时的时间s(i)=s0+a*t(i); end plot(x,y) %comet(x,y);,（e） 结果分析,用求解析解的方法算得的解是最为精确的; 用数值方法计算的结果依赖于迭代终值的设定, 减小迭代终值可以提高计算精度;用计算机仿 真法计算的结果依赖于时间迭代步长的选取和 程序终止条件的设定,修改终止条件的设定和 减小时间迭代步长可以提高计算精度,减小误差
12、.,实验4-常微分方程数值解,4. 刚性现象与刚性方程 刚性现象振动系统或包含电容、电感、电阻的电路系统的数学模型一般为：给定一组参数k=2000.5, r=1000, a=1, b=-1999.5,f(t)=1. 则(*1)的解为 稳态解快瞬态解 慢瞬态解 对快瞬态解: 时间常数t1=1/2000=0.0005,计算到t=10t1=0.005时,该项已衰减到 ; 对慢瞬态解: 时间常数t2=2,计算到t=10t2=20时它才衰减到 .用数值方法求解时, 精度要达到 , 至少要算到t=20, 需要14286步, 这样大的计算量是由快瞬态解和慢瞬态解的衰减速度相差悬殊造成的, 这现象称为刚性现象
13、, 相应的微分方程称为刚性方程.,实验4-常微分方程数值解,4. 刚性现象与刚性方程 刚性方程振动、电路及化学反应中出现刚性现象的方程可表示为：（*2） 其中x, f是n维向量, A是nxn矩阵. 当A的特征根 的实部 均为负数时, 方程通解中对应于 的值大的项为快瞬态解,值小的项为慢瞬态解, 称 为刚性比.s10的方程便可认为是刚性方程, 实际问题中可出现s达 的情况. 刚性是问题本身的性质, 与解法无关. 但正是由于这种性质, 用数值方法求解时需要计算到最慢瞬态解衰减成可忽略的小量为止,使得积分区间很长,而为保证计算的稳定性, 当最快瞬态解的 很大时, 又必须使步长充分小, 这就出现了在大
14、区间上用小步长计算的困难情况.,实验4-常微分方程数值解,4. 刚性现象与刚性方程 MatLab求解Matlab中求解常微分方程的命令ode23, ode45. 由于其步长是按稳定性要求和指定的精度加以调整的, 所以用它们解刚性微分方程时步长会自动变小, 对于大的区间会导致计算时间很长.Matlab中有专门求解刚性方程的命令ode23s, ode15s, 用法与 ode23, ode45相同. 例 解方程特征根 , 刚性比 . 解析解为,实验4-常微分方程数值解,4. 刚性现象与刚性方程 MatLab求解 function dx=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;x2 x0= 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薛要很长很长时间.,谢谢同学们！,
