常微分方程数值计算方法

方法思路：

1. 每一次计算，都使用连个不同的计算方法求解，得出两个值
2. 如果两次求解结果接近，则接受求解结果。
3. 如果两次求解结果的误差超过了指定的容许值，则减小步长。如果误差超过了有效位数，则增加步长。

计算程序(待修改):

1. Runge-Kutta四阶经典算法。

def f(y):
    return 1   y ** 2def rk4(a, b, y_init,M):
    h = (b - a)/M
    y = y_init
    for i in range(M):
        k1 = f(y)
        k2 = f(y   h*k1/2)
        k3 = f(y   h*k2/2)
        k4 = f(y   h*k3)
        y = y   h*(k1   2*k2   2*k3   k4)/6
        print("iteration time {1:2d}: y is {0:.6f}".format(y, i 1))c = rk4(0, 1.4, 0, 14)打印结果:iteration time  1: y is 0.100335iteration time  2: y is 0.202710iteration time  3: y is 0.309336iteration time  4: y is 0.422793iteration time  5: y is 0.546302iteration time  6: y is 0.684137iteration time  7: y is 0.842289iteration time  8: y is 1.029639iteration time  9: y is 1.260159iteration time 10: y is 1.557406iteration time 11: y is 1.964747iteration time 12: y is 2.572072iteration time 13: y is 3.601563iteration time 14: y is 5.791975

2. Runge Kutta Fehlbrg方法

def fsolve(t, y):
    return 1   y ** 2def rkfMethod(a, b, y_init,M):
    h = (b - a)/M
    y = y_init
    t = a
    #-------k coefficient----------
    a1 = 1;     b1 = 1
    a2 = 1/4;   b2 = 1/4
    a3 = 3/8;   b3 = 3/32; c3 = 9/32
    a4 = 12/13; b4 = 1932/2197; c4 = -7200/2197; d4 = 7296/2197
    a5 = 1;     b5 = 439/216;   c5 = -8;         d5 = 3680/513;   e5 = -845/4104
    a6 = 1/2;   b6 = -8/27;     c6 = 2 ;         d6 = -3544/2565; e6 = 1859/4104; f6 = -11/40
    #--------y coefficient---------
    r1 = 16/135; r3 = 6656/12825; r4 = 28561/56430; r5 = -9/50; r6 = 2/55

    for i in range(M):
        k1 = fsolve(t, y)
        k2 = fsolve(t   a2*h, y   b2*k1*h)
        k3 = fsolve(t   a3*h, y   b3*k1*h   c3*k2*h)
        k4 = fsolve(t   a4*h, y   b4*k1*h   c4*k2*h   d4*k3*h)
        k5 = fsolve(t   a5*h, y   b5*k1*h   c5*k2*h   d5*k3*h   e5*k4*h)
        k6 = fsolve(t   a6*h, y   b6*k1*h   c6*k2*h   d6*k3*h   e6*k4*h   f6*k5*h)
        #------------------------------
        y = y   h*(r1*k1   r3*k3   r4*k4   r5*k5   r6*k6)
        print("iteration time {1:2d}: y is {0:.6f}".format(y, i 1))
    return y打印结果：rkfMethod(0, 1.4, 0, 14)iteration time  1: y is 0.100335iteration time  2: y is 0.202710iteration time  3: y is 0.309336iteration time  4: y is 0.422793iteration time  5: y is 0.546303iteration time  6: y is 0.684137iteration time  7: y is 0.842288iteration time  8: y is 1.029639iteration time  9: y is 1.260159iteration time 10: y is 1.557409iteration time 11: y is 1.964762iteration time 12: y is 2.572157iteration time 13: y is 3.602127iteration time 14: y is 5.798128

参考文献：

1. John.H.Mathews. 数值方法(Matlab)第四版.北京：机械工业出版社. 2017.02

2. Ramin S. Esfandiari. Numerical Methods for Engineers and Scientists Using MATLAB. CRC Press. 2017

最早发布于 2018-03-04