1

私は Python の初心者であり、プログラミング言語の知識はまだ初期段階にあるため、ここに示す Runge-Kutta Python スクリプトをコピーし、自分の目的に合わせて変更しました。これが私の現在のスクリプトです：

import numpy as np
from matplotlib import pyplot as plt
a=0
b=np.pi
g=9.8
l=1
N=1000

def RK4(f):
    return lambda t, y, dt: (
            lambda dy1: (
            lambda dy2: (
            lambda dy3: (
            lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6
            )( dt * f( t + dt  , y + dy3   ) )
        )( dt * f( t + dt/2, y + dy2/2 ) )
        )( dt * f( t + dt/2, y + dy1/2 ) )
        )( dt * f( t       , y         ) )

from math import sqrt
dy = RK4(lambda t, y: -y)

t, y, dt = 0., 1., np.divide((b-a),float(N))
i=0
T=np.zeros((N+1,1))
DY=T
Y=T
while t < (b-dt):
    T[i]=t
    DY[i]=dy(t,y,dt)
    t, y = t + dt, y + dy( t, y, dt )
    Y[i]=y
    i=i+1

plt.figure(1)
plt.plot(T,Y)
plt.show()

最初の数行の g 変数と l 変数は無視してかまいません。単純な振り子の問題を解くつもりでしたが、これが 1 次 ODE ソルバーであることを思い出したので、私の ODE はdy/dx=-y. これをIPythonで実行しています。基本的にMATLABTと同等の配列であると予想していました。linspace(0,pi,N+1)したがって、T は 0 と pi の間 (および 0 を含む) の N+1 個の等間隔の値のセットになりますが、代わりにこれはその内容のサンプルです (これは running からの出力として得られますT)。

In [101]: T
Out[101]:
array([[ 0.99686334],
       [ 0.99373651],
       [ 0.9906195 ],
       ...,
       [ 0.04334989],
       [ 0.04321392],
       [ 0.        ]])

（不明確な場合に備えて、ここで言及しているものについてのコンテキストを提供するために、入力行と出力行を含めます）。ああ、ちなみに、なぜこのループでそれを定義する代わりに使用しなかったのか疑問に思っているならT=np.linspace(a,b,num=N+1)、これは同様に珍しい T 配列を与えるからです。

4

1 に答える 1

3

と

DY=T
Y=T

参照をコピーするだけで、それらはすべて同じ配列オブジェクトを指しています。内容は基本的に全てYラストアサインのままとなります。

T.copy()別の配列オブジェクトを取得するために使用します。

于 2016-01-01T00:44:11.047 に答える