如何在代码中根据前一个 ODE 计算第二个 ODE

Question

首先，我将向您展示我的代码：

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

def acceleration(u, t):
  dudt = 3*pow(t,2)
  return dudt

t = np.linspace(0, 86400, 10) #time in s
u0 = 0 #velocity u for t = 0

u = odeint(acceleration, u0, t)

A = np.reshape(acceleration(u,t), (10,1))
print(A, u)

plt.xlabel('Time in s')
plt.ylabel('Velocity in m/s')
plt.plot(t, u)
plt.grid()
plt.show()

plt.xlabel('Time in s')
plt.ylabel('Acceleration in m/s^2')
plt.plot(t, A)
plt.grid()
plt.show()


def velocity(u, t):
   dxdt = u
   return dxdt
x0 = 0 #displacement for t = 0

x = odeint(velocity, x0, t)
print(x)

plt.xlabel('Time in s')
plt.ylabel('Displacement in m')
plt.plot(t, x)
plt.grid()
plt.show()

这是我试图编写的一个程序，用于解决基本物理问题的常微分方程（在本例中为位移、速度、加速度）。该代码一直运行到第二张图的点。当它尝试求解 x 时，我得到一个值为 0 的矩阵和一个包含 x 直线的图形，它等于 0。

我的问题是，如何通过对速度公式 u 积分来计算位移公式 x ？我可能错过了一些我不知道的东西，我尝试先通过 sympy 符号集成 u，但 jupyter 返回此错误：

ImportError: cannot import name 'symbol' from 'sympy' (unknown location)

提前感谢您的回答

Answer 1

如果我理解正确，你的 ODE 如下：

x''(t) = 3*t^2

根据 odeint 文档中的示例，我们需要将这个二阶 ODE 分解为两个一阶 ODE :

v'(t) = 3*t^2
x'(t) = v(t)

要按照您想要的方式使用

odeint

求解此系统，我们需要立即将 整个系统 传递给函数 :

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

# define the system
def acceleration(t):
    return 3*t/1e7

def system(u, t):
    v, x = u # split input variable
    dvdt = acceleration(t)
    dxdt = v 
    return [dvdt, dxdt]

# timeline
t = np.linspace(0, 86400, 10)  # time in s

# initial conditions
x0 = 0  # displacement for t = 0
v0 = 0  # velocity u for t = 0


# solve
sol = odeint(system, (v0, x0), t)
x = sol[:, 1]
v = sol[:, 0]
a = acceleration(t)

# plot
fig, ax = plt.subplots(3, 1)

ax[0].set_xlabel('Time in s')
ax[0].set_ylabel('Acceleration in m/s^2')
ax[0].plot(t, a)
ax[0].grid()

ax[1].set_xlabel('Time in s')
ax[1].set_ylabel('Velocity in m/s')
ax[1].plot(t, v)
ax[1].grid()

ax[2].set_xlabel('Time in s')
ax[2].set_ylabel('Displacement in m')
ax[2].plot(t, x)
ax[2].grid()

plt.show()

---

注意：对于较新的代码，scipy 鼓励使用

solve_ivp

而不是

odeint

。 所需的代码更改很少：

def system(t, u): # invert `u` and `t` args
    v, x = u
    dvdt = acceleration(t)
    dxdt = v
    return [dvdt, dxdt]


# solve
t_span = t.min(), t.max()
sol = solve_ivp(system, t_span, (v0, x0), t_eval=t) # slightly different function call

# `sol` provides much more information than just the solution. Have a look at the docs for a full description
x = sol.y[1, :] 
v = sol.y[0, :]
a = acceleration()
``