14

PythonでODEを数値的に解くにはどうすればよいですか?

検討

解く方程式

\ddot{u}(\phi) = -u + \sqrt{u}

以下の条件で

u(0) = 1.49907

と

\dot{u}(0) = 0

制約付き

0 <= \phi <= 7\pi.

最後に、x 座標と y 座標が u の関数として生成されるパラメトリック プロットを作成します。

問題は、これは 2 次微分方程式であるため、odeint を 2 回実行する必要があることです。初めて実行した後、もう一度実行しようとしましたが、ヤコビアン エラーが返されます。一度に 2 回実行する方法が必要です。

エラーは次のとおりです。

odepack.error: 関数とそのヤコビアンは呼び出し可能な関数でなければなりません

以下のコードが生成します。問題の行は sol = odeint です。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from numpy import linspace


def f(u, t):
    return -u + np.sqrt(u)


times = linspace(0.0001, 7 * np.pi, 1000)
y0 = 1.49907
yprime0 = 0
yvals = odeint(f, yprime0, times)

sol = odeint(yvals, y0, times)

x = 1 / sol * np.cos(times)
y = 1 / sol * np.sin(times)

plot(x,y)

plt.show()

編集

9ページのプロットを構築しようとしています

古典力学テイラー

ここに Mathematica のプロットがあります

プロット

In[27]:= sol = 
 NDSolve[{y''[t] == -y[t] + Sqrt[y[t]], y[0] == 1/.66707928, 
   y'[0] == 0}, y, {t, 0, 10*\[Pi]}];

In[28]:= ysol = y[t] /. sol[[1]];

In[30]:= ParametricPlot[{1/ysol*Cos[t], 1/ysol*Sin[t]}, {t, 0, 
  7 \[Pi]}, PlotRange -> {{-2, 2}, {-2.5, 2.5}}]
4

4 に答える 4

27 
import scipy.integrate as integrate
import matplotlib.pyplot as plt
import numpy as np

pi = np.pi
sqrt = np.sqrt
cos = np.cos
sin = np.sin

def deriv_z(z, phi):
    u, udot = z
    return [udot, -u + sqrt(u)]

phi = np.linspace(0, 7.0*pi, 2000)
zinit = [1.49907, 0]
z = integrate.odeint(deriv_z, zinit, phi)
u, udot = z.T
# plt.plot(phi, u)
fig, ax = plt.subplots()
ax.plot(1/u*cos(phi), 1/u*sin(phi))
ax.set_aspect('equal')
plt.grid(True)
plt.show()

ここに画像の説明を入力

于 2013-04-10T14:56:33.017 に答える
4

あなたの他の質問のコードは、あなたが望むものに本当に近いです。2 つの変更が必要です。

  • 別の ODE を解いていました ( function 内の 2 つの記号を変更したためderiv)
  • 目的のプロットのyコンポーネントは、ソリューションの 1 次導関数の値ではなく、ソリューションの値に由来するため、u[:,0](関数の値) をu[:, 1](導関数) に置き換える必要があります。

これが最終結果です。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def deriv(u, t):
    return np.array([u[1], -u[0] + np.sqrt(u[0])])

time = np.arange(0.01, 7 * np.pi, 0.0001)
uinit = np.array([1.49907, 0])
u = odeint(deriv, uinit, time)

x = 1 / u[:, 0] * np.cos(time)
y = 1 / u[:, 0] * np.sin(time)

plt.plot(x, y)
plt.show()

ただし、unutbuの回答のコードを使用することをお勧めします。これは、自己文書化 ( u, udot = z) でありnp.linspace、np.arange. 次に、これを実行して目的の図を取得します。

x = 1 / u * np.cos(phi)
y = 1 / u * np.sin(phi)
plt.plot(x, y)
plt.show()
于 2013-04-11T21:43:32.067 に答える
3

scipy.integrate.ode を使用できます。dy/dt = f(t,y) を初期条件 y(t0)=y0 で、time=t1 で 4 次ルンゲクッタで解くには、次のようにします。

from scipy.integrate import ode
solver = ode(f).set_integrator('dopri5')
solver.set_initial_value(y0, t0)
dt = 0.1
while t < t1:
    y = solver.integrate(t+dt)
    t += dt

編集:数値積分を使用するには、導関数を最初の順序にする必要があります。これは、たとえば z1=u および z2=du/dt を設定することで実現できます。その後、dz1/dt = z2 および dz2/dt = d^2u/dt^2 が得られます。これらを元の方程式に代入し、1 次のベクトル dZ/dt を反復処理します。

編集 2: 全体のコード例を次に示します。

import numpy as np
import matplotlib.pyplot as plt

from numpy import sqrt, pi, sin, cos
from scipy.integrate import ode

# use z = [z1, z2] = [u, u']
# and then f = z' = [u', u''] = [z2, -z1+sqrt(z1)]
def f(phi, z):
    return [z[1], -z[0]+sqrt(z[0])]


# initialize the 4th order Runge-Kutta solver
solver = ode(f).set_integrator('dopri5')

# initial value
z0 = [1.49907, 0.]
solver.set_initial_value(z0)

values = 1000
phi = np.linspace(0.0001, 7.*pi, values)
u = np.zeros(values)

for ii in range(values):
    u[ii] = solver.integrate(phi[ii])[0] #z[0]=u

x = 1. / u * cos(phi)
y = 1. / u * sin(phi)

plt.figure()
plt.plot(x,y)
plt.grid()
plt.show()
于 2013-04-10T14:49:45.917 に答える