用Python进行数值ODE求解 [英] Numerical ODE solving in Python
问题描述
如何用Python数值求解ODE?
How do I numerically solve an ODE in Python?
考虑
\ddot{u}(\phi) = -u + \sqrt{u}
符合以下条件
u(0) = 1.49907
和
\dot{u}(0) = 0
受约束
0 <= \phi <= 7\pi.
然后,最后,我想生成一个参数图,其中x和y坐标是根据u生成的.
Then finally, I want to produce a parametric plot where the x and y coordinates are generated as a function of u.
问题是,我需要运行odeint两次,因为这是一个二阶微分方程. 我尝试在第一次之后再次运行它,但是它又返回了Jacobian错误.必须有一种方法可以一次运行两次.
The problem is, I need to run odeint twice since this is a second order differential equation. I tried having it run again after the first time but it comes back with a Jacobian error. There must be a way to run it twice all at once.
这是错误:
odepack.error:该函数及其Jacobian函数必须是可调用函数
odepack.error: The function and its Jacobian must be callable functions
下面的代码生成的
.问题所在的行是sol = odeint.
which the code below generates. The line in question is the 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页上绘制情节
I am trying to construct the plot on page 9
这是Mathematica的情节
Here is the plot with 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}}]
推荐答案
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进行数值ODE求解的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!