第15微分方程的解法

/第十五章 y2x2，于是对于用微分方程解决实际问题来说，数值解法就是一个分重要的 dyf(x,

ax在下面的讨论中，我们总假定函数f(x,y)连续，且关于y满足兹(Lipschitz)条件，即存在常数L，使得|f(x,y)f(x,y)|L|yy所谓数值解法，就是求问题（1）的解y(x)在若干点ax0x1x2LxNyn(n1,2,LN)yn(n1,2,LN)称为问题（1）yxn1yxn)y n

y(xn1)y(xn)h

f(

,y(

(ny(xn1)y(xn)hf(xn,y(xnyn1ynhf(xn,yn (n yn1ynhf(xn,yn (n 得到，按式（3）y0y1y2L。式（3）是个离散化的问题，称为差y(xn1)y(xn)

f(x,y( (n 右边的积分用矩形或梯形计算y(xn1)y(xn)hy'(xn)y(xn)hf(xn,y(xn散化的计算yn1ynhf(xn,yn计算。其中的Taylor展开法，不仅可以得到求数值解的，而且容易估计截断§2（Euler）方Euler即由（3）yxnyn(n1,2,L。这组求问题（1）的数值解称为向前Euler。

n

)y(xn1)y(xn)hyn1ynhf(xn1,yn1 (n 用这组求问题（1）的数值解称为向后EulerEulerEulerEuler y(0)yhf(x,ynn nnnnnnn

,y(k)

(k对于向前 （3）我们看到，当n1,2,L 右端的yn都是近似的yn1yn1y(xn)hf(xn,y(xn

)y(

)hy'(xn)

)

n1

2

即局部截断误差是h2 改进的Euler即xn1f(x,y(x))dxh[f(

,y(

))f(

n

,y(

n

yn

h[f(

,

)f(

n

y(0)yhf(x,y y(k1)y

h [f(x,y)f ,

|y(k1)y(k)|hL|y(k)y(k1)n n n

n其中L为Lipschitz常数。因此，当0 2出一种新的方法—改进Euler法。后用梯形校正求得近似值yn1，即yn1ynhf(xnyn yn1

h[f(x,y)f(x ,y )]

ypynhf(xn,ynyqynhf(xnh,yp yn1

(

yq

y(xn1)y(xn)y'( fxy

h),0yxn1yxn)hfxnh,yxnh)) Kfxnhyxnh，称为区间[xnxn1上的平均斜率。可见给出一种（13）EulerfxnynKEuler可理 K，就有可能构造出精度更高的计算。这就是—方法的基本思想yn1ynh(1k12k2k1f(xn,yn kf(xh,

hk1),0, 其中12,为待定系数，看看如何确定它们使（14）yxn1yn1ynyxn，所以（14）k1

f(xn,y(xn))y'(xnk2f(xnh,y(xn)hk1 f(x,y(

))hf(x,y(x

yn

y(

)2

(f

ff

y(xn1)y(xn)hy'(xn)2

3y''(xn)O(h中y'f，y''fff，可见为使误差y( ) O(h3)，只须 n n1，1 12

，1进的 。可以证明，在[xn,xn1]内只取2点 —精度最高为4阶yn1ynh(1k12k23k34k4 f(x,y k2f(xn1h,yn1hk1 f(xh,yhkhk3 3k4f(xn3h,yn4hk15hk26hk3其中待定系数i,i,i共13个，经过与推导2阶 单的一组i,iiyyn

h

2k22k3k4k1k2

f(xh,

kk f(

h,

)2)k4f(xnh,ynhk3 前一步的值yn，单步法的一般形式是yn1ynh(xn,yn, (n0,1,L,N 其中xyhEulerfxyEuler法的(x,y,h)1[f(x,y)f(xh,yhf(x,y))]计算yn1，这就是多步法的基本思想。经常使用的是线性多步法。ny(xn1)y(xn) f(x,y(nfxnyn)fnfxn1yn1fn1fxyx)) 得到（因xxn)，所以是外插f(x,y(x)) xxn1 xnxnn

n

n1

n

1[(xx

此式区[

n

xn1f(x,y(x))dxx

xxnff n

nyn

h(3

fn1

注意到插 得出h3，故（21）的局部截断误差为O(h3)，精度比向前Euler提高1阶。若取r2,3,L可以用类似的方法推导 ，如对于r3有yn

h(55

59

n

37

n

9

n

其局部截断误差为O(h5如果将上面代替被积函数f(x,y(x))用的插值由外插改为内插，可进一步小误差。内插法用的是yn1,yn,L,ynr1，取r1时得到的是梯形 ，取r3时yn

h(9 n

19

5

n

fn

与（22）式相比，虽然其局部截断误差仍为O(h5)，但因它的各项系数（绝对值）大为减小，误差还是小了。当然，（23）fn1Euler一y'ifi(x,y1,y2,L,ym (i1,2,L, y(a) iyy1

,L,ym)T，

,L,ym0)T，

y'f(x,y(a)

ax yf(x,y1)f(x,y2)Ly1那么问题（25）在[abyyx全部用于求解问题（25。设有m阶常微分方程初值问题y(m)f(x,y,y',L,y(m1) axy(a)y,y'(a)y(1),L,y(m1)(a)y y1yy2y',Lymy(m1)，问题（26）y'1 y1(a)

2 y2(a)y0 0

m

ym

(a)y(m0y'mf(x,1ym y(a)y （27dyAy(x)yRmA为mA的特征值i(i1,2,LmRei0 (i1,2,L,m) smax|Rei|/min|Rei

方法最好是对步长h不作任何限制。§7解7.1数值7.1.1ode113，其中ode45采用四五阶RK方法，是刚性常微分方程的首选方法，ode23采用二三阶RK方法，ode113采用的是多步法，效率一般比ode45高。

0y(0

f(x,)ypynhf(xn,ynyqynhf(xnh,yp 1(ypyqyn1我们自己编写改进的Euler方法函数eulerpro.m如下：ifnargin<5,n=50;endfori=1:n例 y'2y2x22x，(0x0.5)，y(0)解编写函数文件doty.m是用yfxy右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。例2用RK方法求解y'2y2x22x，(0x0.5)，y(0)doty.my(n例

f(t,y,y',L,y(n1y'''3y''y'y y(0) y'(0) y''(0)（i）y1yy2yy3yy1' y1(0)y'

y2(0) y'3yy

y3(0) 2初值问题可以写成YF(t,Y),Y(0)Y0Yy1y2y3functiontspan=[t0tfinal]是求解区间，y0是初值列向量。在命令窗口输[T,Y]=ode45('F',[01],[0;1;-例 求vanderPol方y''(1y2)y'y（i）y1yy2yy'1y'(1y2)y 书写M文件（1）vdp1.m:functiondy=vdp1(t,y); 函数。对于初值y(0)2,y'(0)0，解title('SolutionofvanderPolEquation,mu=1');xlabel('timet');SolutionofvanderPol321solutionsolution time例5 vanderPol方程，1000（刚性）解(i)书写M文件vdp1000.m:functiondy=vdp1000(t,y);title('SolutionofvanderPolEquation,mu=1000');xlabel('timet');数时，用Dn表示，D4y表示对变量y4阶导数。 式中diff_equation为待解的常微分方程，第1种格式将以变量t为自变量进行求解，第2种格式则需定义自变量var。例 x2y(x2y)y' symsxcondition1，condition2，即为微分方程的初边值条件。7试求微分方程y'''y''x，y(1)8,y'(1)7,y''(2) 求解常微分方程组令格式为：dsolve('diff_equ1，diff_equ2，…'，例 f''3gsing'f'cosf2)0,f(3)3,g(5)1的解。解编写程序如下：x1 a1nX'AX，XM，A M ann这里的’表示对t求导数。eAtXAXX(t0X00X(teA(tt0)X0292 X'

X，X(0)

0 2 解symsX'AXf(t)，X(t0)X

X(t)eA(tt0)XteA(ts)f(s)ds010 0