如何绘制 Python 中耦合非线性二阶 ODE 的二阶导数？

Question

我是 Python 的新手，编写了这段代码来模拟 spring 钟摆的运动：

import numpy as np
from scipy.integrate import odeint
from numpy import sin, cos, pi, array
import matplotlib.pyplot as plt

init = array([0,pi/18,0,0]) 

def deriv(z, t):
    x, y, dxdt, dydt = z
    dx2dt2=(4+x)*(dydt)**2-5*x+9.81*cos(y)
    dy2dt2=(-9.81*sin(y)-2*(dxdt)*(dydt))/(0.4+x)

    return np.array([dxdt, dydt, dx2dt2, dy2dt2])

time = np.linspace(0.0,10.0,1000)
sol = odeint(deriv,init,time)

plt.xlabel("time")
plt.ylabel("y")
plt.plot(time, sol)
plt.show()

但它给了我 x、dxdt、y 和 dydt 的图表，而不是 dx2dt2 和 dy2dt2（这分别是 x 和 y 的二阶导数）。如何更改我的代码以绘制二阶导数？

Answer 1

odeint 的 return 值是您定义为 z = [x,y,x',y'] 的 z(t) 的解。因此，二阶导数不是 return 由 odeint 编辑的解的一部分。您可以通过取一阶导数的 returned 值的有限差分来近似 x 和 y 的二阶导数。

例如：

import numpy as np
from scipy.integrate import odeint
from numpy import sin, cos, pi, array
import matplotlib.pyplot as plt

init = array([0,pi/18,0,0]) 

def deriv(z, t):
    x, y, dxdt, dydt = z
    dx2dt2=(4+x)*(dydt)**2-5*x+9.81*cos(y)
    dy2dt2=(-9.81*sin(y)-2*(dxdt)*(dydt))/(0.4+x)

    return np.array([dxdt, dydt, dx2dt2, dy2dt2])

time = np.linspace(0.0,10.0,1000)
sol = odeint(deriv,init,time)

x, y, xp, yp = sol.T

# compute the approximate second order derivative by computing the finite
# difference between values of the first derivatives
xpp = np.diff(xp)/np.diff(time)
ypp = np.diff(yp)/np.diff(time)

# the second order derivatives are now calculated at the midpoints of the
# initial time array, so we need to compute the midpoints to plot it
xpp_time = (time[1:] + time[:-1])/2

plt.xlabel("time")
plt.ylabel("y")
plt.plot(time, x, label='x')
plt.plot(time, y, label='y')
plt.plot(time, xp, label="x'")
plt.plot(time, yp, label="y'")
plt.plot(xpp_time, xpp, label="x''")
plt.plot(xpp_time, ypp, label="y''")
plt.legend()
plt.show()

或者，由于您已经有一个函数来计算解的二阶导数，您可以调用该函数：

plt.plot(time, deriv(sol.T,time)[2], label="x''")
plt.plot(time, deriv(sol.T,time)[3], label="y''")

如何绘制 Python 中耦合非线性二阶 ODE 的二阶导数？

How to graph the second derivatives of coupled non-linear second order ODEs in Python?

python

numpy

matplotlib

differential-equations

odeint