打靶法求边值问题教材

本科毕业论文（设计）
论文（设计）题目：打靶法求边值问题
学院：理学院
专业：数学应用数学
班级：091
学号：0907010228
学生姓名：钟玲声
指导教师：汪萌萌
2013 年21日
打靶法求边值问题
目录
摘要： (1)
引言： (2)
第一章常微分方程初值问题的解法 (3)
1.1常微分方程的离散化______________________________________________ 3
1.2 欧拉(Euler)方法________________________________________________ 4
1.3 改进的Euler方法_______________________________________________ 6
1.4 龙格一库塔(Runge—Kutta)方法_________________________________
1.5 4阶龙格一库塔公式______________________________________________ 9
1.6 线性多步法_____________________________________________________ 9第二章边值问题的数值解法 (11)
2.1 打靶法________________________________________________________ 11
2.2 差分法________________________________________________________ 15第三章Matlab数值解 (166)
3.1 常微分方程的解法_____________________________________________ 156
3.2 打靶法的matlab实现_________________________________________ 23致谢： (27)
主要参考文献 (27)
摘要
常微分方程在很多领域都有非常重要的应用，然而很多常微分方程的解是无法用解析解写出的，因而要借助于数值方法。

本文介绍了常微分方程边值问题的常见解法，例如：欧拉法，龙格一库塔法等。

而对于常微分方程边值问题，常见的解法有打靶法、有限差分法和有限元法等。

在本文中，我们重点介绍了打靶法，并给出了相关算法，然后结合实例编写程序进行了上机实验。

关键词：常微分方程，初值问题，边值问题，打靶法
Abstract
Ordinary differential equations play an important role in different areas. However, most equati ons cannot be expressed an alytically. We n eed to use nu merical methods. In this paper, we discuss the methods of sol ving in itial value problem (IVP) , such as Euler method Runge-Kutta method. For boundary value problem (BVP), shooting method, finite differenee method (FDM) and finite element method(FEM) are presented. We mainly discuss shooting method and give the algorithm. Numerical experime nt is prese nted in the end of the paper.
Keywords：Ordi nary differe ntial equati on s, i nitial value problem, boun dary value problem, the shooti ng method
引言
虽然常微分方程理论发展已经有几百年，但目前仍然在发展中。

特别是最近
三十年，常微分方程迎来了发展的高峰。

常微分方程边值问题是常微分方程理论的重要组成部分，在众多科学技术领域中有着特别广泛的应用。

打靶法是求解常微分方程边值问题的一种数值方法，它的基本思想是将微分方程的边值问题转化为初值问题来求解，它的比较突出的特点是精度很高，程序很简单，实用性很强。

