python-3.x - 为什么我不能让这个 Runge-Kutta 求解器随着时间步长的减少而收敛？

Question

由于某些原因，我需要在 PyTorch 中实现 Runge-Kutta4 方法(所以不，我不会使用 scipy.odeint)。我尝试并在最简单的测试用例上得到了奇怪的结果，用 x(0)=1 求解 x'=x(解析解:x=exp(t))。基本上，当我减少时间步长时，我无法降低数值误差。我可以使用更简单的 Euler 方法，但不能使用 Runge-Kutta 4 方法，这让我怀疑这里存在一些浮点问题(也许我遗漏了一些从 double 到单精度的隐藏转换)？

import torch
import numpy as np
import matplotlib.pyplot as plt

def Euler(f, IC, time_grid):
    y0 = torch.tensor([IC])
    time_grid = time_grid.to(y0[0])
    values = y0

    for i in range(0, time_grid.shape[0] - 1):
        t_i = time_grid[i]
        t_next = time_grid[i+1]
        y_i = values[i]
        dt = t_next - t_i
        dy = f(t_i, y_i) * dt
        y_next = y_i + dy
        y_next = y_next.unsqueeze(0)
        values = torch.cat((values, y_next), dim=0)

    return values

def RungeKutta4(f, IC, time_grid):

    y0 = torch.tensor([IC])
    time_grid = time_grid.to(y0[0])
    values = y0

    for i in range(0, time_grid.shape[0] - 1):
        t_i = time_grid[i]
        t_next = time_grid[i+1]
        y_i = values[i]
        dt = t_next - t_i
        dtd2 = 0.5 * dt
        f1 = f(t_i, y_i)
        f2 = f(t_i + dtd2, y_i + dtd2 * f1)
        f3 = f(t_i + dtd2, y_i + dtd2 * f2)
        f4 = f(t_next, y_i + dt * f3)
        dy = 1/6 * dt * (f1 + 2 * (f2 + f3) +f4)
        y_next = y_i + dy
        y_next = y_next.unsqueeze(0)
        values = torch.cat((values, y_next), dim=0)

    return values

# differential equation
def f(T, X):
    return X 

# initial condition
IC = 1.

# integration interval
def integration_interval(steps, ND=1):
    return torch.linspace(0, ND, steps)

# analytical solution
def analytical_solution(t_range):
    return np.exp(t_range)

# test a numerical method
def test_method(method, t_range, analytical_solution):
    numerical_solution = method(f, IC, t_range)
    L_inf_err = torch.dist(numerical_solution, analytical_solution, float('inf'))
    return L_inf_err


if __name__ == '__main__':

    Euler_error = np.array([0.,0.,0.])
    RungeKutta4_error = np.array([0.,0.,0.])
    indices = np.arange(1, Euler_error.shape[0]+1)
    n_steps = np.power(10, indices)
    for i, n in np.ndenumerate(n_steps):
        t_range = integration_interval(steps=n)
        solution = analytical_solution(t_range)
        Euler_error[i] = test_method(Euler, t_range, solution).numpy()
        RungeKutta4_error[i] = test_method(RungeKutta4, t_range, solution).numpy()

    plots_path = "./plots"
    a = plt.figure()
    plt.xscale('log')
    plt.yscale('log')
    plt.plot(n_steps, Euler_error, label="Euler error", linestyle='-')
    plt.plot(n_steps, RungeKutta4_error, label="RungeKutta 4 error", linestyle='-.')
    plt.legend()
    plt.savefig(plots_path + "/errors.png")

结果:

如您所见，Euler 方法收敛(缓慢，正如一阶方法所预期的那样)。但是，随着时间步长越来越小，Runge-Kutta4 方法不会收敛。错误最初下降，然后再次上升。这里有什么问题？

Answer 1

原因确实是 float 精度问题。 torch 默认为单精度，所以一旦 truncation error变得足够小，总误差基本上由 roundoff error 决定，并通过增加步数 <=> 减少时间步进一步减少截断误差不会导致总误差的任何减少。

要解决此问题，我们需要为所有浮点 torch 张量和 numpy 数组强制执行 double 64 位 float 。请注意，执行此操作的正确方法是分别使用 torch.float64 和 np.float64 而不是例如 torch.double 和 np.double，因为 the former are fixed-sized float values, (always 64bit) while the latter depend on the machine and/or compiler .这是固定代码:

import torch
import numpy as np
import matplotlib.pyplot as plt

def Euler(f, IC, time_grid):

    y0 = torch.tensor([IC], dtype=torch.float64)
    time_grid = time_grid.to(y0[0])
    values = y0

    for i in range(0, time_grid.shape[0] - 1):
        t_i = time_grid[i]
        t_next = time_grid[i+1]
        y_i = values[i]
        dt = t_next - t_i
        dy = f(t_i, y_i) * dt
        y_next = y_i + dy
        y_next = y_next.unsqueeze(0)
        values = torch.cat((values, y_next), dim=0)

    return values

def RungeKutta4(f, IC, time_grid):

    y0 = torch.tensor([IC], dtype=torch.float64)
    time_grid = time_grid.to(y0[0])
    values = y0

    for i in range(0, time_grid.shape[0] - 1):
        t_i = time_grid[i]
        t_next = time_grid[i+1]
        y_i = values[i]
        dt = t_next - t_i
        dtd2 = 0.5 * dt
        f1 = f(t_i, y_i)
        f2 = f(t_i + dtd2, y_i + dtd2 * f1)
        f3 = f(t_i + dtd2, y_i + dtd2 * f2)
        f4 = f(t_next, y_i + dt * f3)
        dy = 1/6 * dt * (f1 + 2 * (f2 + f3) +f4)
        y_next = y_i + dy
        y_next = y_next.unsqueeze(0)
        values = torch.cat((values, y_next), dim=0)

    return values

    # differential equation
def f(T, X):
    return X 

# initial condition
IC = 1.

# integration interval
def integration_interval(steps, ND=1):
    return torch.linspace(0, ND, steps, dtype=torch.float64)

# analytical solution
def analytical_solution(t_range):
    return np.exp(t_range, dtype=np.float64)

# test a numerical method
def test_method(method, t_range, analytical_solution):
    numerical_solution = method(f, IC, t_range)
    L_inf_err = torch.dist(numerical_solution, analytical_solution, float('inf'))
    return L_inf_err


if __name__ == '__main__':

    Euler_error = np.array([0.,0.,0.], dtype=np.float64)
    RungeKutta4_error = np.array([0.,0.,0.], dtype=np.float64)
    indices = np.arange(1, Euler_error.shape[0]+1)
    n_steps = np.power(10, indices)
    for i, n in np.ndenumerate(n_steps):
        t_range = integration_interval(steps=n)
        solution = analytical_solution(t_range)
        Euler_error[i] = test_method(Euler, t_range, solution).numpy()
        RungeKutta4_error[i] = test_method(RungeKutta4, t_range, solution).numpy()

    plots_path = "./plots"
    a = plt.figure()
    plt.xscale('log')
    plt.yscale('log')
    plt.plot(n_steps, Euler_error, label="Euler error", linestyle='-')
    plt.plot(n_steps, RungeKutta4_error, label="RungeKutta 4 error", linestyle='-.')
    plt.legend()
    plt.savefig(plots_path + "/errors.png")

结果:

现在，随着我们减少时间步长，RungeKutta4 近似的误差会随着正确率的增加而减少。