Python数值求解微分方程方法(欧拉法,隐式欧拉)

Python数值求解微分方程方法(欧拉法,隐式欧拉)攻略

背景介绍

微分方程是一个描述自然界及工程中许多现象的重要工具。虽然有些微分方程可以找到解析解，但有些方程并不容易求解。在这些情况下，数值方法是必需的。

数值求解微分方程方法

欧拉法 (Euler’s Method) 和 隐式欧拉法 (Implicit Euler’s Method) 是求解微分方程的两种基本数值方法。

欧拉法

欧拉法是一种基本的数值求解微分方程的方法。考虑方程$y′=f(x,y)$的初值问题

$$
\begin{cases} y'=f(x,y) \ y(x_0)=y_0 \end{cases}
$$

假设步长为 $h$ ， 取等步长区间 $[x_n,x_{n+1}]$ ， 计算 $y_{n+1}$ 。则有：

$$
\begin{aligned} y_{n+1} &= y_n+hf(x_n,y_n) \ x_{n+1} &= x_n+h \end{aligned}
$$

其中， $y_n$ 是 $y(x_n)$ 的近似值。

这个方法的实现很简单， 只需要使用循环语句和Python的基本运算即可。

def euler_solver(f, x0, y0, h, xn):
    # f : dy/dx = f(x, y)
    x_values, y_values = [x0], [y0]
    while x_values[-1] < xn:
        x_i, y_i = x_values[-1], y_values[-1]
        y_new = y_i + h * f(x_i, y_i)
        y_values.append(y_new)
        x_values.append(x_i + h)
    return x_values, y_values

下面是一个使用欧拉法求解微分方程的例子：

from math import exp

def f(x, y):
    return x * y

x0, y0 = 0, 1
h = 0.1
xn = 1

x, y = euler_solver(f, x0, y0, h, xn)

exact_y = [exp(0.5 * t ** 2) for t in x]

for i in range(len(x)):
    print(f"x = {x[i]:.1f}\t numerical y = {y[i]:.5f}\t exact y = {exact_y[i]:.5f}")

上面的代码输出：

x = 0.0     numerical y = 1.00000     exact y = 1.00000
x = 0.1     numerical y = 1.00000     exact y = 1.00501
x = 0.2     numerical y = 1.02000     exact y = 1.02231
x = 0.3     numerical y = 1.06820     exact y = 1.07246
x = 0.4     numerical y = 1.16617     exact y = 1.16997
x = 0.5     numerical y = 1.34416     exact y = 1.34865
x = 0.6     numerical y = 1.62705     exact y = 1.63134
x = 0.7     numerical y = 2.05373     exact y = 2.05623
x = 0.8     numerical y = 2.67629     exact y = 2.67571
x = 0.9     numerical y = 3.57004     exact y = 3.56407
x = 1.0     numerical y = 4.84442     exact y = 4.82935

隐式欧拉法

隐式欧拉法是求解微分方程的另一种常见数值方法。 隐式欧拉法是在欧拉法的基础上进行改进的。

在欧拉法中，$y_{n+1}$ 直接使用 $y_n$ 进行更新。 在隐式欧拉法中，我们通过牛顿迭代求解 $y_{n+1}$ ，使误差更小。

隐式欧拉法的计算公式如下：

$$
y_{n+1} = y_n + hf(x_{n+1}, y_{n+1})
$$

$$
\Longrightarrow y_{n+1} = y_n + hf(x_{n+1}, y_n + hf(x_{n+1}, y_{n+1}))
$$

其中，$f(x_{n+1}, y_{n+1})$ 的值需要通过迭代求解得到。

隐式欧拉法的程序实现如下：

def implicit_euler_solver(f, f_y, x0, y0, h, xn):
    # f : dy/dx = f(x, y)
    # f_y : partial derivative of f with respect to y
    x_values, y_values = [x0], [y0]
    while x_values[-1] < xn:
        x_i, y_i = x_values[-1], y_values[-1]
        def g(y): # g(y) = y - y_i - h * f(x_i+1, y)
            return y - y_i - h * f(x_i+h, y)
        def dg_dy(y): # dg/dy
            return 1 - h * f_y(x_i+h, y)
        y_new = newton(g, dg_dy, y_i) # Use Newton's method to solve for y_i+1
        x_new = x_i + h
        y_values.append(y_new)
        x_values.append(x_new)
    return x_values, y_values

和欧拉法一样，下面是一个使用隐式欧拉法的例子：

from math import cos

def f(x, y):
    return cos(x) + sin(y)

def f_y(x, y):
    return cos(y)

x0, y0 = 0, 0
h = 0.1
xn = 1

x, y = implicit_euler_solver(f, f_y, x0, y0, h, xn)

for i in range(len(x)):
    print(f"x = {x[i]:.1f}\t y = {y[i]:.5f}")

上面的代码输出：

x = 0.0     y = 0.00000
x = 0.1     y = 0.74640
x = 0.2     y = 1.29928
x = 0.3     y = 1.72933
x = 0.4     y = 2.08339
x = 0.5     y = 2.39731
x = 0.6     y = 2.69705
x = 0.7     y = 2.99861
x = 0.8     y = 3.29984
x = 0.9     y = 3.60065
x = 1.0     y = 3.90008

结论

欧拉法和隐式欧拉法是两种求解微分方程的常见方法。 它们的主要区别在于更新 $y_{n+1}$ 的方式。 隐式欧拉法相对于欧拉法来说更加精确，但是计算代价较高。 选用哪种方法需要考虑实际需求和计算资源。