常微分方程数值解法
第八章 常微分方程的数值解法
一．内容要点
考虑一阶常微分方程初值问题:?????==0
0)()
,(y x y y x f dx dy
微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节
点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。
用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。 (一)常微分方程处置问题解得存在唯一性定理
对于常微分方程初值问题:?????==0
0)()
,(y x y y x f dx dy
如果：
(1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。 (2) ),(y x f 对于y 满足利普希茨条件，即
2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程?????==0
0)()
,(y x y y x f dx
dy
的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。
定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。
收敛性定理：若一步方法满足： (1)是p 解的.
(2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件.
(3) 初始值y 0是精确的。则),()()(p h O x y kh y =-kh =x -x 0,也就是有
0x y y lim k x x kh 0h 0
=--=→)(
(一)、主要算法 1．局部截断误差
局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~
+k y 的误差y (x k+1)- 1~
+k y 称为局部截断误差。
注意：y k+1和1~
+k y 的区别。因而局部截断误差与误差e k +1=y (x k +1) -y k +1不同。
如果局部截断误差是O (h p+1),我们就说该数值方法具有p 阶精度。
1. 显式欧拉格式：?
??=+=+001)()
,(y x y y x hf y y k k k k k =1,2,?n-1.
显式欧拉格式的特点：
(1)、单步方法； (2)、显式格式；
(3)、局部截断误差)(2h O 因而是一阶精度 。
隐式欧拉格式?
??=+=+++00111)()
,(y x y y x hf y y k k k k
2. 两步欧拉格式：?
??=+=-+001)()
,(y x y y x hf 2y y k k 1k k k =1,2,?n-1.
两步欧拉格式的局部截断误差)(3h O ，因而是二阶精度 3. 梯形方法:
=+1k y []??
???
=+?+=+++00)(),(),(y x y y x f y x f 2h y y k k 1k 1k k 1
k ； 3.改进的欧拉方法：
预测值: ),(^
k k k 1k y x f h y y ?+=+
校正值: =+1k y ??
?
???+?+++),(),(^k k 1k 1k k y x f y x f 2h y 。
或改写为:???
?
???
+=+=+=++)
(21
),(),(11c p k p k k c k k k p y y y y x f h y y y x f h y y
4、梯形方法与改进欧拉方法的截断误差是O(h 3),具有二阶精度。
5、龙格-库塔法的思想
1). 二阶龙格-库塔法计算公式： ???
??λ+λ-+=++==+)
)((),(),(211
121K K 1h y y hK p y h p x f K y x f K k k k k k k
当:2
1
=λp 时，得一簇龙格-库塔公式，其局部截断误差均为O (h 3)都是二阶精度。特别取1p ==
λ,21
,就是改进欧拉公式。 取1p =λ=,2
1
,得二阶龙格—库塔法为：
?
????+=++==+2
1121)2,2(),(hK y y K h y h x f K y x f K k k k
k k k ，称为二阶中点格式。 2)、经典龙格－库塔格式(也称为标准龙格－库塔方法)： ????
?
?
?
????++++=++=++=++==+)22(6),()2,2()2,2(),(432113423121K K K K h y y hK y h x f K K h y h x f K K h y h x f K y x f K k
k k k k
k k
k k k 四阶龙格－库塔方法的截断误差为()
5h O ，具有四阶精度。一般一阶常微分方程初值问题用四阶龙格－库塔方法计算，其精度均满足实际问题精度要求。
3).变步长龙格－库塔方法： 从节点x k 出发，以步长h 据四阶龙格－库塔方法求出一
个近似值)
(1h k y +,然后以步长h/2求出一个近似值)2/(1h k y +，得误差事后估计式：
)()
(1)2/(1)
2/(1
h k h k h k 1k y y 15
1y y ++++-≈- 根据)
(1)2/(1h k h k y y ++-=?来选取步长h 。
4). RKF 格式：变步长龙格－库塔方法，因频繁加倍或折半步长会浪费计算量。Felhberg 改进了传统龙格－库塔方法，得到RKF 格式，较好解决了步长的确定，而且提高了精度与稳定性，为Matlab 等许多数值计算软件采用。
4/5阶RKF 格式：由4阶龙格－库塔方法与5阶龙格－库塔方法结合而成。
?
??
??+-+++=?
?
?
??-+++=++654311^
5431155250956430285611282566561351651410121972565140821625K K K K K h y y K K K K h y y n n n n
??
?
?????
????
????
?
-+-+-+=-+-++=+-++=++
+=++==)401141041859256535442278,2()
410484551336808216439,()219772962197720021971932,1312()32
9
323,83()
4,4(),(443216432153214213121K K K K K y h x f K K K K K y h x f K K K K y h x f K K K y h x f K K h y h x f K y x f K n n n n n n
n n n n n n
Felhberg 得到的最佳步长hs ，其中h 为当前步长，4
/111^
?????
?
?
?-ε
=++n n y y h s .ε为精度要求，
若s<0.75,步长折半；若s>1.5步长加倍。 6.亚当斯方法(Adams) 1).显式Adams 方法： 记：),(k k k y x f f =；
二阶显式Adams 方法：][211-+-+=k k k k f f 3h
y y ;
三阶显式Adams 方法：][2211--++-+=k k k k k f 5f 16f 231h
y y ;
四阶显式Adams 方法：]5559379[24
1231k k k k k k f f f f h
y y +-+-+=---+.
2). 隐式Adams 方法
二阶隐式Adams 方法：][211++-+=k k k k f f 3h
y y ;
三阶隐式Adams 方法：]8[2111-++-++=k k k k k f f f 51h
y y ;
四阶隐式Adams 方法：)9195(24
1121+--+++-+=k k k k k k f f f f h
y y
3）.Adams 预报-校正系统：
先用显式格式作为预测值，再用隐式格式来校正。
预测值： )9375955(24
3211---+-+-+=k k k k k k f f f f h
y y
校正值：
()),9195(24
11121++--+++-+
=k k k k k k k y x f f f f h
y y 4).改进的Adams 预报-校正系统：
预测：)9375955(243211---+-+-+=k k k k k k f f f f h
y p
改进：)(k k 1k 1k p c 270251
p m -+=++ 校正值：
())519,9(24
21111--++++-++
=k k k k k k k f f f m x f h
y c 改进：)(1k 1k 1k 1k p c 270
19
c y ++++--= 7.收敛与稳定性
对于固定的kh x x k +=0，如果数值解y k 当,0h →（同时+∞→n ）时趋向于准确解y (x k )，则称该数值方法方法是收敛的。
如果一种数值方法，在节点值y k 上大小为δ的扰动，于以后各节点值y m (m >k )上产生的偏差均不超过δ，则称数值方法是稳定的.
一般，隐式格式较显式格式有较好的稳定性。 8.刚性方程组：
考虑n 阶常微分方程组：)(x Ay y φ+='，???（*）
若矩阵A 的特征值λ1,λ2,?λn 的实部Re(λi )<0，i=1,2,?,n.，则称)
Re()
Re(i n
i 1i n i 1min max S λλ=
≤≤≤≤为*式
的刚性比，当刚性比很大时称*式为刚性方程组。 刚性方程组需采用稳定性较好的方法。
二. MATLAB 的有关命令
[t ，x]=solver （’f’,ts,x0,options ）
t:输出的自变量值, x: 输出的函数值；
solver :包括ode45、 ode23 、ode113、ode15、sode23s.
非刚性系统：
ode23：用2阶(及3阶)龙格-库塔算法。 ode45：用4阶(及5阶)龙格-库塔-算法。
ode113(Adams-Bashforth-Moulton PECE)多步方法。 刚性系统：
ode15s(Gear 方法)多步方法
ode23s(二阶modified Rosenbroack formula)单步。 ode23t(trapezoidal rule)solve DAEs. ode23tb(TR-BDF2)低精度算法。 f ：由待解方程写成的m-文件名
ts=[t0，tf]，t0、tf 为自变量的初值和终值 x0：函数的初值
options=odeset(’reltol’,rt,’abstol’,at ), rt ，at ：分别为设定的相对误差和绝对误差. (缺省时设定相对误差10-3, 绝对误差10-6) 注意:
1、在解n 个未知函数的方程组时，x0和x 均为n 维向量，m-文件中的待解方程组应以x 的分量形式写成.
2、使用Matlab 软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.
四.重点算法Matlab 程序
Euler 格式、
function[x,y]=naeuler(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0;
for n=1:length(x)-1
y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end
x=x';y=y';
隐式Euler 格式
function[x,y]=naeulerb(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0;
for n=1:length(x)-1
y(n+1)=iter(dyfun,x(n+1),y(n),h);
end
x=x';y=y';
function y=iter(dyfun,x,y,h)
y0=y;e=1e-4;K=1e+4;
y=y+h*feval(dyfun,x,y);
y1=y+2*e;k=1;
while abs(y-y1)>e
y1=y;
y=y0+h*feval(dyfun,x,y);
k=k+1;if k>K,error('迭代发散');end
end
改进Euler格式
function[x,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)=y(n)+h*K1;
K2=feval(dyfun,x(n+1),y(n+1));
y(n+1)=y(n)+h*(K1+K2)/2;
end
x=x';y=y';
4阶经典Runge-Kutta格式
function[x,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),y(n)+h*K3);
y(n+1)=y(n)+h*(K1+2*K2+2*K3+K4)/6; end
x=x';y=y';
变步长4阶经典Runge-Kutta格式
function[x,y]=nark4v(dyfun,xspan,y0,e,h)
if nargin<5,h=(xspan(2)-xspan(1))/10;end
n=1;x(n)=xspan(1);y(n)=y0;
[y1,y2]=comput(dyfun,x(n),y(n),h);
while x(n) if abs(y2-y1)/10>e while abs(y2-y1)/10>e h=h/2; [y1,y2]=comput(dyfun,x(n),y(n),h); end else while abs(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); end x=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/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 不同方法的精度比较 用多种方法解下述初值问题，并与其准确解x e y x +=-进行比较。 ? ??=≤<++-='1)0(6 .00,1y x x y y 解：取步长h=0.1, x k =kh (k=0,1,?,6)。各方法计算结果及绝对误差见表(1)、(2)、(3)。 表（1） x n 欧拉公式 改进的欧拉公式 四阶标准龙格—库塔公式 y n 误差 y n 误差 y n 误差 0.0 1.000000 0 1 0 1 0 0.1 1.000000 1.0050000 1.00483750 0.2 1.010000 1.019025 1.01873090 0.3 1.029000 1.041218 1.04081842 0.4 1.056100 1.070802 1.07032029 0.5 1.090490 1.107076 1.10653093 0.6 1.131441 1.149404 1.14881193 表（2） 四阶阿当姆斯公式 n x n 显示公式隐式公式① y n误差y n误差 3 0.3 取自准确解— 1.04081801 4 0.4 1.07032292 1.07031966 5 0.5 1.10653548 1.10653041 6 0.6 1.14881481 1.14881101 备注y1, y2, y3 取自准确解y1, y2 取自准确解 ①本题关于y为线性，代入隐式公式，可解出y n+1，因此可直接求隐式公式之解。 表（3）四阶阿当姆斯预测—校正公式 n x n 显示公式隐式公式① Y n误差y n误差 4 0.4 1.07031992 1.07032014 5 0.5 1.10653027 1.10653077 6 0.6 1.14881103 1.14881175 备注y1, y2, y3由四阶标准龙格—库塔公式提供 对比以上各表数据知，在相同步长下求解同一问题时，方法阶数越高，解的精度也越高，一阶的欧拉公式精度最低，四阶标准龙格—库塔公式的精度大大高于改进的欧拉公式。四阶阿当姆斯方法、显式的精度略低于隐式。同是阿当姆斯预测—校正系统，带误差修正的精度又高于不带误差修正的。 从计算量看，四阶龙格—库塔法计算量最大，每前进一步要计算4次函数值f，而带误差修正的阿当姆斯预测—校正法与它精度差不多的，每前进一步只计算两次函数值f，所以后者是可取的。但它是四步法，必须用其它方法提供出始值，程序略复杂些。 例2 (步长的计算结果的影响)用欧拉公式求下述初值问题在x = 1处的近似解，并与准确解y (1)=1比较。 解分别取步长h=10-1、h=10-2、h=10-3、h=10-4计算，结果见下表。由表中数据可见， 当时结果完全失真，而取比计算量增加了十倍，但解 的精度却基本一样，可见取太浪费计算量了。 步长h y (1) 的计算值 上述结果差异很大的原因在于欧拉公式的绝对稳定区间为（－2, 0）步长h 应满足 )0,2(-∈??h y f ，对本题，x x y y x f 2)(1000 ),(2+--= 1000-=??y f ，故应取h 满 足 即 。可见取 时欧拉公式是数值不稳 定的，导致结果失真，而取 和 都满足稳定性要求，可用于求解。 此例说明，求解微分方程数值解一定选取注意步长，过大则导致解的失真，过小又会使计算量大增。究竟取多大步长才合适，不仅取决于所采用的数值方法，还决定于待解微分方程本身的特性。 例3：(Matlab 不同命令之差异)取h=0.2分别用不同的Matlab 命令解方程 ? ? ?=='11y 2t/y -y y )( (0 >> [t,y]=ode45(odefun,[0,4],1);n=length(t);e=sqrt(sum((sqrt(1+2*t)-y).^2/n));[n,e] ans =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,e] ans =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,e] ans =18.00000000000000 0.00973277413962 >> [t,y]=ode23t(odefun,[0,4],1);n=length(t);e=sqrt(sum((sqrt(1+2*t)-y).^2/n));[n,e] ans =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,e] ans =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,e] ans =15.00000000000000 0.24308870757312 >> [t,y]=ode15s(odefun,[0,4],1);n=length(t);e=sqrt(sum((sqrt(1+2*t)-y).^2/n));[n,e] ans =22.00000000000000 0.45509043060378 可见ode45精度较高，计算量较大，ode23计算量小，但误差大，ode113适中，用刚性方程组揭发解非刚性问题不适合，特别ode23s 计算量大且但误差大。 例4：解刚性方程组： ? ??=='=-='10z -100z,z 2 1y 99.99z,-01y 0y )()(. 其准确解为:x x x e z e e y 10010001.0,---=+=就h=0.01、h=0.02使用改进欧拉格式及方程组，并在一张 图上绘出其图形。 function[x,y]=naeuler2s(dyfun,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; end x=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=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(-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,?用欧拉显式方法求解：用欧拉法 解方程? ??=='11y xy y 1/3 )( (0 计算到(2)y ，并与精确解3/2 2 ()(2)/3y x x ??=+?? 比较。 MATLAB 程序 t=1; y=1;h=0.1 for k=1:10; y=y+h*t*y^(1/3) t=t+h end 运行结果 h y t 0.1 0.01 0.001 0.0001 y=((t 2+2)/3)3/2 ( 准确解) 1 1 1 1 1 1 1.1 1.1 1.1175 1.1068 1.1068 1.10681660630838 1.2 1.2136 1.2393 1.2277 1.2279 1.22787959356620 1.3 1.3416 1.3762 1.3639 1.3641 1.36413599028836 1.4 1.4849 1.5294 1.5162 1.5165 1.51656453868604 1.5 1.6447 1.6998 1.6857 1.6861 1.68617060118373 1.6 1.8217 1.8883 1.8734 1.8739 1.87398185690257 1.7 2.017 2.0962 2.0804 2.0810 2.08104468957300 1.8 2.2319 2.3243 2.3076 2.3083 2.30842116960842 1.9 2.4671 2.5739 2.5563 2.5571 2.55718653993016 2.0 2.7239 2.8177 2.8274 2.8283 2.82842712474619 欧拉显示方法是收敛的。 例7：实验目的：比较不同算法用于“非光滑”解时的结果。 内容：用欧拉显式方法和2阶龙格?库塔算法求解：sin ,(0)1,y kx y '== 计算(2)y π并在一张图上画出解曲线。 k=6时的程序如下： clear t=0;y=1;h=pi/60; for k=0.1:0.1:2; y=y+h*abs(sin(6*t)) t=t+h plot(t,y,'o'), hold on, end t0=0;tf=1/3*pi;x0=[1];k=pi/60 [t,x]=ode23('ex545f',[t0:k:tf],1) plot(t,x), hold on, for i=0:pi/60:pi/3 if i<=pi/6; y1=7/6-1/6*cos(6*i) else y1=9/6+1/6*cos(6*i) end plot(i,y1,'*'),hold on, end 例8 实验目的：观察欧拉显式方法的数值不稳定性。 内容：取h= 0.05, 0.04, 0.03, ?, 用欧拉显式方法求解： '=-=并画出曲线，同样用欧拉隐式方法重复求解上述问题。y y y 50,(0)100, 解：欧拉显式方法程序略。 '是y的线性函数故用欧拉隐式方法可直接解出y。程序如下 因50y y- = clear; h=0.01; t=0;y=100; for k=0:h:1; y=y/(1+50*h); t=t+h; plot(t,y,'b *'), hold on, end gtext('隐t=0.01') 运行结果： 图中标示为“隐00x ”为隐式欧拉格式结果，其余是欧拉格式的运行结果，可见隐式欧拉格式较欧拉格式稳定。 例9：解2阶微分方程?? ???? ?==π-=+....,)(0 y 00y 2t 1y y 2设初始时间t0=0，终止时间tf=3π， 解：将方程改写为：?????π -+-==2t 1x x x x 2 122 1.. ， 写为矩阵表达式：??? ? ??=???? ?????? ??π-???? ??+???? ?????? ? ?-=???? ??000x 0x 2t 110x x 0110x x 2122121)()(;..。 本题解析解：2 2 2 t cont 12 1y ππ- -+ =))(( MATLAB 程序 function xdot=ex545f(t,x) u=1-(t.^2)/(pi^2); xdot=[0 1;-1 0]*x+[0 1]'*u; %建立ex545f.m 函数程序。 clf,t0=0;tf=3*pi;x0t=[0;0]; %给出初始值。 [t,x]=ode23('ex545f',[t0,tf],x0t) %调用MATLAB 中的数值积分2阶3阶龙—库公式。 y=x(:,1); %y 是x 的第2列。 for I=1:length(t); y2(I)=(1+2/(pi^2))*(1-cos(t(I)))-t(I)^2/(pi^2) end %在数值积分的序列点上计算解析解的值。 u=1-(t.^2)/(pi^2); clf,plot(t,y,'-',t,u,'+',t,y2,'o') % 绘出数值解、解析解的图形。 legend('数值解','输入量','解析解') 运行结果解析解与数值解看不出差别，是因为本数值积分是按默认精度10-3,如果用长格式输出y 、y2可看出它们的微小误差。 例10 解微分方程组. ?? ?? ??? ===-=-==1 )0(,1)0(,0)0(51.0'''321213312321y y y y y y y y y y y y 解 1、建立m-文件rigid.m 如下： function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2); 2、取t0=0，tf=12，输入命令： >> [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') 3、结果如图. 【例11】采用最简洁格式的ODE 文件和解算指令，研究围绕地球旋转的卫星轨道。 （1）问题的形成 轨道上运动的卫星，在Newton 第二定律 22dt r d m ma F == ，和万有引力定律 程序见数学建模高教出版社数学建模与 数学实验To Matlab （ff5） r r m M G F E 3-=作用下，有r r M G dt r d a E 322-==。即3r x GM a E x -=，3r y GM a E y -=，而 22y x r += 。这里)/(10672.62211kg m N G ??=-是引力常数，)(1097.524kg M E ?=是地球 的质量。又假定卫星以初速度)/(4000)0(s m v y =在)(102.4)0(7m x ?-=处进入轨道。 （2）构成一阶微分方程组 令[ ][ ][]T T y x T y x y x v v y x y y y y Y ''===43 21 ，则 ? ??? ????? ??? ??? ?? ?+?-+?-=????????????''''='2 /322 2 2/322 14 34 3 21 )()()(y x y GM y x y GM y y y y y y t Y E E (5.14.2.1-5) 初始条件为 [] T y v x Y )0(00)0()0(= (5.14.2.1-6) （1）根据式(5.14.2.1-5)编写最简洁格式的ODE 文件 [DYDt2.m] function Yd=DYDt2(t,Y,flag,G,ME) % flag 按ODE 文件格式规定，必须是第三输入宗量。对它的赋值由ode45指令自动产生。 第4、5宗量是被传递的参数 switch flag case '' %按规定：这里必须是空串。在此为"真"时，完成以下导数计算。 X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V;-G*ME*X/r^3]; otherwise error(['Unknown flag ''' flag '''.']); end （2）对微分方程进行解算 G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e 7;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 on axis('image') %保证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.m] function varargout=DYDt3(t,Y,flag,G,ME,tspan,Y0) % DYDt3.m 供主程序调用的ODE函数文件 % 本文件自带三个子函数：f,fi,fev。 % t, Y 分别是自变量和一阶函数向量，是最基本的输入宗量。 % flag 第三输入宗量，它专供解算指令（如ode45）作调用通知。 % 在运行中，解算指令会根据需要向flag发不同的字符串。 % varargout 是"变"输出宗量。它由变维的元胞数组构成。每个元胞中可以存放指令 % 所产生的任意形式的数据。 switch flag case '' % 必须用空串符。情况为"真"时，计算导数 dY/dt = f(t,Y)。 varargout{1} = f(t,Y,G,ME); %输出为一个元胞，容纳f子函数的一个输出Yd case 'init' % 必须用'init'，情况为"真"时，传送计算区间、初值、设置参数。 [varargout{1:3}] = fi(tspan,Y0); %输出为三个元胞，容纳fi子函数的三个输出量 case 'events' % 必须用'events'，情况为"真"时，设置事件性质。 [varargout{1: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/r^3]; % ------------------------------------------------------------------ function [ts,y0,options] = fi(tspan,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 = [dDSQdt; dDSQdt]; %定义两个穿越0的事件 direction = [1; -1]; %第一事件：以渐增方式穿越0。第二事件：以渐减方式穿越0 isterminal = [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); % <3> X=YY(:,1);Y=YY(:,2); plot(X,Y,'b','Linewidth',2);hold on text(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,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;Y E=RE*YE;ZE=0*ZE; mesh(XE,YE,ZE) text(1e7,1e7,'地球','Color','r'), hold off %产生红色文字注释 例13 范例：状态转移方程组模型 计算机网络可靠性分析问题随着计算机通信网络系统特别是Internet 网络日益广泛的应用，提高系统的可靠性意义重大．研究和分析具有实用性的高可靠计算机通信网络系统，是国际上非常活跃的一个研究方向． （1） 问题及假设 我们将研究的计算机通信网络系统称为无冗余的防火墙协议系统．计算机随时可能发生三个事件一无故障、间歇故障和永久故障．因此，计算机一般处于三种工作状态： 无故障工作、带故障工作和不工作．这三种状态之间的转移过程如图3．9所示．要求建立该系统的状态转移模型，并进行可靠性分析． (2) 分析与模型 该问题属于状态转移问题，利用马尔科夫状态转移原理，用)(1t P ,)(2t P 和)(3t P 分别表示系统处于无故障工作、带故障工作和不工作三种状态的概率，则有以下状态转移方程组： ???? ?? ? ? ???+++=++-=++-=+++-=)()()()2/() ()()]([)()2/() ()()()()()(]2/)2/[() (21321221211t P t P dt t dP t P t P dt t dP t P t P t P t P dt t dP t p t p t p t t p t t p λλλλλλγλγλλγλλλ 初始条件]0,0,1[)]0(),0(),0([321=P P P ,参数取值1.0~01.0:,10~10:,10~10:3445γλλ----t p 这是一个带参数的微分方程组模型． (3) 模型求解 求解这个带参数的微分方程组模型有两种方法，一是用特征根法求解析解，另一个是用数值解法求数值解．分别求解如下： 假定模型中的参数取下限值，即 01.0,10,104 5===--γλλt p ①.解析解法 利用常系数线性微分方程组的特征根法． 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); （0．7053 -0．7089 0．0035）； （-0．7053 -0．0035 0．7089） 根据特征值和特征向量写出其解析解,利用初始条件]0,0,1[)]0(),0(),0([321=P P P 得到特解 t t e e t P 0001.00102.0199503724.0004937.0)(--+= t t e e t P 000.00102.020049378.00049623.0)(--+-= t t e e t P 0001.00102.0300011612.10000245.01)(---+= ②．数值解法 MATLAB 程序如下 首先编写m 文件（eqs0．m ） function xdot=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(:,)); xlabe l(‘时间t(小时)’);ylabel(‘可靠度R(t)’); 050010001500200025003000-2.5 -2-1.5 -1-0.5 0.5 11.5 2title(‘参数取值lp=0.00001;lt=0.00001;gm=0.01’);grid 输出结果如图3.10 计算表明：当t=1000小时情况下，].,.,.[)](),(),([058400047093690t p t p t p 321= 显示系统工作的概率为94.18%，从图3.10中还可以看出系统可靠性变化更加细致的情况 , 何时 可靠度达到99%，何时达到98%等等。 ③.敏感性分析 当参数取值在以下范围 1.0~01.0:,10~10:,10~10:3445γλλ----t p 内变化时，系统的可靠性情况如何呢？留给读者自己分析 例14 ?? ???===---0)0(';2)0(0 )1(1000222x x x dt dx x dt x d 解: 令 y 1=x ，y 2=y 1’ 则微分方程变为一阶微分方程组： ?? ???==--==0 )0(,2)0()1(1000 ''211221221y y y y y y y y 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',[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))，乙舰位于),1(0t v Q 由于导弹头始终对准乙舰，故此时直线PQ 就是导弹的轨迹曲线弧OP 在点P 处的切线， 即有 x y t v y --=1'0 即 y y x t v +-=')1(0 (1) 又根据题意，弧OP 的长度为AQ 的5倍， 即 t v dx y x 00 25'1= +? (2) 由(1),(2)消去t 整理得模型: 21 (1)"1' (3)5 x y y -=+ 初值条件为: 0)0(=y 0)0('=y 解即为导弹的运行轨迹: 24 5 )1(125)1(8556 54+-+--=x x y 当1=x 时245=y ，即当乙舰航行到点)24 5 ,1(处时被导弹击中. 被击中时间为:0 0245 v v y t = =. 若v 0=1, 则在t=0.21处被击中. 轨迹图见程序chase1 解法二(数值解) 令y1=y,y2=y1’，将方程（3）化为一阶微分方程组。 21(1)''15x y y -=+ ? 12 2 21'1'1/(1)5y y y y x =?? ?=+-?? 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，xf=0.9999，建立主程序ff6.m 如下： x0=0，xf=0.9999 [x,y]=ode15s('eq1',[x0 xf],[0 0]); plot(x,y(:,1),’b.') hold on y=0:0.01:2; plot(1,y,’b*') 结论: 导弹大致在（1，0.2）处击中乙舰 程序数学建模高教出版社数学建模与数学实验To Matlab(ff6) 六.习题 1. 实验目的：利用实际例子筛选有效算法。 内容：将不同算法用于下面的方程，选出你最满意的方法。也可以自己设计新方法，看是否比我们所讲方法好。 1).22()/,(1)1y y xy x y '=+=, 计算到x=2. 2). 21 4,(2)12 y x y x y '=+-=-, 计算到x=4. 3). 1(sin ),(2/)/2y y y y x y ππ-'=-=, 计算到x=2. 4).sin ,(1)1y y x y '==, 计算到x=2. 提示：先确定“满意”的尺度和标准，可以比较算法的效率、精确度等。 2. 实验目的：检查各种算法的长期行为。 内容：给定方程组()(),()(),(0)0,(0). x t ay t y t bx t x y b '=??'=-? ?=??=?的解是x ?y 平面上一个椭圆，利用你已经知道的算法，取足够 小的步长，计算上述方程的轨道，看那种算法能够保持椭圆轨道不变。(计算的时间步要足够多)。 3.实验目的：初步认识刚性微分方程。 内容：任选一种你知道的显式方法，取不同步长，求解 第8章 常微分方程边值问题的数值解法 引 言 第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。 只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为 则边值问题(8.1.1)有唯一解。 推论 若线性边值问题 ()()()()()(),, (),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤?? ==? (8.1.2) 满足 (1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。 求边值问题的近似解，有三类基本方法： (1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解； (2) 有限元法(finite element method)； (3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。 差分法 8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法 设二阶线性常微分方程的边值问题为 (8.2.1)(8.2.2) ()()()(),,(),(), y x q x y x f x a x b y a y b αβ''-=< 第 一 章 一阶微分方程的解法的小结 ⑴、可分离变量的方程： ①、形如 )()(y g x f dx dy = 当0)(≠y g 时，得到 dx x f y g dy )() (=，两边积分即可得到结果； 当0)(0=ηg 时，则0)(η=x y 也是方程的解。 例1.1、 xy dx dy = 解：当0≠y 时，有 xdx y dy =，两边积分得到)(2ln 2为常数C C x y += 所以)(112 12 C x e C C e C y ±==为非零常数且 0=y 显然是原方程的解； 综上所述，原方程的解为)(12 12 为常数C e C y x = ②、形如0)()()()(=+dy y Q x P dx y N x M 当0)()(≠y N x P 时，可有 dy y N y Q dx x P x M ) () ()()(=，两边积分可得结果； 当0)(0=y N 时，0y y =为原方程的解，当0(0=） x P 时，0x x =为原方程的解。 例1.2、0)1()1(2 2 =-+-dy x y dx y x 解：当0)1)(1(2 2 ≠--y x 时，有 dx x x dy y y 1 122-=-两边积分得到 )0(ln 1ln 1ln 22≠=-+-C C y x ，所以有)0()1)(1(22≠=--C C y x ； 当0)1)(1(2 2 =--y x 时，也是原方程的解； 综上所述，原方程的解为)()1)(1(2 2 为常数C C y x =--。 ⑵可化为变量可分离方程的方程： ①、形如 )(x y g dx dy = 解法：令x y u =，则udx xdu dy +=，代入得到)(u g u dx du x =+为变量可分离方程，得到 常微分方程解题方法总结 来源：文都教育 复习过半， 课本上的知识点相信大部分考生已经学习过一遍 . 接下来， 如何将零散的知 识点有机地结合起来， 而不容易遗忘是大多数考生面临的问题 . 为了加强记忆， 使知识自成 体系，建议将知识点进行分类系统总结 . 著名数学家华罗庚的读书方法值得借鉴， 他强调读 书要“由薄到厚、由厚到薄”，对同学们的复习尤为重要 . 以常微分方程为例， 本部分内容涉及可分离变量、 一阶齐次、 一阶非齐次、 全微分方程、 高阶线性微分方程等内容， 在看完这部分内容会发现要掌握的解题方法太多， 遇到具体的题 目不知该如何下手， 这种情况往往是因为没有很好地总结和归纳解题方法 . 下面以表格的形 式将常微分方程中的解题方法加以总结，一目了然，便于记忆和查询 . 常微分方程 通解公式或解法 ( 名称、形式 ) 当 g( y) 0 时，得到 dy f (x)dx ， g( y) 可分离变量的方程 dy f ( x) g( y) 两边积分即可得到结果； dx 当 g( 0 ) 0 时，则 y( x) 0 也是方程的 解 . 解法：令 u y xdu udx ，代入 ，则 dy 齐次微分方程 dy g( y ) x dx x u g (u) 化为可分离变量方程 得到 x du dx 一 阶 线 性 微 分 方 程 P ( x)dx P ( x) dx dy Q(x) y ( e Q( x)dx C )e P( x) y dx 伯努利方程 解法：令 u y1 n，有 du (1 n) y n dy ， dy P( x) y Q( x) y n（n≠0,1）代入得到du (1 n) P(x)u (1 n)Q(x) dx dx 求解特征方程：2 pq 三种情况： 二阶常系数齐次线性微分方程 y p x y q x y0 二阶常系数非齐次线性微分方程 y p x y q x y f ( x) (1)两个不等实根：1, 2 通解： y c1 e 1x c2 e 2x (2) 两个相等实根：1 2 通解： y c1 c2 x e x (3) 一对共轭复根：i , 通解： y e x c1 cos x c2 sin x 通解为 y p x y q x y 0 的通解与 y p x y q x y f ( x) 的特解之和. 常见的 f (x) 有两种情况： x （ 1）f ( x)e P m ( x) 若不是特征方程的根，令特解 y Q m ( x)e x；若是特征方程的单根，令特 解 y xQ m ( x)e x；若是特征方程的重根， 令特解 y*x2Q m (x)e x； （2）f (x) e x[ P m ( x) cos x p n ( x)sin x] 淮北师范大学 2013届学士学位论文 常微分方程数值解法的误差分析 学院、专业数学科学学院数学与应用数学 研究方向计算数学 学生姓名李娜 学号 20091101070 指导教师姓名陈昊 指导教师职称讲师 年月日 常微分方程数值解法的误差分析 李娜 （淮北师范大学数学科学学院，淮北，235000） 摘要 自然界与工程技术中的很多现象，往往归结为常微分方程定解问题。许多偏微分方程问题也可以化为常微分方程问题来近似求解。因此，研究常微分方程的数值解法是有实际应用意义的。数值解法是一种离散化的数学方法，可以求出函数的精确解在自变量一系列离散点处的近似值。随着计算机计算能力的增强以及数值计算方法的发展，常微分方程的数值求解方法越来越多，比较成熟的有Euler 法、后退Euler法、梯形方法、Runge—Kutta方法、投影法和多步法，等等．本文将对这些解的误差进行分析,以求能够得到求解常微分数值解的精度更好的方法。 关键词：常微分方程, 数值解法, 单步法, 线性多步法, 局部截断误差 Error Analysis of Numerical Method for Solving the Ordinary Differential Equation Li Na (School of Mathematical Science, Huaibei Normal University, Huaibei, 235000) Abstract In nature and engineering have many phenomena , definite solution of the problem often boils down to ordinary differential equations. So study the numerical solution of ordinary differential equations is practical significance. The numerical method is a discrete mathematical methods, and exact solution of the function can be obtained in the approximation of a series of discrete points of the argument.With the enhanced computing power and the development of numerical methods,ordinary differential equations have more and more numerical solution,there are some mature methods. Such as Euler method, backward Euler method, trapezoidal method, Runge-Kutta method, projection method and multi-step method and so on.Therefore, numerical solution of differential equation is of great practical significance. Through this paper, error of these solutions will be analyzed in order to get a the accuracy better way to solve the numerical solution of ordinary differential. Keywords:Ordinary differential equations, numerical solution methods, s ingle ste p methods, l inear multi-step methods, local truncation error i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。 常微分方程初值问题数值解法 朱欲辉 （浙江海洋学院数理信息学院, 浙江舟山316004） [摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法.然而在生产实 际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的Euler法、Runge－Kutta法以及线性多步法中的Adams显隐式公式和预测校正 公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法计算出近似解, 从得出的结果对比各种方法的优缺点. [关键词]：常微分方程;初值问题; 数值方法; 收敛性; 稳定性; 误差估计 Numerical Method for Initial-Value Problems Zhu Yuhui (School of Mathematics, Physics, and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang 316004) [Abstract]:In the course about ordinary differential equations, the methods for analytic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can’t be e xpressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis. [Keywords]:Ordinary differential equation; Initial-value problem; Numerical method; Convergence; Stability;Error estimate 各类微分方程的解法 1.可分离变量的微分方程解法 一般形式:g(y)dy=f(x)dx 直接解得∫g(y)dy=∫f(x)dx 设g(y)及f(x)的原函数依次为G(y)及F(x),则G(y)=F(x)+C为微分方程的隐式通解 2.齐次方程解法 一般形式:dy/dx=φ(y/x) 令u=y/x则y=xu,dy/dx=u+xdu/dx,所以u+xdu/dx=φ(u),即du/［φ(u)-u］=dx/x 两端积分,得∫du/［φ(u)-u］=∫dx/x 最后用y/x代替u,便得所给齐次方程的通解 3.一阶线性微分方程解法 一般形式:dy/dx+P(x)y=Q(x) 先令Q(x)=0则dy/dx+P(x)y=0解得y=Ce－ ∫P(x)dx,再令y=u e－∫P(x)dx代入原方程解得u=∫Q(x) e∫P(x)dx dx+C,所以y=e－∫P(x)dx［∫Q(x)e∫P(x)dx dx+C］ 即y=Ce－∫P(x)dx +e－ ∫P(x)dx∫Q(x)e∫P(x)dx dx为一阶线性微分方程的通解 4.可降阶的高阶微分方程解法 ①y(n)=f(x)型的微分方程 y(n)=f(x) y(n-1)= ∫f(x)dx+C1 y(n-2)= ∫［∫f(x)dx+C1］dx+C2 依次类推,接连积分n次,便得方程y(n)=f(x)的含有n个任意常数的通解②y”=f(x,y’) 型的微分方程 令y’=p则y”=p’,所以p’=f(x,p),再求解得p=φ(x,C1) 即dy/dx=φ(x,C1),所以y=∫φ(x,C1)dx+C2 ③y”=f(y,y’) 型的微分方程 令y ’=p 则y ”=pdp/dy,所以pdp/dy=f(y,p),再求解得p=φ(y,C 1) 即dy/dx=φ(y,C 1),即dy/φ(y,C 1)=dx,所以∫dy/φ(y,C 1)=x+C 2 5.二阶常系数齐次线性微分方程解法 一般形式:y ”+py ’+qy=0,特征方程r 2+pr+q=0 6.二阶常系数非齐次线性微分方程解法 一般形式: y ”+py ’+qy=f(x) 先求y ”+py ’+qy=0的通解y 0(x),再求y ”+py ’+qy=f(x)的一个特解y*(x) 则y(x)=y 0(x)+y*(x)即为微分方程y ”+py ’+qy=f(x)的通解 求y ”+py ’+qy=f(x)特解的方法: ① f(x)=P m (x)e λx 型 令y*=x k Q m (x)e λx ［k 按λ不是特征方程的根,是特征方程的单根或特征方程的重根依次取0,1或2］再代入原方程,确定Q m (x)的m+1个系数 ② f(x)=e λx ［P ｌ(x)cos ωx+P n (x)sin ωx ］型 令y*=x k e λx ［Q m (x)cos ωx+R m (x)sin ωx ］［m=max ﹛ｌ,n ﹜,k 按λ+i ω不是特征方程的根或是特征方程的单根依次取0或1］再代入原方程,分别确定Q m (x)和R m (x)的m+1个系数 第八章 常微分方程数值解法 考核知识点： 欧拉法，改进欧拉法，龙格-库塔法，单步法的收敛性与稳定性。 考核要求： 1． 解欧拉法，改进欧拉法的基本思想；熟练掌握用欧拉法，改进欧拉法、求微 分方程近似解的方法。 2． 了解龙格-库塔法的基本思想；掌握用龙格-库塔法求微分方程近似解的方 法。 3． 了解单步法的收敛性、稳定性与绝对稳定性。 例1 用欧拉法，预估——校正法求一阶微分方程初值问题 ? ??=-='1)0(y y x y ，在0=x （0，1）0.2近似解 解 （1）用1.0=h 欧拉法计算公式 n n n n n n x y y x y y 1.09.0)(1.01+=-+=+，1.0=n 计算得 9.01=y 82.01.01.09.09.02=?+?=y （2）用预估——校正法计算公式 1,0)(05.01.09.0)0(111)0(1=???-+-+=+=++++n y x y x y y x y y n n n n n n n n n 计算得 91.01=y ，83805.02=y 例2 已知一阶初值问题 ???=-='1 )0(5y y y 求使欧拉法绝对稳定的步长n 值。 解 由欧拉法公式 n n n n y h y h y y )51(51-=-=+ n n y h y ~)51(~1-=+ 相减得01)51()51(e h e h e n n n -==-=-Λ 当 151≤-h 时，4.00≤ 第八章 常微分方程的数值解法 一．内容要点 考虑一阶常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节 点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。 用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。 (一)常微分方程处置问题解得存在唯一性定理 对于常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 如果： (1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。 (2) ),(y x f 对于y 满足利普希茨条件，即 2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程?????==0 0)() ,(y x y y x f dx dy 的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。 定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。 收敛性定理：若一步方法满足： (1)是p 解的. (2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件. (3) 初始值y 0是精确的。则),()()(p h O x y kh y =-kh =x -x 0,也就是有 0x y y lim k x x kh 0h 0 =--=→)( (一)、主要算法 1．局部截断误差 局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~ +k y 的误差y (x k+1)- 1~ +k y 称为局部截断误差。 注意：y k+1和1~ +k y 的区别。因而局部截断误差与误差e k +1=y (x k +1) -y k +1不同。 如果局部截断误差是O (h p+1),我们就说该数值方法具有p 阶精度。 常微分方程数值解法 【作用】微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步： 1. 根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。 2. 找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。 3. 运用这些规律列出方程和定解条件。基本模型 1. 发射卫星为什么用三级火箭 2. 人口模型 3. 战争模型 4. 放射性废料的处理通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来” 的于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 1. 改进Euler 法： 2. 龙格—库塔( Runge—Kutta )方法： 【源程序】 1. 改进Euler 法： function [x,y]=eulerpro(fun,x0,x1,y0,n);%fun 为函数，(xO, x1)为x 区间，yO 为初始值，n 为子 区间个数 if nargin<5,n=5O;end h=(x1-xO)/n; x(1)=xO;y(1)=yO; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(fun,x(i),y(i)); y2=y(i)+h*feval(fun,x(i+1),y1); y(i+1)=(y1+y2)/2; end 调用command 窗口 f=i nlin e('-2*y+2*x A2+2*x') [x,y]=eulerpro(f,O,,1,1O) 2 x +2x , (0 < x < , y(0) = 1 求解函数y'=-2y+2 2. 龙格—库塔( Runge—Kutta )方法： [t,y]=solver('F',tspan ，y0) 这里solver为ode45, ode23, ode113，输入参数F是用M文件定义的微分方程y'= f (x, y)右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。 注：ode45和ode23变步长的，采用Runge-Kutta算法。 ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差，是一种自适应步长(变步长)的常微分方程数值解法，其整体截断误差为(△ 口人5解 决的是Nonstiff(非刚性)常微分方程。 实验4 常微分方程的数值解 【实验目的】 1．掌握用MATLAB软件求微分方程初值问题数值解的方法； 2．通过实例用微分方程模型解决简化的实际问题； 3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。 【实验内容】 题3 小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。 模型及其求解 火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。 在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式： a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8 dh/dt=v 又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下： function [ dx ] = rocket( t,x ) a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8; dx=[x(2);a]; end ts=0:1:60; x0=[0,0]; [t,x]=ode45(@rocket,ts,x0); h=x(:,1); v=x(:,2); a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a]; 数据如下： t h v a 0 0 0 13.06 1.00 6.57 13.19 13.30 2.00 26.44 26.58 1 3.45 3.00 59.76 40.06 13.50 4.00 106.57 53.54 13.43 5.00 16 6.79 66.89 13.26 6.00 240.27 80.02 12.99 7.00 326.72 92.83 12.61 8.00 425.79 105.22 12.15 9.00 536.99 117.11 11.62 10.00 659.80 128.43 11.02 11.00 793.63 139.14 10.38 12.00 937.85 149.18 9.71 13.00 1091.79 158.55 9.02 14.00 1254.71 167.23 8.33 15.00 1425.93 175.22 7.65 16.00 1604.83 182.55 6.99 17.00 1790.78 189.22 6.36 18.00 1983.13 195.27 5.76 19.00 2181.24 200.75 5.21 20.00 2384.47 205.70 4.69 21.00 2592.36 210.18 4.22 22.00 2804.52 214.19 3.79 23.00 3020.56 217.79 3.41 24.00 3240.08 221.01 3.07 25.00 3462.65 223.92 2.77 26.00 3687.88 226.56 2.50 27.00 3915.58 228.97 2.27 -293- 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如 22x y dx dy +=。于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 §1 常微分方程的离散化 下面主要讨论一阶常微分方程的初值问题，其一般形式是 ?????=≤≤=0 )() ,(y a y b x a y x f dx dy （1） 在下面的讨论中，我们总假定函数),(y x f 连续，且关于y 满足李普希兹(Lipschitz)条 件，即存在常数L ，使得 |||),(),(|y y L y x f y x f ?≤? 这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。 所谓数值解法，就是求问题（1）的解)(x y 在若干点 b x x x x a N =<<<<=L 210 处的近似值),,2,1(N n y n L =的方法，),,2,1(N n y n L =称为问题（1）的数值解， n n n x x h ?=+1称为由n x 到1+n x 的步长。今后如无特别说明，我们总取步长为常量h 。 建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： （i ）用差商近似导数 若用向前差商h x y x y n n ) ()(1?+代替)('n x y 代入（1）中的微分方程，则得 )1,,1,0())(,() ()(1?=≈?+N n x y x f h x y x y n n n n L 化简得 ))(,()()(1n n n n x y x hf x y x y +≈+ 如果用)(n x y 的近似值n y 代入上式右端，所得结果作为)(1+n x y 的近似值，记为1+n y ， 则有 )1,,1,0() ,(1?=+=+N n y x hf y y n n n n L （2） 这样，问题（1）的近似解可通过求解下述问题 ?? ?=?=+=+) () 1,,1,0(),(01a y y N n y x hf y y n n n n L （3） 得到，按式（3）由初值0y 可逐次算出N y y y ,,,21L 。式（3）是个离散化的问题，称为差分方程初值问题。 i.常微分方程初值问题数值解法 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法--差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<<<<=L （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-L 方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t L 上的差分解1,,N u u L 。 二阶微分方程： 时为非齐次 时为齐次，0)(0)()()()(22≠≡=++x f x f x f y x Q dx dy x P dx y d 二阶常系数齐次线性微分方程及其解法： 2 122,)(2,,(*)0)(1,0(*)r r y y y r r q pr r q p qy y p y 式的两个根、求出的系数； 式中的系数及常数项恰好是，，其中、写出特征方程：求解步骤： 为常数； ，其中?'''=++?=+'+''式的通解：出的不同情况，按下表写、根据(*),321r r 二阶常系数非齐次线性微分方程 型 为常数； 型，为常数 ，]sin )(cos )([)()()(,)(x x P x x P e x f x P e x f q p x f qy y p y n l x m x ωωλλλ+===+'+'' 二阶常系数非齐次线性微分方程的一般形式是 ''+'+=y py qy f x () （1） 其中p q ,是常数。 方程（1）的通解为对应的齐次方程 0=+'+''qy y p y （2） 的通解Y 和方程（1）的一个特解*y 之和。即 *y Y y +=.我们已解决了求二阶常系数齐 次线性方程通解的问题，所以，我们只需讨论求二阶常系数非齐次线性微分方程的特解* y 的方法。 下面我们只介绍当方程（1）中的)(x f 为如下两种常见形式时求其特解*y 的方法。 一、 f x e P x x m ()()=?λ型 由于方程（1）右端函数f x ()是指数函数e x λ?与m 次多项式P x m ()的乘积，而指数 函数与多项式的乘积的导数仍是这类函数，因此，我们推测： 方程（1）的特解应为 y e Q x x *?=λ()( Q x ()是某个次数待定的多项式 ) y e Q x e Q x x x *??'=+'λλλ()() y e Q x Q x Q x x *?"=?+'+''λλλ[()()()]22 代入方程（1），得 e Q x p Q x p q Q x e P x x x m λλλλλ???''++'+++≡?[()()()()()]()22 消去e x λ?，得 ''++'+++≡Q x p Q x p q Q x P x m ()()()()()()22λλλ （3） 讨论 01、如果λ不是特征方程 r pr q 20++=的根。 即 02≠++q p λλ 由于P x m ()是一个m 次的多项式，欲使（3）的两端恒等，那未Q x ()必为一个m 次多项式，设为 Q x b x b x b x b m m m m m ()=++++--0111Λ 将之代入（3），比较恒等式两端x 的同次幂的系数，就得到以b b b b m m 01 1,,,,Λ-为未知数的m +1个线性方程的联立方程组，解此方程组可得到这m +1个待定的系数，并得到特解 y e Q x x m *?=λ() 02、如果λ是特征方程 r pr q 20++=的单根。 即 λλ20++=p q ，但 20λ+≠p 欲使（3）式的两端恒等，那么'Q x ()必是一个m 次多项式。 因此，可令 Q x x Q x m ()()=? 并且用同样的方法来确定)(x Q 的系数b b b b m m 0 11,,,,Λ-。 03、如果λ是特征方程 r pr q 20++=的二重根。 即 λλ20++=p q ，且 20λ+=p 。 欲使（3）式的两端恒等，那么''Q x ()必是一个m 次多项式 因此， 可令 Q x x Q x m ()()=?2 并且用同样的方法来确定)(x Q 的系数b b b b m m 011,,,,Λ-。 目录 1 引言 (1) 2 二阶常系数常微分方程的几种解法 (1) 2.1 特征方程法 (1) 2.1.1 特征根是两个实根的情形 (2) 2.1.2 特征根有重根的情形 (2) 2.2 常数变异法 (4) 2.3 拉普拉斯变化法 (5) 3 常微分方程的简单应用 (6) 3.1 特征方程法 (7) 3.2 常数变异法 (9) 3.3 拉普拉斯变化法 (10) 4 总结及意义 (11) 参考文献 (12) 二阶常微分方程的解法及其应用 摘要：本文通过对特征方程法、常数变易法、拉普拉斯变换法这三种二阶常系数常微分方程解法进行介绍，特别是其中的特征方程法分为特征根是两个实根的情形和特征根有重根的情形这两种情况，分别使用特征值法、常数变异法以及拉普拉斯变换法来求动力学方程，现今对于二阶常微分方程解法的研究已经取得了不少成就，尤其在二阶常系数线性微分方程的求解问题方面卓有成效。应用常微分方程理论已经取得了很大的成就，但是，它的现有理论也还远远不能满足需要，还有待于进一步的发展，使这门学科的理论更加完善。 关键词：二阶常微分方程;特征分析法；常数变异法；拉普拉斯变换 METHODS FOR TWO ORDER ORDINARY DIFFERENTIAL EQUATION AND ITS APPLICATION Abstract:This paper introduces the solution of the characteristic equation method, the method of variation of parameters, the Laplasse transform method the three kind of two order ordinary differential equations with constant coefficients, especially the characteristic equation method which is characteristic of the root is the two of two real roots and characteristics of root root, branch and don't use eigenvalue method, method of variation of constants and Laplasse transform method to obtain the dynamic equation, the current studies on solution of ordinary differential equations of order two has made many achievements, especially in the aspect of solving the problem of two order linear differential equation with constant coefficients very fruitful. Application of the theory of ordinary differential equations has made great achievements, however, the existing theory it is still far from meeting the need, needs further development, to make the discipline theory more perfect. Keywords:second ord er ordinary differential equation; Characteristic analysis; constant variation method; Laplasse transform 1 引言 数学发展的历史告诉我们，300年来数学分析是数学的首要分支，而微分方程常微分方程边值问题的数值解法
一阶常微分方程解法总结
常微分方程解题方法总结.doc
常微分方程数值解法的误差分析教材
常微分方程数值解法
常微分方程初值问题数值解法
各类微分方程的解法大全
常微分方程数值解法
常微分方程数值解法
常微分方程数值解法
常微分方程的数值解
15第十五章 常微分方程的解法
(整理)常微分方程数值解法
(完整版)专题一(二阶常微分方程解法)
二阶常微分方程的解法及其应用