微分方程的边值问题

微分方程边值问题的数值方法

本部分内容只介绍二阶常微分方程两点边值问题的的打靶法和差分法。 二阶常微分方程为

(,,),y f x y y a x b '''=≤≤

(1.1)

当(,,)f x y y '关于,y y '为线性时，即(,,)()()()f x y y p x y q x y r x ''=++，此时(1.1)变成线性微分方程

()()(),y p x y q x y r x a x b '''--=≤≤

(1.2)

对于方程(1.1)或(1.2)，其边界条件有以下3类： 第一类边界条件为

(),()y a y b αβ==

(1.3)

当0α=或者0β=时称为齐次的，否则称为非齐次的。 第二类边界条件为

(),()y a y b αβ''==

(1.4)

当0α=或者0β=时称为齐次的，否则称为非齐次的。 第三类边界条件为

0101()(),()()y a y a y b y b ααββ''-=+=

(1.5)

其中00000,0,0αβαβ≥≥+>，当10α=或者10β=称为齐次的，否则称为非齐次的。微分方程(1.1)或者(1.2)附加上第一类，第二类，第三类边界条件，分别称为第一，第二，第三边值问题。

1 打靶法介绍

下面以非线性方程的第一类边值问题(1.1)、(1.3)为例讨论打靶法，其基本原理是将边值问题转化为相应的初值问题求解。

【原理】假定()y a t '=，这里t 为解()y x 在x a =处的斜率，于是初值问题为

(,,)

()()y f x y y y a y a t α

'''=⎧⎪

=⎨⎪'=⎩

(1.6)

令z y '=，上述二阶方程转化为一阶方程组

(,,)()()y z

z f x y z y a z a t

α'=⎧⎪'=⎪

⎨

=⎪⎪=⎩ (1.7)

原问题转化为求合适的t ，使上述初值问题的解(,)y x t 在x b =的值满足右端边界条件

(,)y b t β=

(1.8)

这样初值问题(1.7)的解(,)y x t 就是边值问题(1.1)、(1.3)的解。而对给定的t ，求(1.7)的初值问题可以用欧拉方法、龙格-库塔方法等初值问题的数值解法求解。

理论上(,)y x t 是隐含t 的连续函数，如果(,)y x t 已知，要使得(1.8)成立，可以通过求非线性方程(1.8)的零点来得到合适的t ，这可用任何方程求根的方法，例如牛顿法、或者其它迭代法。

实际上，(,)y x t 是很难找到的，因此必须寻找满意的离散解数值解。

下面叙述打靶法的计算过程：（这里ε为允许误差，t 的修改使用线性插值方法） Step 1：先设0t t =，求解初值问题(1.7)，得到00(,)y b t β=；

若0||ββε-≤，则0(,)(0,1,

,)j y x t j n =为问题(1.7)的满意的离散解，结束；

Step2： 若0||ββε->时，令1t t =，求解初值问题(1.7)，得到11(,)y b t β=； 若1||ββε-≤，则1(,)(0,1,

,)j y x t j n =为问题(1.7)的满意的离散解，结束；

否则转Step3；

Step3：由线性插值得到一般计算公式

1(,)11(,)(,)

(),1,2,

k k k y b t k k k k y b t y b t t t t t k β--+--=-

-

= (1.9)

Step4： 令1k t t +=，求解初值问题(1.7)，得到11(,)k k y b t β++=；

若1||k ββε+-≤，则1

(,)(0,1,)j k y x t j n +=为问题(1.7)的满意的离散解，结束；

否则转Step3。

这个过程好比打靶，k t 为子弹发射率，()y b β=为靶心，当||k ββε-≤时则得到解，故称打靶法。

【例1】用打靶法求解非线性两点边值问题

335

34216(2)8

,23(3)y yy x y x y '''⎧+=+⎪

=≤≤⎨⎪=⎩

要求误差6

0.510ε-≤⨯。精确解为 28

()x

y x x =+

。

【解】：首先将原问题化成初值问题

3

42

4

(2)8(2)yz x k y z

z y z t

'=⎧⎪'=-++⎪⎨

=⎪⎪=⎩

对每个k t ，使用4阶RK 方法求解上述问题，即利用公式

112346(22)h n n y y k k k k +=++++

其中 112122(,),(,)h

n n n n k f t y k f t h y k ==++，

1324322(,),(,)h

n n n n k f t h y k k f t h y hk =++=++

计算，取步长为 h=0.02。

Step1 选择0 1.5t =，求得0(3,)11.4889y t =，3503|(3,)|0.1777y t ε-=>； Step2 选择1 2.5t =，求得1(3,)11.8421y t =，3513|(3,)|0.0755y t ε-=>；

Step3 根据01,t t 以及0(3,)11.4889y t =和1(3,)11.8421y t =，利用公式(1.9)，计算得到

110(3,)35/3

2110(3,)(3,)() 2.0032251

y t y t y t t t t t --=--= Step4 对2t ，利用RK 方法求解，计算得到2(3,)11.6678y t =，3523|(3,)|y t ε->，

转Step3。重复Step3和Step4，可求得

331.999979,(3,)11.66659t y t ==； 442.000000,(3,)11.66666667t y t ==，

满足要求，此时解4(,)(0,1,,)j y x t j m =即为所求。

对于第二类、第三类边值问题也可以作类似处理。例如，对第二类边值问题，它可以转化为以下边值问题