边值问题：对n阶常微分方程
y(= f(x,y,y ,...,y()),
如果能在不同的两点a和b处，唯一地刻画n个附加条件，并且在区间
a乞t汕上求解，则称此为边值问题。

在微分方程中，所谓的边值问题就是我们给定的一个微分方程和一组被我们称之为边界条件的约束条件。

边值问题的解一般情况下是符合特定的约束条件的微分方程的解。

我们在求解这个微分方程时，除了给出方程的本身，往往还需要提供一定的定解条件。

最常见的就是给出初值问题，也就是说给出的定解条件为初始条件；但是也有一些情况，定解条件要求我们考虑所讨论区域的边界，比如说在一个给定区间讨论时，把定解条件在区间的两个端点给出，给定的这种定解条件就被我们称之为边界条件，与之相应的定解问题我们就称之为边值问题。

第一章常微分方程组初值问题的解法
1.1 常微分方程的离散化
下面主要讨论一阶常微分方程的初值问题，它的一般形式是
¥ = f(x, y) a <x<b
dx ⑴
.y(a)二y o
在下面的讨论中，总假定函数f (x, y)连续，且满足Lipschitz条件，也就是存在常数L ,使得
I f(x,y) -f(x,y)匡L|y-y|
那么，根据常微分方程理论知，初值问题(1)的解存在并且唯一.
所谓数值解法，就是求问题(1)的解y(x)在若干点
a =x0：：：% ：：：X2 ：：：…：::X N=b
处的近似值y n(n=1,2, ,N)的方法，%(n=1,2, ,N)称为问题(1)的数值解,
h n二X n 1 -X n称为由X.到焉1的步长•今后如无特别说明，我们总取步长为常量h.
建立数值解法，首先要将微分方程离散化，一般采用以下几种方法：
1.1.1用差商近似导数
如果用向前差商y(Xn A y(X n)代替y'(x n)代入(1)中的微分方程，则得h
y(X n 十)—y(X n)上/ / 、、/ 01…、
f (X n，y(X n)) (n -0,1,)
h
化简得
y(X n 1) ： y(X n) hf (X n, y(X n))
如果用y(X n)的近似值y n代入上式右端，所得结果作为y(X n 1)的近似值，记为
y n 1，那么有
y n 1 *n hf(X n,y n) (n =0,1, ) (2)这样，问题(1)的近似解可通过求解下面的问题
；%卑=y n+hf (X n,y n) (n =0,1，…)(3)
yo= y(a)
得到，按式(3)由初值y o可逐次算出力」2,….式子(3)是个离散化的问题，称为差分方程初值冋题.
需要说明的是，用不一样的差商近似导数，将得到不一样的计算公式
1.1.2 用数值积分方法
将问题(1)的解表成积分形式，用数值积分方法离散化.例如，对微分方程两端积分，得到
x n -f-
y(X n 1) -y(X n)二f (x, y(x))dx (n =0,1, )
X n
⑷
右边的积分用矩形公式或梯形公式计算
1.1.3 Taylor 多项式近似
将函数y(X)在X n处展开，取一次Taylor多项式近似，则得
y(X n 1)： y(X n) hy'(X n)二y(X n) hf(X n, y(X n))
再将y(X n)的近似值y n代入上式右端，所得结果作为y(X n 1)的近似值y n .1，得到
离散化的计算公式
y n 1 二y n hf (X n,y n)
上面的三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式.其中的Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差.
1.2欧拉(Euler)方法
贵州大学本科毕业论文(设计)
第5页
1.2.1 Euler 方法
Euler 方法就是用差分方程初值问题(3)的解来近似微分方程初值问题(1)
的解，即由公式(3)依次算出y(x n )的近似值y n (n =1,2,…)。

这组公式求问题(1) 的数值解称为向前Euler 公式. 如果在微分方程离散化时，用向后差商代替导数，也就是 y'(x n
.1)坐』炷，则得计算公式
h
Tn* = y n +hf(X n 卅 y n*)(门=0,1，…) 丿
(5)
』0 =y(a) 用这组公式求问题(1)的数值解称为向后Euler 公式.
向后Euler 法与Euler 法形式上相似，但实际计算时却复杂得多.向前Euler
公式是显式的，可直接求解.向后Euler 公式的右端含有ym ，因此是隐式公式，
般要用迭代法求解，迭代公式通常为
1.2.2 Euler 方法的误差估计
对于向前Euler 公式(3)我们看到，因为n =1,2,…时公式右端的y 都是近 似的，所以用它计算的y n 1会有累积误差，分析累积误差比较复杂，这里先讨论 比较简单的所谓局部截断误差•
假设用(3)式时右端的y n 没有误差，即y n =y(x n )，那么由此算出
y n 1 =y(X n ) hfb n ’yd n )) (7)
局部截断误差指的是，按(7)式计算由X n 到X n1这一步的计算值y n ・1与精确值 y(X nd 之差y(X n1)-y n1.为了估计它，根据Taylor 展开得到的精确值y(X n J 是 h 2 3
y(X n 1)= y(X n ) hy'(X n ) y''(X n ) O(h ) (8)
2 (7)、(8)两式相减(注意到y'=f(x,y))得到
y n* =yn +hf (X n ,y n )
(k 卅)
J n ^ =yn hf (x n 1, y n k)1) (k =0,1,2,)
(6)
也就是局部截断误差是h 2阶的，而数值算法的精度定义为:
如果一种算法的局部截断误差为 0(h p 1)，则称该算法具有p 阶精度.
显然p 越大，方法的精度越高。

式(9)说明，向前Euler 方法是一阶方法, 因此它的精度不高。

1.3 改进的Euler 方法
1.3.1梯形公式
利用数值积分方法将微分方程离散化的时候，如果用梯形公式计算式 (4)中 之右端积分，即
x n 1 h
X f (x,y(x))dx 拓二[f (X n , y(X n )) + f (X n*, y(X nG )]
X n
并用Y n .Y n 1代替y(X n ),y(X n 1)，那么得计算公式
h
y n ^y n ■-[ f (X n ,y n ) f (X n 1, y n 1)]
2
这就是求解初值问题(1)的梯形公式.
直观上容易看出，用梯形公式计算数值积分要比矩形公式好.梯形公式为二 阶方法。

梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为
匕¥ =y n +hf (X n ,y n )
“ y n 「= y n +£[ f (X n , y n ) + f (X n*, 丫补站]
(10)
2
Jk =0,1,2,…) 由于函数f (x, y)关于y 满足Lipschitz 条件，容易看出
I (k 1) (k) |.... hL (k) (2) 1
| y n 1 - y n 1 | 2 | y n 1 一 y n 1 | 其中L 为Lipschitz 常数.因此，当0 ：：：匹::1时，迭代收敛.但是这样做计算量
2
较大.如果实际计算时精度要求不太高，用公式(10)求解时，每步可以只迭代 一次，由此导出一种新的方法一改进 Euler 法.
y(X n 1
h 2 )—y n i — 2 y ''(X n ) O (h 3) ： O (h 2) (9)
1.3.2 改进Euler 法
按照式(5)计算问题(1)的数值解时，如果每步只迭代一次，相当于将
Euler公式与梯形公式结合使用：先用Euler公式求1的一个初步近似值y n 1，
称为预测值，然后用梯形公式校正求得近似值y n 1，即
y n ^y n hf(X n,y n) 预测
{ h ( 11)
Y n 1 = Y n (X n』n) f(X n n 1 )]
式(11)称为由Euler公式和梯形公式得到的预测一校正系统，也叫改进Euler 法.
为了便于编制程序上机，式子(11)常改写成
L
y p = y n +hf(X n,y n)
* y q = y n +hf (X n +h,y p) (12)
1
% 卅=2(y p +y q)
改进Euler法是二阶方法.
1.4龙格一库塔(Runge-Kutta )方法
回到Euler方法的基本思想一用差商代替导数.实际上，按照微分中值定理
应有
h
注意到方程f (x, y)就有
y(X n 1) =y(X n) hf (X n F,y(X n 巾)) (13)
不妨记K =f(X n •击,y(X n Th))，称为区间[X n,X n1〕上的平均斜率.可见给出一
种斜率K，(13)式就对应地导出一种算法.
向前Euler公式简单地取f(X n，y n)为K，精度自然很低.改进的Euler公式
可理解为K取f(X n，y n) , f (X n J 1)的平均值，其中?n ^n - hf(X n，Y n)，这
种处理提高了精度•
如上分析启示我们，在区间［X n，X n』内多取几个点，将它们的斜率加权平均
作为K，就有可能构造出精度更咼的计算公式.这就是龙格一库塔方法的基本思想.
首先不妨在区间［X n,X n 1］内仍取2个点，仿照(13)式用以下形式试一下
y n 1 =yn h( i k i —2k2)
k i = f (X n, y n) ( 14)
k2 二f (X n —h y n ：hk i), 0 ：：:，: ：：1
其中1, -2^'，1为待定系数，看看如何确定它们使(14)式的精度尽量高.为此我
们分析局部截断误差y(X n 1)- y n・1，因为y n =y(X n)，所以(14)可以化为
y n+ = y(X n) +h(X j k1 + 丸2k2)
k1 = f (X n, y(X n)) = y'(X n)
^2= f (X n +ah,y(X n) + Bhk1)(15)
=f(X n, y(X n)) +Ghf x(X n, y(X n))
[+0hk1f y(X n,y(X n)) +O(h2)
其中k2在点(X n, y(X n))作了Taylor展开.(15)式又可表为
2 0 3
Y n 1 =y(X n) ( ‘1 ’2)hy'(X n) ‘2： h ( f x ff y) O(h )
a
注意到
h23
y(X n1)= y(X n) hy'(X n) y''(X n) O(h )
2
中f，y'^ f x ff y，可见为使误差y(X n .1)- y n .1 二O(h3)，只须令
r , 1 P
'「‘2=1，‘ —，- 1 ( 16)
2 :
待定系数满足(16)的(15)式称为2阶龙格一库塔公式.由于(16)式有4个
未知数而只有3个方程，所以解不唯一.不难发现，若令= 2 =—，〉=：=1，
2
即为改进的Euler公式.可以证明，在[X n，X n』内只取2点的龙格一库塔公式精
度最高为2阶.
1.5 四阶龙格一库塔公式
要进一步提高精度，必须取更多的点，如取 4点构造如下形式的公式：
= yn
+h 仏 & +扎2k 2 "冰3 +入4k 4)
k i = f (X n , y n )
*2 =
f(X n +%h,y n +Rhk i )
( 17)
k 3 = f (X n +吋,y n + %hk i +^3hk 2) k 4 二 f(X n ： 3h,y n Lhk i Jk ?
：6hk 3
)
其中待定系数i 〉i 「i 共13个，经过与推导2阶龙格一库塔公式类似、但更复杂 的计算，得到使局部误差y(X n .1)- y n .1 =0(h 5)的11个方程•取既满足这些方程、 又较简单的一组，可得
y n 1 =6你1 2k 2 2k 3 k 4) k^ f (X n ,y n )
h hk 1
弦2 = f (X^-, y n + 二-)
2
2
k 3 = f (X n 弓 y n +导)
2 2 k^ f (X n h, y n hk 3)
这就是常用的4阶龙格一库塔方法(简称RK 方法) 1.6线性多步法
以上所介绍的各种数值解法都是单步法，这是因为它们在计算 y n.1时，都只
用到前一步的值y n ，单步法的一般形式是
Y n
Y n h (X n .Y n .h) (n =0,1, ,N -1)
其中:(x, y, h)称为增量函数，例如Euler 方法的增量函数为f (x,y)，改进Euler 法的增量函数为
1
(X, y,h) [ f (x, y) f (x h, y hf (x, y))]
2
怎么样通过较多地利用前面的已知信息，如 y n ,y n v,…，y n _r ，来构造高精度
(18)
(19)
的算法计算y n.1，这就是多步法的基本思想.经常使用的是线性多步法.
从用数值积分方法离散化方程的（4）式
X
n +
y （X n 1） —y （X n ） = . f （x, y （x ））dx
x
n
出发，记 f (X n , Y n )二 f n ， f (X n J, Y n 」)二 f （X n 」，f n 」），（X n ,f n ）的插值公式得到（因X _ X .），所以是外插.
此式在区间［X i ，X n 』上积分可得
X
n 1
3h x f (x, y(x))dx f
X n
2
于是得到
注意到插值公式（20）的误差项含因子（x -x n 」）（x - X n ），在区间［x n ,x n 上积分 后得出h 3，故公式（21）的局部截断误差为O （h 3），精度比向前Euler 公式提高
若取r =2,3;可以用类似的方法推导公式，如对于
—3有
h
Y n 1 = Y n 丁（55仁-59f n 」37 f
24
其局部截断误差为O （h 5）.
如果将上面代替被积函数f （x,y （x ））用的插值公式由外插改为内插，可进
r =3时可以得到
Y n 1 二 Y n 丁旧仁 1
19
齢-5f
24
与（22）式相比，虽然其局部截断误差仍为O （h 5），但因它的各项系数（绝对值） 大为减小，误差还是小了 •当然，(23)式右端的花.1未知，需要如同向后Euler 公式一样，用迭代或校正的办法处理•
让我们试验一下r =1，即利用 Y n , Y n J 计算Y n 1的情况.
式中被积函数f （x,y （x ））用二节点
X -X n 」
f (X, Y (X )) = f n
——f X n - X 」
_ 1 _ h
[(X -X nj ) f n X — X n
n J
X nJ -X n
-(X -X n )f n 」] (20)
(21)
(22)
步减小误差.内插法用的是
Y n 1, Y n ,
，Yn..1，取r = 1时得到的是梯形公式，取
仁）
(23)
贵州大学本科毕业论文(设计) 第11页
第二章边值问题的数值解法
2.1打靶法
解两点边值问题
y = f (x,y, y ), a z x^b,—：： :: y ：：：：：, ( 1)
y(a) =；;, y(b)=:
的打靶法实质上是把边值问题化为初值问题来解，我们设法确定y (a)的值t，使
得初值问题
y = f(x,y,y),a Ex^b,i： ::: y < ::, (2)
y(a) = : , y(b)=:
的解y(x,t)在X =b的值y(b,t)满足
y(b,t)二-
或者是
y(b,t)-—；
其中；为我们允许的误差界限.这样，我们把y(b,t)看作边值问题(1)的近似解. 为此，可以采用逐次逼近的方法来实现.
假设y(x)为边值问题的解，我们估计y (a)的值为t。

之后，求解初值问题
y = f(x,y, y ), a :: y ：：二，
y(a)二a, y (a)二t。

(在(2)中令t =1。

).这样得到的解为y(x,t。

)，并且计算得到y b to)^#。

. 一般地，
P。

式3 .但是如果P。

= B或0。

-円v乞，则把y(x,t0)作为边值问题(1)的近td
似解；否则的话，贝M故适当的调整to，例如取11。

，则在(2)中令ti，
贵州大学本科毕业论文（设计）第12页再求解此初值问题•设计算得到它的解为y（x,tj , y（b,tj = X，若-或
上面的积分曲线（y = y（x,t。

））的咕过大，下面的曲线（y = y（x,tj）的打过小.
我们已经介绍了解常微分方程初值问题的数值方法.现在的问题是怎么样确定参
数t k.当t二t k时，从初值问题计算得到解y（x,tQ.自然地，我们希望
k imybtQ = y（b）= E .
因此，确定t k的问题可以归纳为求方程
y（b,t）「J o （3）
的近似根•（3）是一个非线性方程.二分法，割线法和牛顿法都可以用来求解方程（3）.
假如用弦割法来求解方程（3），我们需要选取初始近似t o，t i，由公式
t t(yQt kj -0)仇」-t k/) k 23...
I k 二I k 丄,k =厶3,
- y(b,t kj - y(b,t k/)
贵州大学本科毕业论文(设计) 第13页
(4)
生成序列｛t｝按照(4)的方式求t k，直到y(b,t k)- 1「为止，其中■:为允许的误差界.
令h=(b_a)/N必=a ih(i =1, ,N _1,N).给定初值t。

，误差容限TOL，最大迭代次数m.
取t二t o，解初值问题(2)得解y(x,t o)在x = x(i =1,…，N)的近似值：
y(X1,t°), ,y(X N」,t o),y(X N,t o)(二y(b,t°)).
如果|y(X N,t。

)-円＜TOL，则输出
y(X1,t°), ,y(X N」,t o),y(X N,t°)
作为初值问题(1)的解在％，…,X M d, X M的近似值；停机
B
令t =t1t0，解初值冋题(2)得解y(X,t1)在x = X j(i=1,…,N)的近似值: y(X N,t o)
y(X1 ,tj, , y(X N “tj, y(X N ,t1)(= y(b, tj)
如果y(X N,tj-円＜TOL，则输出
yd」),』区小)，丫区我1)
作为初值问题(1)的解在％，…,X N 4, X N的近似值；停机.
对k =2,3,…,m重复做。

由(4)公式计算t k;
令t =t k，解初值问题(2)的解y(X,t k)在x = X j(i =1/ , N)的(近似)值:
y(x1,t k), ,y(x N 4,t k),y(X N,t k)(二y(b,t k))
如果y(x ,tj -円＜TOL，则输出
丫(也),,y(X N」k),y(X N,t k)
作为初值问题(1)的解在X1,…,X N4,X N的近似解；停机.
贵州大学本科毕业论文(设计)
第14页
对第三边值问题，我们同样可以用打靶法求解•设第三边值问题为
/ = f(x,y,y),
P i y(a) qy (a)「，口| 书=0, (6) P 2y(b)+q 2y"(b) =
P 2 +|q 2 式0・
(7)
打靶法的基本过程如下：
选取参数t 0，令y (a)二t 0.设p<■- 0，则可以由(6)式确定y(a)，从而得到 初值问题
”y"=f (x,y, y),a 兰xEb,
1
y(a) (:-q i t o
), ( 8)
P i
y (a) =t°.
求初值问题(8)的解y(x,t o )代入(7)式左端得
二 P 2y(b,t o ) q 2y(b,t o ).
如果飞「，则y(x,t 。

)为所求解的边值问题的近似解；如果 ,=1 ,贝冷
设y(a)=t i (以t i 取代t 。

)解初值问题(8)，求得解y(x,t i )后在代入(7) 式左端得
：i =P 2y(b,t i ) q 2y (b,t i ).
从t 0,t i 「0「i 出发，由割线法
产生t 2，以t 2取代t 0解决初值问题(8)确定匕，如此继续进行下去，求得序列｛t k ｝，
直到
P 2y(b,t k ) +q 2y"(b,t k ) — 円 <TOL
为止，TOL 为误差容限.
如果P i=P 2=0，则(5)，(6)，( 7)为第二边值问题，处理方法完全类似.
(5) t k 二t k4
(：
k4 - -)(tk^ -tk^)
Wk 」-h)
k = 2,3「. (9)
贵州大学本科毕业论文（设计）第15页2.2差分法
在这里我们介绍一种求解（1）的简单方法，我们以最简单的二阶常微分方程为例来说明边值问题的数值解法•
解法可以用数值微分公式替代导数，将其变成代数方程，然后我们再求解，这种求解的方法我们通常把它叫做差分法•
y （X
k）（h2）, （即中心差商）
2h
y（Xkr y（X2）—2y（X k）n .°（h2）. （即二阶差商）
h
令
y（X k J y（xj =y k,y（X k 1）=y k 1,
并略去y（xQ和y”（xj用差商表示中h的二次项，然后代入（1），这样我们就得到：
y k卅一2y k +y k4 =h2f （X k,y k, ** %」），
2h
* k =1,2，…，n -1, （10）y（a）" =y°,
y（b）「二y n.
b _ a
其中h =，兀=x0 kh（k =1,2, ,n），Xo=a.
n
从（10）中解出y1,y2 / ,y n，即为所求的常微分方程的计算解.一般情况下，（10）是一个非线性方程组，求解也比较困难.但是，如果微分方程（1）中第一式的右端f为y, y的线性函数时，则方程组（10）就可以转化成线性方程组，求解问题
也就简单多了.
第三章Matlab 数值解
3.1常微分方程的解法
3.1.1利用Matlab求解一阶方程初值问题
在Matlab的工具箱中，给我们提供了很多求解常微分方程的功能函数.如ode45, ode23, ode113,其中ode45采用四五阶RK方法，是解非刚性常微分方程的首选方法，ode23采用二三阶RK方法，ode113采用的是多步法，效率一般比ode45高.
Matlab的工具箱中没有Euler方法的功能函数.
对简单的一阶方程的初值问题
(y'=f(x, y)
MX。

)= y。

改进的Euler公式为
y p =y n +hf (X n,y n)
* y q = y n +hf (X n +h,y p)
1
y n 卅「(y p +y q)
一阶微分方程初值问题的求解在Matlab中非常简单明了并且容易求解的方法，而且由于Matlab中内嵌的微分方程初值问题的求解方法是变步长的龙格库塔法，故而求解的结果是值得我们相信的.
= f (x, y)
八"(a兰x兰b)
对于微分方程j y(a)=y0，直接利用Matlab语言的内嵌函数ode45,按照调用格式[x,y]=ode45( ‘ f' , x, y0)即可求出最终的结果.其中
指函数名或所建立的M文件的文件名，x为已知离散节点向量，y。

为初值.
例：求解一阶微分方程初值问题：
2x
y (0乞x乞1)
M0) =1 (1)
先编写M文件fun.m:
Function f=fun (x,y)
f=y-2*x/y;
命令窗口写入:
x=0:0.2:1;
[x,y]=ode45( ‘ fun ' ,x,1)
求解结果为：
[x,y]=0.0000 1.0000
0.2000 1.1832 0.4000 1.3416
0.6000 1.4832
0.8000 1.6125
1.0000 1.7321
这个结果与四阶龙格库塔法求解结果相同，比改进欧拉法求解结果要精确得多.运用内嵌函数直接求解最大的优点是简洁方便，只要依照函数的调用格式调用即可，但是对于理解函数值解法的思想没有太大的帮助，故而要体现数值解法的具体思想，可以通过编写程序来实现.
3.1.2根据数值解法的思想编程实现
对于一阶微分方程的初值问题
y " = -xy2
y (0兰X
兰5)
畀(0)=2
(2)
运用改进欧拉法和四阶龙格库塔法求解程序如下:
(1)改进欧拉法
建立如下M文件：
a=0;
b=5;
n=120;
xn=0;
yn=2;
h=(b-a)/n;
for i=1: n
贵州大学本科毕业论文（设计）第18页yp=yn-h*x n*yn*yn;
xn=xn+h;
yc=yn-h*(x n*yp*yp);
yn=(yp+yc)/2;
end
迭代次数n=120时求解结果：yn=0.07694733409450;若n=5000时,求解结果yn=0.07692309052101.
标准四阶龙格库塔法
对于上述初值问题建立如下M文件：
a=0;
b=5;
n=20;
xn=0;
yn=2;
h=(b-a)/h;
for i=1: n ；
k1=-h*x n*yn*yn;
xn=xn+h/2;
k2=-h*x n*(y n+k1/2)*(y n+k1/2);
k3=-h*x n*(y n+k2/2)*(y n+k2/2);
xn=xn+h/2;
k4=-h*x n*(y n+k3/2)*(y n+k3);
yn=y n+(k1+2*k2+2*k3+k4)/6;
end
当迭代次数n=120的时候，获得计算结果：yn=0.07692668513821，比运用改进欧拉法迭代120步的效果都要好的多；儿当迭代次数n=120时，迭代结果
yn=0.07692307941295基本已经达到精确解.
3.1.3结果分析
对问题(1)的求解简单方便，但难于体现求解的具体思想：对于问题(2)
的求解运用变成实现，从实验角度帮助理解理论知识，不仅可以深入掌握求解方法的思想，而且还可以提高实践动手能力.
问题(2)中，初值问题的精确解为yn=0.07692307…，改进欧拉法在迭代到
120步的时候误差还较大，迭代到5000步的时候误差较小，但是与精确解的误差
仍然有出入，而且从迭代过程来看，迭代次数n越大，收敛速度越慢，效果
不明显.而四阶龙格库塔法迭代20步的结果远远好于改进欧拉法120步的结果, 当迭
代到120步时已经跟精确解非常接近.故而龙格库塔法迭代不仅得到的解比
改进欧拉法准确得多，而且收敛速度非常快，确实是一种很好的算法.但是唯
不足的是计算量大，这一矛盾很难解决.
3.2 打靶法的Matlab实现
3.2.1打靶法算法
我们考虑两个边界条件的二阶常微分方程
y 二 f (x, y, y ), a ：： x ：：： b
«y(a)=a ，
L y(b) = P
其中a,b「, ■-是给定的常数，y是关于x的未知的函数，f是一个函数的微
分方程的指定者，实现了y与x之间的关系.这是一个两点之间的边值问题.初始
值问题需要两个条件被给定在相同的x值下.例如，y(a)二:和y(a)二.因为两
个独立的边界条件，上述两点边值问题是比较难解决的问题.
“打靶法”是一个初始值问题代替上述问题的基本思想.当然，我们并不知
道在x二a时y的导数值，但是我们可以猜测，然后进一步完善猜测的迭代.更确
切的说是，我们对y(a)为未知，并使用割线法或者牛顿法(或者其他的方法求
解非线性方程组)来确定y (a).
我们引入一个函数u，这是x的函数，但是它取决于一个参数t.即u=u(x,t).
我们使用「和u ”表示u相对于x的偏导数.我们希望y是准确的，如果t是正确的
选择.U 通过以下定义：
u — f (x,u,u )
* u(a,t) =a
u "(a,t) =t.
当给定t 的值，然后就可以求解上述关于u 的初始值问题.一般u 与y 是不等 的，由于u(a)=t=y(a).但是如果t 等于y (a)，则我们认为u = y.因为我们不知 道y(a)，我们确定它在x 二b 边界条件.即,我们要求解t 使得：
(t)二 u(b,t)— 1 = 0.
如果我们找到了 t 使得(t^0，这意味着u(b,t)=2.因此，在x 二a 和x 二b ， u 满足两个相同的边界条件，换句话说，
u 二y.因此，()=0的解t 必须满足
t 二 y(a). 如果我们可以写出u 初始值问题(任意t )的解析式，我们可以写一个公式
(t)二u(b,t) -二当然，这一般是不可能的.然而，没有一个解析公式，我们仍然 可以求解(t) =0的数值.对任何t ，初始值问题u 的数值方法可以用来找到一个 近似解u(b,t)(因此(t)).最简单的方法是使用割线法.
t j 一― t j 円一(tj-g)(t j )，jg ，3,
为了达到这个目的，我们需要两个初始猜测：t 0和1.我们也可以使用牛顿
法:
我们需要一种方法来计算导数 (t)，由于"t)二u(b,t)- 1 ,我们有
专(b,t). :t
如果我们定义了 V (X,t)二-：＜ ；：
t ，我们有一下的初始值问题V :
f u (x,u,u)v+ f u ・(x,u,『)v
«v(a,t) =0
V (a,t) =1. t j1 t j
(t j ) 0,1,2；
(t)二
在这里，对于x , v和v“是V的第一和第二阶偏导数•方程的上面设置是从
偏导数相对于u系统x获得.链式法则是用来获得V的微分方程•现在，我们设(t) =v(b,t).这是打靶法的方法，其中包括牛顿法求解的算法(tH0 : t o =初始猜测值y (a).
for j=0,1,2,...
求解如下的微分方程系统，求解区间从x=a到x=b
u =f(x,u, u )
U
x』=t j
v= f u (x,U,U)v f u(X,U,U)V
Vx£ =0
V x」=1.
V|x出
如果我们想使用在前面解决了的上述系统中的两个第二阶方程U和V的方法,我们需要引入一个向量z=(u,u,v,v)T和微分方程为z'F(x,z)的向量F.初始条件是z(a) =(：,t j,0,1)T.
打靶法也使用于特征值问题:
y"二f(x,y, y ：■), a ：：x :: b, y(a) =0, y(b) =0.
在f满足条件f(x,0,0, ')=0或者跟普遍，f在y是均匀的，即,
f (x, cy,cy , ■) =cf (x, y, y ,').
贵州大学本科毕业论文(设计) 第22页对于任何常数c.注意y = 0总是一个上述边值问题的解.事实上，一个特征值问题是一个特殊的边界值问题满足y=0始终是一个解决方案，有一个称为■的方程参数(或边界条件).特征值问题是确定非零解为•的特殊价值而存在的. 特征值问题的解决方案是对・，y(x)1,其中•是y特征值问题的特征.通常，有许多(可能是无限的)的特征值和本征函数.利用打靶法，我们考虑初始值问题
u = f (x,u,u , ■), x a,
u(a, •) =0,
u(a, ■) =1,
其中■为参数.由于该解决方案依赖于■，我们用符号u二u(x, •)，但u和u 代表第一和第二阶导数取决于x.对于任何给定的■，我们可以解决上述的初始值问题.现在假设■满足条件
(，)=u(b,，) = 0,
则y(x) =u(x, ■)是我们正在寻找的特征函数，•是相应的特征值.因此，我们只需要使用割线法或者牛顿法求解■的方程(y^0.如果使用割线法，我们只需要解决不同的迭代器的■初值问题.如果用牛顿法，我们必须得到关于y的
(y).因此，我们需要
v(x, ■ H —(x,').
C/L.
我们需要一个初始值问题的V.这可以通过偏导数相对于■的出数值问题，得到u.我们有
V = f u(x,u,u , )v f u (x,u,u , ■ )v f .(x,u,u , ■ ),x a,
v(a, )=0,
v(a, ■) = 0.
请注意，我们已经使用链式法则(微积分)为v得到方程.现在，你可以解
决初始值问题(和u —起)，然后评估 (■)为任何给定的•.
3.2.2实例
(1) 用打靶法算法编写程序上机求解下面的实例：
y" = (1+x2)y,—1 兰x 兰
1, *1—1.
解编写算法如下(我们用D2i表示u的2阶导数，用Du表示u的一阶导数；用D2v
表示v的2阶导数，用Dv表示v的一阶导数)：
to=o.
for j=0,1,2，…
solve the following system numerically from x=-1 to x=1
D2u=f(x,u,Du)
u(x=-1)=1
Du(x=-1)=t(j)
D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv
v(x=-1)=0
Dv(x=-1)=1.
set
t(j+1)=t(j)-(u(x=1)-1)/v(x=1).
plot(t,y)
用上面的程序上机进行测试，就会得到一个很接近精确解的一个解
(2) 用打靶法算法编写程序上机求解下面的实例:
3T
y 4y = cosx,0 咗x ,
4
y(0) =0, y(-H0.
解编写算法如下(我们用D2i表示u的2阶导数，用Du表示u的一阶导数；用D2v 表示v的2阶导数，用Dv表示v的一阶导数)：
t0=0.
for j=0,1,2，…
solve the followi ng system nu merically from x=0 to x= n /4
D2u=f(x,u,Du)
u(x=O)=O
Du(x=0)=t(j)
D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=O)=O
Dv(x=0)=1.
set
n /4)-0)/v(x= n /4).
t(j+1)=t(j)-(u(x=
Plot(t,y)
(3) 用打靶法算法编写程序上机求解下面的实例：
jy” + (x+1)yJ2y = (1_x2)e」,0兰x 兰1, y(0)=y(1)=0.
解编写算法如下(我们用D2i表示u的2阶导数，用Di表示u的一阶导数；用D2v
表示v的2阶导数，用Dv表示v的一阶导数)：
t0=0.
for j=0,1,2，…
solve the followi ng system nu merically from x=0 to x=1
D2u=f(x,u,Du)
u(x=0)=0
Du(x=0)=t(j)
D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv
v(x=0)=0
Dv(x=0)=1.
set
t(j+1)=t(j)-(u(x=1)-0)/v(x=1).
Plot(t,y)
(4) 用打靶法算法编写程序上机求解下面的实例:
y —y=x,x [0,1],
y(0) =0,
y(1)
解编写算法如下(我们用D2i表示u的2阶导数，用Du表示u的一阶导数；用D2v 表示v的2阶导数，用Dv表示v的一阶导数)：
t0=0.
for j=0,1,2,...
solve the followi ng system nu merically from x=0 to x=1 D2u=f(x,u,Du)
u(x=0)=0
Du(x=0)=t(j)
D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=0)=0
Dv(x=0)=1.
set
t(j+1)=t(j)-(u(x=1)-1)/v(x=1).
plot(t,y)
(5) 用打靶法算法编写程序上机求解下面的实例：
'(1+x2)y“_xy，_3y =6x_3,x^[0,1],
“ y(0) =1,
J(1) =2.
解编写算法如下(我们用D2i表示u的2阶导数，用Du表示u的一阶导数；用D2v 表示v的2阶导数，用Dv表示v的一阶导数)：
t0=0.
for j=0,1,2,...
solve the followi ng system nu merically from x=0 to x=1 D2u=f(x,u,Du)
u(x=0)=1
Du(x=0)=t(j)
D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=1)=0
Dv(x=1)=1.
set
t(j+1)=t(j)-(u(x=1)-2)/v(x=1).
Plot(t,y)
(6) 用打靶法算法编写程序上机求解下面的实例：
y -2y = x,x [1,2],
y(1) =0,
y(2)=2.
解编写算法如下(我们用D2i表示u的2阶导数，用Di表示u的一阶导数；用D2v 表示v的2阶导数，用Dv表示v的一阶导数)：
t0=0.
for j=0,1,2,...
solve the followi ng system nu merically from x=1 to x=2 D2u=f(x,u,Du)
u(x=1)=0
Du(x=1)=t(j)
D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=1)=0
Dv(x=1)=1.
set
t(j+1)=t(j)-(u(x=2)-2)/v(x=2).
Plot(t,y)
(7) 用打靶法算法编写程序上机求解下面的实例:
“ 4 "丄2 2ln x
y =——y +py ——— x x x
1
y(1) ,y(2) Rn2.
L 2
to=o.
for j=0,1,2,...
solve the followi ng system nu merically from x=1 to x=2
D2u=f(x,u,Du)
u(x=1)=-0.5
Du(x=1)=t(j)
D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv
v(x=1)=0
Dv(x=1)=1.
贵州大学本科毕业论文(设计) 第28页set
t(j+1)=t(j)-(u(x=2)-l n2)/v(x=2).
Plot(t,y)
致谢：
在本论文完成之际，我要由衷感谢汪萌萌老师在课题设计和论文写作上的精
心指导，同时对所有帮助过我的老师以及同学致以最真诚的感谢。

主要参考文献：
[1] 施吉林，刘淑珍，陈桂芝•计算机数值方法(M).北京.高等教育出版
社.2009.165.
[2] 安乐.打靶法在常微分方程边值问题中的一些应用[J].科技广
场.2011年.05期.
[3] 夏必腊；王金山.自共轭常微分方程边值问题与变分问题的等价性
[J].大学数学.2011(03)
[4] 王国栋.常微分方程中数学建模思想的教学探讨[J].中国科教创新
导刊.2011年.22期.
[5] 王国栋.常微分方程中数学思想方法教学探讨[J].科教新报(教育科
研).2011 年.22 期.
[6] 朱美玲.常微分方程在数学建模中的应用[J].太原城市职业技术学
贵州大学本科毕业论文（设计）第28页院学报.2011.05期.
[7] 吴德林.常微分方程教学探索[J].读与写（教育教学刊）.2011年04 期.
[8] 张建龙.常微分方程教学改革的探索[J].才智.2011.16期.
[9] 李培峦.几类非线性微分方程边值问题解的存在性及多解性研究[D].中南大学.2010年.
[10] 连海强.几类非线性微分方程边值问题的解及应用[D].曲阜师范大学.2010年.。