python解常微分方程龙格库_数值常微分方程-欧拉法与龙格-库塔法

大三时候在跳蚤市场闲逛，从一位数学院的学长那里买了一些闲书，最近翻出来刚好有李荣华、刘播老师的《微分方程数值解法》和王仁宏老师的《数值逼近》，结合周善贵老师的《计算物理》课程，整理一下笔记。

本文整理常微分方程数值求解的欧拉法与龙格-库塔法。

一般地，动力学系统的时间演化可以用常微分方程的初值问题来描述，例如设一维简谐运动的回复力：

，有则运动方程：

。令

，可以将二阶微分方程转化为一阶微分方程组：

因此本文主要整理一阶常微分方程初值问题的数值解法。

一阶常微分方程初值问题

设

在区域

：

上连续，对于一个给定的常微分方程

及初值

，求解

。为了保证解

存在、唯一且连续依赖初值

，要求

满足Lipschitz条件：

存在常数L，使得

对所有

和

成立。

假设

总满足上述条件。常用的近似解法有级数解法等近似解析方法，以及下文整理的数值方法：欧拉法与龙格-库塔法。

欧拉法

将区间

作N等分，每一小区间长度

称为步长，

称为节点。根据初值

，代入微分方程可直接解出

的导数值

。

推导

1、根据泰勒展开式：

略去二阶小量，得：

以此类推，得到递推公式：

2、数值积分推导

由

可得：

，使用左矩形积分得：

以此类推，可得到：

为了提高精度，可以使用梯形积分代替矩形积分，即：

以此类推，得到改进的欧拉法：

Python计算实例

以

为例，其精确解为

，使用欧拉法求解的Python代码如下：

import math

from matplotlib import pyplot as plt

t_0 = 0

y_0 = 1

tau = 0.1

i = 1

solve = []

Euler = []

t = []

while i < 100:

if i == 1:

y_n = y_0

t_n = t_0

Euler.append(y_n)

solve.append(math.exp(t_n))

t.append(t_n)

func = y_n

y_n = y_n + tau * func

t_n = t_n + tau

i += 1

plt.plot(t, Euler, c='green', label=' Euler method')

plt.plot(t, solve, c='red', label=' accuracy')

plt.fill_between(t, solve, Euler, facecolor='blue', alpha=0.2)

plt.title('Euler method', fontsize=19)

plt.xlabel('t', fontsize=19)

plt.ylabel('y', fontsize=19)

plt.legend()

plt.show()

作图可以看到，当迭代步数较多后，欧拉法的结果逐渐落后于精确指数解的增长速度。下面分析欧拉法的误差来源。

在

中，略去了高阶小量

，因此在每一步的递推中，都有局部截断误差

，其阶为

。

在计算中，我们更关心精确解和数值解之间的误差

，称为整体误差，其满足

根据Lipschitz条件，可得：

由

且

，可得：

局部截断误差

是

的二阶量，设

，得：

，整体误差是

的一阶量。同理可得，改进的欧拉法局部截断误差

是

的三阶量

，整体误差是

的二阶量

。

稳定性分析

如果计算的初值不能精确给定，例如存在测量、舍入误差等，在计算过程中，每一步传递的误差连续依赖于初始误差，则称算法稳定，否则该算法不稳定。

对于不同的初值

和

，有

两式相减，得：

根据Lipschitz条件，可得：

连续依赖于初始误差，欧拉法稳定。同理，改进的欧拉法也稳定。

龙格-库塔法

龙格库塔法的主要思想：在

点的附近选取一些特定的点，然后把这些点的函数值进行线性组合，使用组合值代替泰勒展开中

点的导数值。

将

泰勒展开：

令

，根据多元函数求导法则有：

以此类推，可以得到：

同时，我们可以写出泰勒展开的形式解:

其中：

通项为：

基本思路是，利用当前点的函数值

可以计算出

，然后引入参数

和步长

可以计算出

，之后使用

和步长

计算

，以此类推，直到

。

现在把

展开：

将

代入得：

将

代入

可得：

与泰勒展开式

相比较，可知：

2个方程有3个未知数，因此有无穷多个解，可采用

，

，

，则：

令

，

，可以改写为：

此即为二阶龙格-库塔法。

与上一节的欧拉法公式对比：

，因此二阶龙格-库塔法取参数

，

，

时，即为改进的欧拉法。

Python计算实例

仍以

为例，其精确解为

，使用二阶龙格-库塔法求解的Python代码如下：

import math

from matplotlib import pyplot as plt

t_0 = 0

y_0 = 1

z_0 = 1

tau = 0.1

i = 1

j = 1

solve = []

Euler = []

R_K = []

t = []

while i < 100:

if i == 1:

y_n = y_0

t_n = t_0

R_K.append(y_n)

solve.append(math.exp(t_n))

t.append(t_n)

func_n = y_n

func_m = y_n + tau * func_n

y_n = y_n + 0.5 * tau * (func_n + func_m)

t_n = t_n + tau

i += 1

t = []

while j < 100:

if j == 1:

z_n = z_0

t_n = t_0

Euler.append(z_n)

t.append(t_n)

func = z_n

z_n = z_n + tau * func

t_n = t_n + tau

j += 1

plt.scatter(t, R_K, marker='^', c='blue', s=70, label=' R-K method')

plt.plot(t, Euler, c='green', label=' Euler method')

plt.plot(t, solve, c='red', label=' accuracy')

plt.fill_between(t, solve, Euler, facecolor='yellow', alpha=0.2)

plt.title('Euler method & R-K method', fontsize=19)

plt.xlabel('t', fontsize=19)

plt.ylabel('y', fontsize=19)

plt.legend()

plt.show()

黄色部分表示数值解和精确解的偏离，可以看到，二阶龙格-库塔法(改进的欧拉法)精确度得到了很大的提升。

二阶龙格-库塔法中，泰勒展开到了

阶，通过与泰勒展开系数进行对比，可以得到含3个未知数的2个方程。依次类推，如果泰勒展开到了

阶，对比

可以得到

级

阶龙格-库塔法。常用经典四阶龙格-库塔法：

Reference：

1、周善贵，《计算物理》课程讲义

2、李荣华，刘播，《微分方程数值解法》

3、王仁宏，《数值逼近》