RK4给出错误的结果

计算科学 Python 龙格库塔
2021-12-05 00:41:49

我正在尝试数值求解一个简单的二阶微分方程我使用了一个新变量,所以我有两个方程。虽然看起来很简单,但它以某种方式产生的结果与正确结果相差甚远。x=xx=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
1个回答

正如在https://stackoverflow.com/a/29947165/3088138中发现的那样,该行中有一个简单的错字x4=...应该没有f