-1

私の2つの一次微分は次のとおりです

y1' = -sin(y0) + (gamma)*cos(y0)sin(beta * x)

と

y0' = y1

どこ

(theta)'' = y, (theta)' = y1, theta = y0

私の元の方程式は

(((d^2)*theta)/dt^2)=-sin(theta)+(gamma)cos(theta)sin(Bx)

時間の関数としてシータを解き、t=0 から t=40 までプロットするにはどうすればよいですか。システムは で静止状態から始まりtheta = 0 and d(theta)/dt = 0ます。

4

1 に答える 1

0

振動する外力で強制物理振り子をシミュレートしているようです。

theta''+sin(theta) = gamma * cos(theta)*sin(beta*t)

あなたが正しく認識しているように、これは次数 1 のシステムに変換され、ベクトル値の ODE 関数としてパッケージ化される必要があります。

def odefunc(t,y):
    y0, y1 = y
    return y1, -sin(y0)+gamma*cos(y0)*sin(beta*t)

定数がグローバル変数の場合、またはパラメーターとしての場合

def odefunc(t,y,params):
    y0, y1 = y
    beta, gamma = params
    return y1, -sin(y0)+gamma*cos(y0)*sin(beta*t)

次に、ラムダ式を使用して、引数を ODE インテグレーターの標準形式に減らす必要があります。

于 2015-12-07T10:34:54.227 に答える