常微分方程的数值解法

1、第八章 常微分方程的数值解法一内容要点 考虑一阶常微分方程初值问题: 微分方程的数值解：设微分方程的解y(x)的存在区间是a,b，在a,b内取一系列节点a= x0 x1 xn =b,其中hk=xk+1-xk ；(一般采用等距节点，h=(b-a)/n称为步长)。在每个节点xk求解函数y(x)的近似值:yky(xk)，这样y0 , y1 ,.,yn称为微分方程的数值解。用数值方法，求得f(xk)的近似值yk，再用插值或拟合方法就求得y(x)的近似函数。(一)常微分方程处置问题解得存在唯一性定理对于常微分方程初值问题:如果：(1) 在的矩形内是一个二元连续函数。(2) 对于y满足利普希茨条件，即则在

2、上方程的解存在且唯一，这里C=min(A-x0),x0+B/L),L是利普希茨常数。定义：任何一个一步方法可以写为，其中称为算法的增量函数。收敛性定理：若一步方法满足：(1)是p解的.(2) 增量函数对于y满足利普希茨条件.(3) 初始值y0是精确的。则kh=x-x0,也就是有(一)、主要算法1局部截断误差局部截断误差：当y(xk)是精确解时，由y(xk)按照数值方法计算出来的的误差y(xk+1)- 称为局部截断误差。注意：yk+1和的区别。因而局部截断误差与误差ek+1=y(xk+1) -yk+1不同。如果局部截断误差是O(hp+1),我们就说该数值方法具有p阶精度。1. 显式欧拉格式： k

3、=1,2,n-1.显式欧拉格式的特点：(1)、单步方法；(2)、显式格式；(3)、局部截断误差因而是一阶精度 。隐式欧拉格式2. 两步欧拉格式： k=1,2,n-1.两步欧拉格式的局部截断误差，因而是二阶精度3. 梯形方法: ；3.改进的欧拉方法：预测值: 校正值: 。或改写为: 4、梯形方法与改进欧拉方法的截断误差是O(h3),具有二阶精度。5、龙格-库塔法的思想 1). 二阶龙格-库塔法计算公式： 当:时，得一簇龙格-库塔公式，其局部截断误差均为O(h3)都是二阶精度。特别取,就是改进欧拉公式。取,得二阶龙格库塔法为：，称为二阶中点格式。2)、经典龙格库塔格式(也称为标准龙格库塔方法)：四

4、阶龙格库塔方法的截断误差为，具有四阶精度。一般一阶常微分方程初值问题用四阶龙格库塔方法计算，其精度均满足实际问题精度要求。3).变步长龙格库塔方法： 从节点xk出发，以步长h据四阶龙格库塔方法求出一个近似值,然后以步长h/2求出一个近似值，得误差事后估计式：根据来选取步长h。4). RKF格式：变步长龙格库塔方法，因频繁加倍或折半步长会浪费计算量。Felhberg改进了传统龙格库塔方法，得到RKF格式，较好解决了步长的确定，而且提高了精度与稳定性，为Matlab等许多数值计算软件采用。4/5阶RKF格式：由4阶龙格库塔方法与5阶龙格库塔方法结合而成。Felhberg得到的最佳步长hs，其中h为

5、当前步长，.e为精度要求，若s1.5步长加倍。6.亚当斯方法(Adams)1).显式Adams方法：记：；二阶显式Adams方法：;三阶显式Adams方法：;四阶显式Adams方法：.2). 隐式Adams方法二阶隐式Adams方法：;三阶隐式Adams方法：;四阶隐式Adams方法：3）.Adams预报-校正系统：先用显式格式作为预测值，再用隐式格式来校正。预测值： 校正值： 4).改进的Adams预报-校正系统：预测：改进：校正值： 改进：7.收敛与稳定性对于固定的，如果数值解yk 当（同时）时趋向于准确解y(xk)，则称该数值方法方法是收敛的。如果一种数值方法，在节点值yk上大小为d的扰

