1

我正在尝试数值求解一个简单的二阶 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
4

1 回答 1

1

这是一个简单的复制粘贴错误。在该行x4=...中应该没有f，因为它在x1=..., x2=..., 中x3=...。

然后一切看起来都很好。

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

顺便说一句，如果您希望最后一点有 time 2*pi，您应该（在两种方法中）迭代range(n)以具有n集成步骤和T=n*(2*pi/n)=2*pi.

于 2015-04-29T14:47:30.590 回答