python - Python中的数值ODE求解

如何在python中数值求解一个ode？
考虑

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

具备以下条件

u(0) = 1.49907

和

\dot{u}(0) = 0

有约束的

0 <= \phi <= 7\pi.

最后，我要生成一个参数图，其中x和y坐标是作为u的函数生成的。
问题是，我需要运行两次，因为这是一个二阶微分方程。
我试着让它在第一次运行后再次运行，但它回来时出现了雅可比错误。一定有办法一次运行两次。
错误如下：
odepack.error:函数及其Jacobian必须是可调用函数
下面的代码生成的。有问题的行是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页上画出这个情节。
Classical Mechanics Taylor
这是有数学的图

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}}]

最佳答案

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()

关于python - Python中的数值ODE求解，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/15928750/