常微分方程初值问题数值解法

常微分方程初值问题数值求解

学生：张玉娟 学号：3070942232 指导老师：梁鹏

摘 要：数值求解是指在没有办法知道未知函数的解析表达式的情况下，近似计算未知函数在其定义域

中的某些离散点上的函数值。本文利用欧拉方法、梯形方法和龙格_库塔方法三种单步法基于Matlab 求解常微分方程初值问题。

关键词：常微分初值问题；数值求解；欧拉方法；梯形方法；龙格_库塔方法；Matlab 实现

1 引言

很多科学技术和工程问题常用微分方程的形式建立数学模型，因此微分方程的求解是很有意义的。建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些典型的方程，而对于绝大多数的微分方程问题，很难或者根本不可能得到它的解析解，实际问题终归结出来的微分方程主要靠数值解法。因此，研究微分方程求解的数值方法是非常有意义的。本文介绍了欧拉法、梯形法和四阶龙格_库塔方法三种单步法，通过Matlab 的平台运行。

2 建模

常微分方程初值问题的数学模型是：求y ()y x =，使之满足

0'()(,)()y x f x y x b y a y =≤≤⎧⎨=⎩，(a ),,

(1)

其中0y ，a ，b 是已知的常数，(,)f x y 是已知函数，且满足条件：

(,)(,')'f x y f x y L y y -≤-，

式中的L 是不依赖于y ，'y 的常数，称为利普希茨（Lipschitz ）常数。

由常微分方程理论知识可知，上述问题存在唯一解()y x 。现在的目标就是计算区间[],a b 上等节点i x a ih =+处该未知函数的函数值()i y x ，其中()/h b a N =-，0,1,2,...,i N =。为此，可以用求解常微分方程问题的单步法，即欧拉方法、梯形方法和龙格_库塔方法求解。

3 算法求解和程序设计

3.1 欧拉方法

将微分方程离散化，用向前商

1()()

n n y x y x h

+-代替微分()n y x ，代入（1）中的微分方程，可得

1()()

(,())1,2,3,...n n n n y x y x f x y x n h

+-==，（）

化简可得

1()()(,())1,2,3,...n n n n y x y x f x y x h n +=+=，（）

如果用n y 近似()n y x 代入上式便可得到1()n y x +的近似值1n y +，计算式为：

1(,)1,2,3,...n n n n y y f x y h n +=+=，（）

(2) 这样问题(1)的近似解课通过求解下面的差分初值问题：

10(,)1,2,3,...()n n n n y y f x y h n y a y +=+=⎧⎨

=⎩

，（）

(3) 得到，按(3)式由初值0y 可逐次求出12y ,,...y 。

欧拉方法就是用差分方程初值问题(3)的解来近似微分方程初值问题(1)的解。即由公式(3)算出

()n y x 的近似值。这组公式求问题(1)的数值解就是著名的欧拉(Euler )公式。欧拉法是数值求解微分

方程最简单的方法，它特别适合快速编程，尽管精度不是很高。

Euler 方法的局部截断误差为

11()()(,())n n n n n T y x y x hf x y x ++=--

()()'()n n n y x h y x hy x =+-- 2

3

2

''()/2()()n h y x h h οο=+≈

即局部截断误差是2h 阶的。而数值算法的精度定义为：若一种算法的局部截断误差为1()p h ο+，则称该算法具有p 阶精度。显然p 越大，方法的精度越高。上式说明，欧拉方法是一阶方法，因此它的精度不高。

Euler 方法的整体截断误差为

()n n n e y x y =-

((1)1)/((1)1)n

n e hL hL T ≤+-+-

()

(1)/()()L b a e

hL T h ο-≤-=

即整体截断误差是h 阶的。 欧拉法流程图：

程序代码见附录。 3.2 梯形方法

在一阶微分方程111()()'()(,())n n n

n

x x n n x x y x y x y x dx f x y x dx +++-=

=

⎰

⎰

右端积分用梯形求积公式近

似，并用n y 代替()n y x ，1n y +代替1()n y x +

111((,)(,))2

n n n n n n h y y f x y f x y +++=+

+

梯形法的迭代公式为：

(0)1(1)()

111(,)((,)(,))2

n n n k k n n n n n n y f x y h h y y f x y f x y +++++⎧=⎪

⎨=++⎪⎩，(k=0,1,2,...) 梯形方法也是隐式方法，要迭代求解，迭代收敛的条件是

h 12

L <。

梯形法的误差估计的局部截断误差是3

()h ο，整体截断误差是2

()h ο。即梯形方法具有2阶精度，是

梯形法的流程图：

程序代码见附录。 3.3 龙格_库塔方法

由Lagrange 微分中值定理，

11()()'()()()(,())n n n n n y x y x y x x y x hf y ξξξ++=+-=+

记*

(,())k hf y ξξ=，则得到

*

1()()n n y x y x k +=+

这样，给出*

k 的一种算法，就得到求解微分方程初值问题的一种计算公式。