Python 中的数值 ODE 求解

Numerical ODE solving in Python

提问by dustin

How do I numerically solve an ODE in Python?

如何在 Python 中数值求解 ODE？

Consider

考虑

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

with the following conditions

具有以下条件

u(0) = 1.49907

and

和

\dot{u}(0) = 0

with the constraint

与约束

0 <= \phi <= 7\pi.

Then finally, I want to produce a parametric plot where the x and y coordinates are generated as a function of u.

最后，我想生成一个参数图，其中 x 和 y 坐标是作为 u 的函数生成的。

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.

问题是，我需要运行 odeint 两次，因为这是一个二阶微分方程。我尝试在第一次之后再次运行它，但它返回一个雅可比错误。必须有一种方法可以一次运行两次。

Here is the error:

这是错误：

odepack.error: The function and its Jacobian must be callable functions

odepack.error: 函数及其雅可比行列式必须是可调用函数

which the code below generates. The line in question is the sol = odeint.

下面的代码生成。有问题的行是 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()

Edit

编辑

I am trying to construct the plot on page 9

我正在尝试构建第 9 页的情节

Classical Mechanics Taylor

经典力学泰勒

Here is the plot with Mathematica

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

采纳答案by unutbu

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

回答by Bitwise

scipy.integrate()does ODE integration. Is that what you are looking for?

scipy.integrate() 进行ODE 集成。这就是你要找的吗？

回答by HenriV

You can use scipy.integrate.ode. To solve dy/dt = f(t,y), with initial condition y(t0)=y0, at time=t1 with 4th order Runge-Kutta you could do something like this:

您可以使用 scipy.integrate.ode。要求解 dy/dt = f(t,y)，初始条件为 y(t0)=y0，在 time=t1 和 4 阶 Runge-Kutta，您可以执行以下操作：

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

Edit: You have to get your derivative to first order to use numerical integration. This you can achieve by setting e.g. z1=u and z2=du/dt, after which you have dz1/dt = z2 and dz2/dt = d^2u/dt^2. Substitute these into your original equation, and simply iterate over the vector dZ/dt, which is first order.

编辑：你必须让你的导数达到一阶才能使用数值积分。这可以通过设置例如 z1=u 和 z2=du/dt 来实现，之后你有 dz1/dt = z2 和 dz2/dt = d^2u/dt^2。将这些代入原始方程，然后简单地迭代向量 dZ/dt，这是一阶。

Edit 2: Here's an example code for the whole thing:

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

回答by jorgeca

The code from your other questionis really close to what you want. Two changes are needed:

您另一个问题中的代码非常接近您想要的。需要做两个改动：

• You were solving a different ODE (because you changed two signs inside function deriv)
• The ycomponent of your desired plot comes from the solution values, not from the values of the first derivative of the solution, so you need to replace u[:,0](function values) for u[:, 1](derivatives).

• 您正在求解不同的 ODE（因为您更改了函数内部的两个符号deriv）
• y您所需绘图的组成部分来自解值，而不是来自解的一阶导数的值，因此您需要替换u[:,0]（函数值）为u[:, 1]（导数）。

This is the end result:

这是最终结果：

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

However, I suggest that you use the code from unutbu's answer because it's self documenting (u, udot = z) and uses np.linspaceinstead of np.arange. Then, run this to get your desired figure:

不过，我建议你使用的代码从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()