python - RK4 给出错误结果

Question

我正在尝试对一个简单的二阶 DE x'' = -x 进行数值求解。我使用了一个新变量 x'=v，所以我有两个方程。虽然看起来很简单，但它以某种方式产生的结果与正确结果相差甚远。

def f(x):
    return -x

def rk4(f=f,h=2*pi/100,x0=1,v0=0,t0=0,n=100):
    '''RK4'''
    v=[v0]
    x=[x0]
    for i in range(n-1):
        v1= h*f(x[-1])
        x1= h *(v[-1])
        v2= h*f(x[-1]+1/2*x1)
        x2= h *(v[-1]+1/2*v1)
        v3= h*f(x[-1]+1/2*x2)
        x3= h *(v[-1]+1/2*v2)
        v4= h*f(x[-1]+x3)
        x4= h *f(v[-1]+v3)
        v.append(v[-1]+1/6 *(v1+2*v2+2*v3+v4))
        x.append(x[-1]+1/6 *(x1+2*x2+2*x3+x4))
    return x,v

奇怪的是，如果我使用 RK45 系数，一切都会正常。知道可能出了什么问题吗？

def rk4(f=f,h=2*pi/100,x0=1,v0=0,t0=0,n=100):
    '''RK4'''
    v=[v0]
    x=[x0]

    for i in range(n-1):
        v1=h*f(x[-1])
        x1=h *v[-1]
        v2=h*f(x[-1]+1/4*x1)
        x2=h *(v[-1]+1/4*v1)
        v3=h*f(x[-1]+9/32*x2+3/32*x1)
        x3=h *(v[-1]+9/32*v2+3/32*v1)
        v4=h*f(x[-1]+7296/2197*x3 - 7200/2197 * x2 + 1932 / 2197 * x1)
        x4=h *(v[-1]+7296/2197*v3 - 7200/2197 * v2 + 1932 / 2197 * v1)
        v5=h*f(x[-1]-845/4104 * x4 +3680/513*x3 - 8 * x2 + 439 / 216 * x1)
        x5=h *(v[-1]-845/4104 * v4 +3680/513*v3 - 8 * v2 + 439 / 216 * v1)
        v.append(v[-1]+25/216*v1+1408/2565*v3+2197/4104*v4-1/5*v5)
        x.append(x[-1]+25/216*x1+1408/2565*x3+2197/4104*x4-1/5*x5)
    return x,v

Answer 1

这是一个简单的复制粘贴错误。在 x4=... 行中不应有 f，就像在 x1=...、x2= 中一样...，x3=...。

然后一切看起来都很好。

这是尽可能接近食谱实现、使用向量值函数和向量算术的原因之一。

顺便说一句，如果您希望最后一个点有时间 2*pi，您应该迭代(在两种方法中)range(n) 以获得n 集成步骤和 T=n*(2*pi/n)=2*pi。