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

【问题标题】：How to graph the second derivatives of coupled non-linear second order ODEs in Python?如何在 Python 中绘制耦合非线性二阶 ODE 的二阶导数？
【发布时间】：2019-05-15 03:04:50
【问题描述】：

我对 Python 非常陌生，并且编写了这段代码来模拟弹簧钟的运动：

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 的二阶导数） .如何更改我的代码以绘制二阶导数？

【问题讨论】：

您能否包含您想要求解的原始微分方程以及如何将其转化为一阶微分方程组？我怀疑答案很简单，要获得 d2xdt2 你想调用plot((time[1:] + time[:-1])/2,np.diff(sol[:,1])/np.diff(time)) 和plot((time[1:] + time[:-1])/2,np.diff(sol[:,3])/np.diff(time))。
@user545424 我想求解一个耦合 ODE 系统。原始方程是 x'' = (0.18+x)*(y')^2-51x+9.81*cos(y) 和 (0.18+x)y''+2x'y'=-10*sin(y)，其中 x 和 y 都是关于时间的。

标签： python numpy matplotlib differential-equations odeint

【解决方案1】：

odeint 的返回值是您定义为z = [x,y,x',y'] 的z(t) 的解。因此二阶导数不是odeint 返回的解的一部分。您可以通过对一阶导数的返回值进行有限差分来近似 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''")

【讨论】：