1

単純な 2 次 DE x'' = -x を数値的に解こうとしています。新しい変数 x'=v を使用したので、2 つの方程式ができました。シンプルに見えますが、どうにか正しい結果とはかけ離れた結果​​が得られます。

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=...は はありません。fx1=...x2=...x3=...

その後、すべてが正常に見えます。

これは、クックブックの実装に可能な限り近づけ、ベクトル値関数とベクトル演算を使用する理由の 1 つです。

ところで、最後の点に時間を持たせたい場合は、2*pi(両方の方法で) 反復range(n)してn積分ステップとを作成する必要がありますT=n*(2*pi/n)=2*pi。

于 2015-04-29T14:47:30.590 に答える