6、动，于以后各节点值ym(mk)上产生的偏差均不超过d，则称数值方法是稳定的. 一般，隐式格式较显式格式有较好的稳定性。8.刚性方程组：考虑n阶常微分方程组：，（*）若矩阵A的特征值l1,l2,ln的实部Re(li)e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1;if kK,error(迭代发散);endend改进Euler格式functionx,y=naeulerg(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 K1=feval(dyfun,x(n),y(n); y(n+1)=

7、y(n)+h*K1; K2=feval(dyfun,x(n+1),y(n+1); y(n+1)=y(n)+h*(K1+K2)/2;endx=x;y=y;4阶经典Runge-Kutta格式functionx,y=nak4(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 K1=feval(dyfun,x(n),y(n); K2=feval(dyfun,x(n)+h/2,y(n)+h/2*K1); K3=feval(dyfun,x(n)+h/2,y(n)+h/2*K2); K4=feval(dyfun,x(n+1)

8、,y(n)+h*K3); y(n+1)=y(n)+h*(K1+2*K2+2*K3+K4)/6;endx=x;y=y;变步长4阶经典Runge-Kutta格式functionx,y=nark4v(dyfun,xspan,y0,e,h)if nargin5,h=(xspan(2)-xspan(1)/10;endn=1;x(n)=xspan(1);y(n)=y0;y1,y2=comput(dyfun,x(n),y(n),h);while x(n)e while abs(y2-y1)/10e h=h/2; y1,y2=comput(dyfun,x(n),y(n),h); endelse while a

9、bs(y2-y1)/10=e h=2*h; y1,y2=comput(dyfun,x(n),y(n),h); end h=h/2;h=min(h,xspan(2)-x(n); y1,y2=comput(dyfun,x(n),y(n),h);end n=n+1;x(n)=x(n-1)+h;y(n)=y2; y1,y2=comput(dyfun,x(n),y(n),h); endx=x;y=y;function y1,y2=comput(dyfun,x,y,h); y1=rk4(dyfun,x,y,h);y21=rk4(dyfun,x,y,h/2);y2=rk4(dyfun,x+h/2,y21,h

10、/2);function y=rk4(dyfun,x,y,h); K1=feval(dyfun,x,y); K2=feval(dyfun,x+h/2,y+h/2*K1); K3=feval(dyfun,x+h/2,y+h/2*K2); K4=feval(dyfun,x+h,y+h*K3); y=y+h*(K1+2*K2+2*K3+K4)/6;五.试验例题例1 不同方法的精度比较 用多种方法解下述初值问题，并与其准确解进行比较。 解：取步长h=0.1, xk=kh (k=0,1,6)。各方法计算结果及绝对误差见表(1)、(2)、(3)。表（1）xn欧拉公式改进的欧拉公式四阶标准龙格库塔公式yn误

11、差yn误差yn误差0.01.000000010100.11.0000001.00500001.004837500.21.0100001.0190251.018730900.31.0290001.0412181.040818420.41.0561001.0708021.070320290.51.0904901.1070761.106530930.61.1314411.1494041.14881193表（2） 四阶阿当姆斯公式nxn显示公式隐式公式yn误差yn误差30.3取自准确解1.0408180140.41.070322921.0703196650.51.106535481.106530416

12、0.61.148814811.14881101备 注y1, y2, y3 取自准确解y1, y2 取自准确解 本题 关于y为线性，代入隐式公式，可解出yn+1，因此可直接求隐式公式之解。表（3） 四阶阿当姆斯预测校正公式nxn显示公式隐式公式Yn误差yn误差40.41.070319921.0703201450.51.106530271.1065307760.61.148811031.14881175备 注y1, y2, y3由四阶标准龙格库塔公式提供 对比以上各表数据知，在相同步长下求解同一问题时，方法阶数越高，解的精度也越高，一阶的欧拉公式精度最低，四阶标准龙格库塔公式的精度大大高于改进的欧

13、拉公式。四阶阿当姆斯方法、显式的精度略低于隐式。同是阿当姆斯预测校正系统，带误差修正的精度又高于不带误差修正的。 从计算量看，四阶龙格库塔法计算量最大，每前进一步要计算4次函数值f，而带误差修正的阿当姆斯预测校正法与它精度差不多的，每前进一步只计算两次函数值f，所以后者是可取的。但它是四步法，必须用其它方法提供出始值，程序略复杂些。 例2 (步长的计算结果的影响)用欧拉公式求下述初值问题在x = 1处的近似解，并与准确解y (1)=1比较。 解 分别取步长h=10-1、h=10-2、h=10-3、h=10-4计算，结果见下表。由表中数据可见，当 时结果完全失真，而取 比 计算量增加了十倍，但解

14、的精度却基本一样，可见取 太浪费计算量了。步长hy (1) 的计算值 上述结果差异很大的原因在于欧拉公式的绝对稳定区间为（2, 0）步长h应满足 ，对本题， ，故应取h满足即 。可见取 时欧拉公式是数值不稳定的，导致结果失真，而取 和 都满足稳定性要求，可用于求解。此例说明，求解微分方程数值解一定选取注意步长，过大则导致解的失真，过小又会使计算量大增。究竟取多大步长才合适，不仅取决于所采用的数值方法，还决定于待解微分方程本身的特性。例3：(Matlab不同命令之差异)取h=0.2分别用不同的Matlab命令解方程 (0t t,y=ode45(odefun,0,4,1);n=length(t);

15、e=sqrt(sum(sqrt(1+2*t)-y).2/n);n,eans =45.00000000000000 0.00020257808813 t,y=ode23(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n);n,eans =13.00000000000000 0.19048360113013 t,y=ode113(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n);n,eans =18.00000000000000 0.00973277413962 t,y=o

16、de23t(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n);n,eans =18.00000000000000 0.03919668151270 t,y=ode23s(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n);n,eans =81.00000000000000 2.54374270001824 t,y=ode23tb(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n);n,eans =15.000

17、00000000000 0.24308870757312 t,y=ode15s(odefun,0,4,1);n=length(t);e=sqrt(sum(sqrt(1+2*t)-y).2/n);n,eans =22.00000000000000 0.45509043060378可见ode45精度较高，计算量较大，ode23计算量小，但误差大，ode113适中，用刚性方程组揭发解非刚性问题不适合，特别ode23s计算量大且但误差大。例4：解刚性方程组： 其准确解为:就h=0.01、h=0.02使用改进欧拉格式及方程组，并在一张图上绘出其图形。functionx,y=naeuler2s(dyfun

18、,xspan,y0,h)x=xspan(1):h:xspan(2);y=zeros(length(y0),length(x);y(:,1)=y0(:);for n=1:length(x)-1 K1=feval(dyfun,x(n),y(:,n); y(:,n+1)=y(:,n)+h*K1; K2=feval(dyfun,x(n+1),y(:,n+1); y(:,n+1)=y(:,n)+h*(K1+K2)/2;endx=x;y=y; %用欧拉格式解常微分方程组的m文件function f=dyfun(x,y)f(1)=-0.01*y(1)-99.99*y(2);f(2)=-100*y(2);f=

19、f(:); %函数m文件clear; x,y=naeuler2s(dyfun,0 500,2,1,0.01); plot(x,y),axis(-50 500 -0.5 2) gtext(y(x) h=0.01), gtext(z(x) h=0.01) %步长为0.01时解方程并画图 clear; x,y=naeuler2s(dyfun,0 500,2,1,0.02); hold on plot(x,y),gtext(y(x) h=0.02), gtext(z(x) h=0.02) %步长为0.02时解方程并画图。clear; x=0:0.01:500; y=exp(-0.01*x)+exp(-

20、100*x);z=exp(-100*x); hold on plot(x,y,r,x,z,r),gtext(y(x) 真根), gtext(z(x) 真根) %画出真根图形运行结果真根与h=0,01时结果重合，h=0.02时结果是错误的。故刚性方程组应采用稳定性较好方法。例6：数值实验（摘自蔡）观察欧拉显示方法的收敛性。内容：取h=0.1, 0.01,用欧拉显式方法求解：用欧拉法解方程 (0x1). 计算到，并与精确解比较。MATLAB程序t=1; y=1;h=0.1 for k=1:10; y=y+h*t*y(1/3) t=t+hend运行结果hyt0.10.010.0010.0001y=(

21、t2+2)/3)3/2(准确解)1111111.11.11.11751.10681.10681.106816606308381.21.21361.23931.22771.22791.227879593566201.31.34161.37621.36391.36411.364135990288361.41.48491.52941.51621.51651.516564538686041.51.64471.69981.68571.68611.686170601183731.61.82171.88831.87341.87391.873981856902571.72.0172.09622.08042.0

22、8102.081044689573001.82.23192.32432.30762.30832.308421169608421.92.46712.57392.55632.55712.557186539930162.02.72392.81772.82742.82832.82842712474619欧拉显示方法是收敛的。例7：实验目的：比较不同算法用于“非光滑”解时的结果。内容：用欧拉显式方法和2阶龙格库塔算法求解： 计算并在一张图上画出解曲线。k=6时的程序如下：cleart=0;y=1;h=pi/60; for k=0.1:0.1:2; y=y+h*abs(sin(6*t) t=t+h plo

23、t(t,y,o), hold on,endt0=0;tf=1/3*pi;x0=1;k=pi/60t,x=ode23(ex545f,t0:k:tf,1) plot(t,x), hold on,for i=0:pi/60:pi/3 if i T,Y=ode45(rigid,0 12,0 1 1); plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+)3、结果如图.【例11】采用最简洁格式的ODE文件和解算指令，研究围绕地球旋转的卫星轨道。 （1）问题的形成轨道上运动的卫星，在Newton第二定律，和万有引力定律作用下，有。即，而。这里是引力常数，是地球的质量。又假定卫星以初

24、速度在处进入轨道。（2）构成一阶微分方程组令，则(5.14.2.1-5)初始条件为(5.14.2.1-6)（1）根据式(5.14.2.1-5)编写最简洁格式的ODE文件DYDt2.mfunction Yd=DYDt2(t,Y,flag,G,ME) % flag按ODE文件格式规定，必须是第三输入宗量。对它的赋值由ode45指令自动产生。 第4、5宗量是被传递的参数switch flagcase %按规定：这里必须是空串。在此为真时，完成以下导数计算。X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.2);Yd=V;-G*ME*X/r3;otherwise error(Unknown

25、 flag flag .);end（2）对微分方程进行解算G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;tspan=t0,tf;%指定解算微分方程时间区间Y0=x0;0;0;vy0;%按式(5.14.2.1-6)给定初值向量t,YY=ode45(DYDt2,tspan,Y0,G,ME);%第4宗量取缺省设置,第5、6宗量是被传递参数X=YY(:,1);%输出Y第一列是位移数据x(t)Y=YY(:,2);%输出Y第二列是位移数据y(t)plot(X,Y,b,Linewidth,2); hold onaxis(image)

26、%保证x、y轴等长刻度，且坐标框恰包容图形XE,YE,ZE = sphere(10);%产生单位球面数据RE=0.64e7;%地球半径XE=RE*XE;YE=RE*YE;ZE=0*ZE;%坐标纸上的地球平面数据mesh(XE,YE,ZE),hold off%绘地球示意图 【*例12】带事件设置的ODE文件及主程序编写演示。本例将以较高精度计算卫星经过近地点和远地点的时间，并在图上标志。（1）ODE文件的编写DYDt3.mfunction varargout=DYDt3(t,Y,flag,G,ME,tspan,Y0)% DYDt3.m 供主程序调用的ODE函数文件%本文件自带三个子函数：f,fi

27、,fev。%t, Y分别是自变量和一阶函数向量，是最基本的输入宗量。%flag第三输入宗量，它专供解算指令（如ode45）作调用通知。%在运行中，解算指令会根据需要向flag发不同的字符串。% varargout 是变输出宗量。它由变维的元胞数组构成。每个元胞中可以存放指令%所产生的任意形式的数据。switch flagcase % 必须用空串符。情况为真时，计算导数 dY/dt = f(t,Y)。 varargout1 = f(t,Y,G,ME);%输出为一个元胞，容纳f子函数的一个输出Ydcase init % 必须用init，情况为真时，传送计算区间、初值、设置参数。 varargout

28、1:3 = fi(tspan,Y0);%输出为三个元胞，容纳fi子函数的三个输出量case events% 必须用events，情况为真时，设置事件性质。 varargout1:3 = fev(t,Y,Y0);%输出三个元胞，容纳fev子函数的三个输出量otherwise error(Unknown flag flag .);end% -function Yd = f(t,Y,G,ME)%计算导数子函数，被父函数DYDt3调用。X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.2);Yd=V; -G*ME*X/r3;% -function ts,y0,options = fi(ts

29、pan,Y0)%设置时间区间、初值、算法参数子函数，被父函数DYDt3调用。ts=tspan;y0 = Y0;options = odeset(Events,on,Reltol,1e-5,Abstol,1e-4);%开动ode45的事件判断功能，设置相对误差和绝对误差。% -function value,isterminal,direction = fev(t,Y,Y0)%事件判断子函数，被父函数DYDt3调用。dDSQdt = 2 * (Y(1:2)-Y0(1:2) * Y(3:4);%dDSQdt是离初始点的二乘距离u的时间导数du/dt，而u=(x-x0)2+(y-y0)2 。value

30、 = dDSQdt; dDSQdt;%定义两个穿越0的事件direction = 1; -1;%第一事件：以渐增方式穿越0。第二事件：以渐减方式穿越0isterminal = 1; 0;%第一事件发生后，终止计算；而第二事件发生后，继续计算。（2）运行以下主程序G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;tspan=t0,tf;Y0=x0;0;0;vy0;t,YY,Te,Ye,Ie=ode45(DYDt3,G,ME,tspan,Y0);%X=YY(:,1);Y=YY(:,2);plot(X,Y,b,Linewidth

31、,2);hold ontext(0,6e7,轨道,Color,b)%产生蓝色文字注释axis(image);%保证x、y轴等长刻度，且坐标框恰包容图形在三个事件发生点上画标记plot(Ye(1,1),0.4e7+Ye(1,2),r,MarkerSize,10)plot(Ye(2,1),0.4e7+Ye(2,2),bv,MarkerSize,10)plot(Ye(3,1),-0.4e7+Ye(3,2),b,MarkerSize,10)%把轨道的半周期和全周期标在图上text(0.8*Ye(3,1),-2e7+Ye(3,2),t3= sprintf(%6.0f,Te(3)text(0.8*Ye(2

32、,1),1.5e7+Ye(2,2),t2= sprintf(%6.0f,Te(2) %在x-y坐标上画地球 XE,YE,ZE= sphere(10);RE=0.64e7;XE=RE*XE;YE=RE*YE;ZE=0*ZE;mesh(XE,YE,ZE)text(1e7,1e7,地球,Color,r),hold off%产生红色文字注释 例13 范例：状态转移方程组模型 计算机网络可靠性分析问题随着计算机通信网络系统特别是Internet网络日益广泛的应用，提高系统的可靠性意义重大研究和分析具有实用性的高可靠计算机通信网络系统，是国际上非常活跃的一个研究方向（1） 问题及假设 我们将研究的计算机通

33、信网络系统称为无冗余的防火墙协议系统计算机随时可能发生三个事件一无故障、间歇故障和永久故障因此，计算机一般处于三种工作状态：无故障工作、带故障工作和不工作这三种状态之间的转移过程如图39所示要求建立该系统的状态转移模型，并进行可靠性分析 (2) 分析与模型该问题属于状态转移问题，利用马尔科夫状态转移原理，用,和分别表示系统处于无故障工作、带故障工作和不工作三种状态的概率，则有以下状态转移方程组： 初始条件,参数取值这是一个带参数的微分方程组模型(3) 模型求解 求解这个带参数的微分方程组模型有两种方法，一是用特征根法求解析解，另一个是用数值解法求数值解分别求解如下： 假定模型中的参数取下限值，

34、即 .解析解法 利用常系数线性微分方程组的特征根法MATLAB程序； lp=10(-5);lt=10(-4);gm=0.01; A=-(lp+lt),gm,0;lt/2,-(gm+lp+lt),0;lp+lt/2,lp+lt,0;l,v=eig(A); c=inv(1)*1;0;0输出结果：特征值（0，-0.0102,-0.0001）与之对应的特征向量(0 0 1); （07053 -07089 00035）； （-07053 -00035 07089）根据特征值和特征向量写出其解析解,利用初始条件得到特解 数值解法 MATLAB程序如下首先编写m文件（eqs0m）function xdot=

35、eqs0(t,p,flag,lp,lt,gm)a=-(lp=lt),gm,0;lt/2,-(gm+lp+lt),0;lp+lt/2,lp+lt,0;p=p(1);p(2);p(3)xdot=a*p;在工作空间执行以下程序：ts=0 1000;p0=1;0;0;lp=10(-5);lt=10(-4);gm=0.01;t,p=ode23(eqs0,ts,p0,lp,lt,gm)plot(t,1-p(:,);xlabel(时间t(小时);ylabel(可靠度R(t);title(参数取值lp=0.00001;lt=0.00001;gm=0.01);grid输出结果如图3.10计算表明：当t=1000

36、小时情况下， 显示系统工作的概率为94.18%，从图3.10中还可以看出系统可靠性变化更加细致的情况 , 何时可靠度达到99%，何时达到98%等等。 .敏感性分析 当参数取值在以下范围 内变化时，系统的可靠性情况如何呢？留给读者自己分析例14 解: 令 y1=x，y2=y1 则微分方程变为一阶微分方程组： 1、建立m-文件vdp1000.m如下： function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)2)*y(2)-y(1); 2、取t0=0，tf=3000，输入命令： T,Y=ode15s(vdp1000,

37、0 3000,2 0); plot(T,Y(:,1),-)3、结果如图程序见数学建模高教出版社数学建模与数学实验To Matlab（ff4）导弹追踪问题设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中？解法一（解析法）假设导弹在t时刻的位置为P(x(t), y(t)，乙舰位于由于导弹头始终对准乙舰，故此时直线PQ就是导弹的轨迹曲线弧OP在点P处的切线，即有 即 (1)又根据题意，弧OP的长度为的5倍，即 (2)由(1),(2)消去t整理得模型: 初值条件为: 解即为导弹的运行轨迹: 当时，即当乙舰航行到点处时被导弹击中. 被击中时间为:. 若v0=1, 则在t=0.21处被击中.轨迹图见程序chase1解法二(数值解)令y1=y,y2=y1，将方程（3）化为一阶微分方程组。 1.建立m-文件eq1.m function dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1/5*sqrt(1+y(1)2)/(1-x); 2. 取x0=0