Как построить график вторых производных связанных нелинейных ОДУ второго порядка в Python?

Question

Я очень новичок в 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 соответственно).Как я могу изменить свой код, чтобы отобразить вторые производные вместо этого?

Answer 1

Возвращаемое значение odeint является решением z(t), которое вы определили как z = [x,y,x',y'].Следовательно, вторая производная не является частью решения, возвращаемого 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''")