lunge-kutta法matlab,数值计算方法—微分方程(组).doc

数值计算方法—微分方程(组)

科学计算—理论、方法

及其基于MATLAB的程序实现与分析

微分方程(组)数值解法

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

微分方程(组)是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton第二定律：

(１)

和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题：

(２)

的解为：

(３)

但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：

１、绝大多数的常微分方程和常微分方程组得不到(有限形式的)解析解；

２、实际应用中往往只需要知道常微分方程(组)的解在(人们所关心的)某些点处的函数值(可以是满足一定精度要求的近似值)；

如果只需要常微分方程(组)的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程(组)的初值问题的数值解法，就可以达到这一目的。

一般的一阶常微分方程(组)的初值问题是指如下的一阶常微分方程(组)的定解问题：

(７)

其中

(８)

(９)

常微分方程(组)的初值问题通常是对一动态过程(动态系统、动力系统)演化规律的描述，求解常微分方程(组)的初值问题就是要了解和掌握动态过程演化规律。

§１.1 常微分方程(组)的Cauch问题数值解法概论

假设要求在点(时刻)，处初值问题(７)的解的(近似)值,如果已求得时刻的值或它的近似值(如时刻的值)，那么将式(７)的两端在区间上积分

(10)

可得

(11)

或

(12)

显然,为了利用式(11)或(12)求得的精确值(近似值)，必须计算右端的积分，这是问题的关键也是难点所在，如前所述，一般得不到精确的公式解，因此需要采用数值积分的方法求其近似解,

可以说,不同的式值积分方法将给出不同的Cauch问题的数值解法。

§1.2　最简单的数值解法——Euler 方法

假设要求在点(时刻)，，处初值问题(７)的解的近似值。首先对式(７)的两端积分，得

(13)

对于式(13)的右边，如果用积分下限处的函数值代替被积函数作积分(从几何上的角度看，是用矩形面积代替曲边梯形面

积)，则有

(14)

进而得到下式给出的递推算法—Euler 方法

(15)

例１　用Euler 方法解如下初值问题，取，

解：由(15)得

结果如下：

open Euler_Method.m

如果取，其结果如下图所示：

Euler_Method

§1.3　改进的Euler 方法

对于(15)的右边，如果被积函数用积分限和处的函数值的算术平均值代替(几何上，是用梯形面积代替曲边梯形面积)，则有

(16)

进而得到下式给出的递推算法：

(17)

通常算法(17)比Euler 方法(15)的精度高，但是，按算法(17)求时要解(非线性)方程(组)，这是算法(17)不如Euler 方法的方面，为了

尽可能地保持算法(17)精度高的优点；

尽可能地利用Euler 方法计算简单的长处；

人们采取了如下的称之为改进的Euler 方法的折衷方案：

预测　 　 　　　　　　(18)

修正　 　　　 (19)

例２　Euler 方法与改进的Euler 方法的比较

下图是当时比较的结果：

open Improved_Euler_Method.m

§1.4　Euler 方法和改进的Euler 方法的误差分析

由Taylor 公式

(19)

说明Euler 方法的截断误差是，类似地，由

(20)

(21)

以及

(22)

让式(20)的两端减式(21)的两端，可得

(23)

从上述推导Euler 方法、改进的Euler 方法的过程以及例１、例２容易看出，改进的Euler 方法Euler 方法的精度高，其原因在于：

在推导Euler 方法时，我们是用待求解函数在一点处的变化

率代替在区间上的平均变化率：

(24)

2 而在推导改进的Euler 方法时，我们是用待求解函数在两点处变化率的平均值代替在区间上的平均变化率；

显然，通常比更接近于在区间上的平均变化率。由此启发人们：适当地选取区间上函数若干点处的变化率，用它们加权平均值代替在区间上