常微分方程初值问题

常微分方程描述了不同变量之间的变化关系，通过该关系我们可以确定变量之间的具体函数关系，但是微分方程不一定总是可解的，并且有些求解起来相当困难。给出积分曲线初始位置的状态，求解需要预测点的状态，为常微分方程的初值问题，我们可以在求解微分方程即的情况下对结果做出复合精度要求的预测。

以$f(x) = \sqrt{1+2x}$为例

有$y’=\frac{dy}{dx}=y-\frac{2x}{y}$

本文给出几种算法的原理与实现。

Euler方法及其改进

从初始点开始，根据不同离散点的导数值对曲线进行预测，导数值可以通过对ODE的化简求解，做出一条折线图，最终曲线会逼近预测值。

该方法有明显的缺点就是在每一步做出抉择时，只考虑了当前的状态，并没有考虑后面的状态，因此必然会造成较大的误差，因此采用下一个结点的导数值进行修正。

具体流程如下

选取一定的步长h，一般为区间的n等分

预测下一个值$\overline{y}_{i+1} = y_i + hf(x_i, y_i)$

校正预测值$y_{i+1} = y_i + \frac{h}{2}[f(x_i, y_i)+f(x_{i+1}, \overline{y}_{i+1})]$

改写成平均化形式

$y_p = y_i +hf(x_i, y_i)$

$y_c = y_i + hf(x_{i+1}, y_p)$

$y_{i+1} = \frac{1}{2} (y_p+y_c), i =0,1,2,3,…,n-1$

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37


import numpy as np

# 常微分方程化简
def ode(x, y):
    return y - 2 * x / y

def f(x):
    return np.sqrt(1+2*x)

def Euler(bound, start):
    '''
    bound : 预测边界值
    start : 起始点
    '''
    h = 0.1
    # 对步长h等分
    n = int(bound/h)
    x = np.zeros(n)
    y = np.zeros(n)
    sample = np.zeros(n)
    for i in range(0, n):
        # 获取均分后的所有点
        x[i] = 0+i*h
        # 生成标准值做对比
        sample[i] = f(x[i])
    y[0] = start
    for i in range(1, n):
        # 预测值
        yp = y[i-1] + h*ode(x[i-1], y[i-1])
        yc = y[i-1] + h*ode(x[i], yp)
        # 使用预测值求加权平均值
        y[i] = (yp+yc)/2.0
    return y, sample

answer,sample = Euler(2, 1)
print(answer)
print(sample)

运行结果：

标准值

[1. 1.09544512 1.18321596 1.26491106 1.34164079 1.41421356 1.4832397 1.54919334 1.61245155 1.67332005 1.73205081 1.78885438 1.84390889 1.8973666 1.94935887 2. 2.04939015 2.0976177 2.14476106 2.19089023]

预测值

[1. 1.09590909 1.18409657 1.26620136 1.34336015 1.41640193 1.4859556 1.55251409 1.61647478 1.67816636 1.7378674 1.79581974 1.8522386 1.90732042 1.96124939 2.01420304 2.06635728 2.11789132 2.16899248 2.21986124]

Runge-Kuta法

Runge-Kuta法在思想上与之前的方法大致相同，都是通过若干个斜率的加权平均值来做出下一步选择。此处只给出常用的四阶Runge-Kuta的处理方式。

$K_1=f(x_i,y_i)$

$K_2 =f(x_i + \frac{h}{2}, y_i + \frac{h}{2}K_1)$

$K_3=f(x_i + \frac{h}{2}, y_i + \frac{h}{2}K_2)$

$K_4=f(x_i + h, y_i + hK_3)$

$y_{i+1}= y_i + \frac{h}{6}(K_1+K_2+K_3+K_4)$

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38


import numpy as np
from matplotlib import pyplot as plt

def ode(x, y):
    return y - 2 * x / y

def f(x):
    return np.sqrt(1+2*x)

def RungeKutta(bound, start):
    '''
    bound : 计算边界
    start : 初始值f(0)
    '''
    h = 0.1
    n = int(bound/h)
    x = np.zeros(n)
    y = np.zeros(n)
    sample = np.zeros(n)
    for i in range(0, n):
        x[i] = 0+i*h
        sample[i] = f(x[i])
    y[0] = start
    for i in range(1, n):
        # 求各个点的斜率
        k1 = ode(x[i-1], y[i-1])
        k2 = ode(x[i-1] + h/2, y[i-1] + h/2*k1)
        k3 = ode(x[i-1] + h/2, y[i-1] + h/2*k2)
        k4 = ode(x[i-1] + h, y[i-1] + h * k3)
        # 对斜率求平均值并进行预测
        y[i] = y[i-1] + h/6*(k1 + 2 * k2 + 2 * k3 + k4)
    
    plt.show()
    return y, sample

y, sample = RungeKutta(2, 1)
print(y)
print(sample)

运行结果：预测值

[1. 1.09544553 1.18321675 1.26491223 1.34164235 1.41421558 1.48324222 1.54919645 1.61245535 1.67332466 1.73205637 1.78886107 1.84391692 1.89737622 1.9493704 2.00001382 2.0494067 2.09763752 2.14478481 2.1909187 ]

标准值

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

常微分方程初值问题

Euler方法及其改进

Runge-Kuta法