问题描述

我正在实现一个RKF4(5)积分器，但我无法确定我的代码是否工作正常，并且我不理解本地截断错误，或者我的代码是否不工作。

对于代码块的大小，我深表歉意，但在这种情况下，最小可重现示例相当大。

import numpy as np

def RKF45(state, derivative_function, h):
    """
    Calculate the next state with the 4th-order calculation, along with a 5th-order error
    check.

    Inputs:
    state: the current value of the function, float
    derivative_function: A function which takes a state (given as a float)
                         and returns the derivative of a function at that point
    h: step size, float
    """
    k1 = h * derivative_function(state)
    k2 = h * derivative_function(state + (k1 / 4))
    k3 = h * derivative_function(state + (k1 * (3/32)) + (k2 * (9/32)))
    k4 = h * derivative_function(state + (k1 * (1932/2197)) + (k2 * (-7200/2197)) + (k3 * (7296/2197)))
    k5 = h * derivative_function(state + (k1 * (439/216)) + (k2 * (-8)) + (k3 * (3680/513)) + (k4 * (-845/4104)))
    k6 = h * derivative_function(state + (k1 * (-8/27)) + (k2 * (2)) + (k3 * (-3544/2565)) + (k4 * (1859/4104)) + (k5 * (-11/40)))
    y1 = state + ((25/216) * k1) + ((1408/2565) * k3) + ((2197/4101) * k4) - ((1/5)*k5)
    y2 = state + ((16/135) * k1) + ((6656/12825) * k3) + ((28561/56430) * k4) - ((9/50) * k5) + ((2/55) * k6)
    return(y1, y2)


def integrate_RKF45(t0, tmax, tol, h_init, x_0, df, verbose = False):
    """
    integrate a function whose derivative is df from t0 to tmax
    t0: starting time
    tmax: end time
    h_init: initial timestep
    x_0: starting position
    df: a function which takes x and returns the derivative of a function at x
    """
    h = h_init
    x_i = x_0
    t = t0
    while t < tmax:
        h = min(h, tmax - t)
        y1, y2 = RKF45(x_i, df, h)
        err_i = np.abs(y1 - y2)
        R = 2 * err_i / h
        delta = (tol/R)**(1/4)
        if verbose:
            print(f"t: {t:0.2e}, dt: {h:0.2e}, x: {x_i:0.2e}, err: {err_i:0.2e}")
        if err_i < tol:
            t += h
            x_i = y1
        elif err_i > tol:        
            h *= delta
    return(x_i)


def exponential(x_0, k=1):
    """
    A simple test function, this returns the input, so it'll integrate to e^x.
    """
    return(k * x_0)

if __name__ == "__main__":
    integrate_RKF45(t0 = 0., 
                    tmax = 0.15,
                    tol = 1e-4, 
                    h_init = 1e-2, 
                    x_0 = 1.,
                    df = exponential,
                    verbose=True)

所以，这个代码的作用是，它返回我给它的任何函数的积分的近似值。然而，局部截断误差似乎太大了。运行以上代码将输出：

t: 0.00e+00, dt: 1.00e-02, x: 1.00e+00, err: 3.95e-06
t: 1.00e-02, dt: 1.00e-02, x: 1.01e+00, err: 3.99e-06
t: 2.00e-02, dt: 1.00e-02, x: 1.02e+00, err: 4.03e-06
t: 3.00e-02, dt: 1.00e-02, x: 1.03e+00, err: 4.07e-06
t: 4.00e-02, dt: 1.00e-02, x: 1.04e+00, err: 4.11e-06
t: 5.00e-02, dt: 1.00e-02, x: 1.05e+00, err: 4.16e-06
t: 6.00e-02, dt: 1.00e-02, x: 1.06e+00, err: 4.20e-06
t: 7.00e-02, dt: 1.00e-02, x: 1.07e+00, err: 4.24e-06
t: 8.00e-02, dt: 1.00e-02, x: 1.08e+00, err: 4.28e-06
t: 9.00e-02, dt: 1.00e-02, x: 1.09e+00, err: 4.32e-06
t: 1.00e-01, dt: 1.00e-02, x: 1.11e+00, err: 4.37e-06
t: 1.10e-01, dt: 1.00e-02, x: 1.12e+00, err: 4.41e-06
t: 1.20e-01, dt: 1.00e-02, x: 1.13e+00, err: 4.46e-06
t: 1.30e-01, dt: 1.00e-02, x: 1.14e+00, err: 4.50e-06
t: 1.40e-01, dt: 1.00e-02, x: 1.15e+00, err: 4.55e-06

其中err值是四阶方法和五阶方法之间的差异。我的印象是，n^th阶方法有O(dt^(n+1))阶的局部截断误差，这意味着上面的积分误差应该是1e-9左右，而不是1e-6。

那么是我的代码错误还是我的理解错误？ 谢谢！

推荐答案

请参阅https://math.stackexchange.com/questions/2701385/adaptive-step-size-in-rk45-for-second-order-ode/2701678#2701678，您似乎对方法系数使用了相同的损坏来源。

y1中的分母4101错误，必须为4104。

增量因子应该稍微缓和一点，delta = (tol/R)**(1/5)或delta = (tol/R)**(1/6)，每一步都要应用，成功的也要应用。

本地错误err_i的参考误差为tol*h，这就是在R中除以h的原因。

这将在径向较少的迭代步骤中产生测试场景

t: 0.000000e+00, dt: 1.00e-02, x: 1.000000e+00, err: 1.28e-13
t: 1.000000e-02, dt: 1.40e-01, x: 1.010050e+00, err: 6.60e-08
t: 1.500000e-01, dt: 3.88e-01, x: 1.161834e+00

或稍长一点的时间间隔以查看步长控制器实际正在工作

t: 0.000000e+00, dt: 1.00e-02, x: 1.000000e+00, err: 1.28e-13
t: 1.000000e-02, dt: 2.27e-01, x: 1.010050e+00, err: 7.18e-07
t: 2.372490e-01, dt: 4.31e-01, x: 1.267757e+00, err: 2.02e-05
t: 6.680839e-01, dt: 4.76e-01, x: 1.950509e+00, err: 5.03e-05
t: 6.680839e-01, dt: 4.47e-01, x: 1.950509e+00, err: 3.73e-05
t: 1.115525e+00, dt: 3.84e-01, x: 3.051213e+00, err: 2.81e-05
t: 1.500000e+00, dt: 3.89e-01, x: 4.481769e+00

​ 所有提到的更正都给出了RKF45中的新循环

    while t < tmax:
        h = min(h, tmax - t)
        y1, y2 = RKF45(x_i, df, h)
        err_i = np.abs(y1 - y2)
        R = err_i / h
        delta = 0.95*(tol/R)**(1/5)
        if verbose:
            print(f"t: {t:0.6e}, dt: {h:0.2e}, x: {x_i:0.6e}, err: {err_i:0.2e}")
        if R < tol:
            t += h
            x_i = y1
        h *= delta
    if verbose:
        print(f"t: {t:0.6e}, dt: {h:0.2e}, x: {x_i:0.6e}")

这篇关于理解自适应龙格库塔积分器的局部截断误差的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持编程学习网！