在Python中绘制和求解三个相关的ODE

Question

我有一个问题可能更数学化，但关键是我想用 python 解决它并绘制它。我有三个 ODE，它们通过以下方式与一个另一个相关：

x''(t)=b*x'(t)+c*y'(t)+d*z'(t)+e*z(t)+f*y(t)+g*x(t)
y''(t)=q*x'(t)+h*y'(t)+i*z'(t)+p*z(t)+l*y(t)+m*x(t)
z''(t)=a*x'(t)+w*y'(t)+v*z'(t)+u*z(t)+o*y(t)+n*x(t)

我将如何解决它们，以便通过它们的加速度将它们绘制在 3D 图表中。 我知道，我必须解决它们，困难在于它们是二阶 ODE，并且通过对一个个体的依赖而相互关联。

对于一些额外的信息，这里是代码（它实际上并没有那么好，请随意以不同的方式尝试）：

from sympy import symbols, Function, Eq
import numpy as np
from scipy.integrate import solve_ivp

t = symbols('t')
x = Function('x')(t)
y = Function('y')(t)
z = Function('z')(t)

b, c, d, e, f, g = symbols('b c d e f g')
q, h, i, p, l, m = symbols('q h i p l m')
a, w, v, u, o, n = symbols('a w v u o n')

eqx = b*x.diff(t) + c*y.diff(t) + d*z.diff(t) + e*z + f*y + g*x
eqy = q*x.diff(t) + h*y.diff(t) + i*z.diff(t) + p*z + l*y + m*x
eqz = a*x.diff(t) + w*y.diff(t) + v*z.diff(t) + u*z + o*y + n*x

#First I replaced the derivatives to turn it into a first order ODE
eqx = eqx.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})
eqy = eqy.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})
eqz = eqz.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})

#to solve them nummerically I started to lambdify them:
freex = eqx.free_symbols
freey = eqy.free_symbols
freez = eqz.free_symbols

Xl = sp.lambdify(freex , eqx , 'numpy')
Yl = sp.lambdify(freey , eqy , 'numpy')
Zl = sp.lambdify(freez , eqz , 'numpy')

#from here on I tried to solve them, but had trouble with the dependencies and the arguments
#(so here is only the line for x)
Xo = [0]
ExampleValues = np.array([0, 0, 0, 1, 9, 0, 1, 2, 4])
space= np.linspace(0, 10, 50)
Solx = solve_ivp(XSL, (0, 10), Xo, t_eval=sace, args=ExampleValues)

感谢您的提前答复！

Answer 1

由于

XSL

未定义，我无法确定问题是什么。这是代码的工作版本（带有描述更改的注释）。

from sympy import symbols, Function, Eq
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import sympy as sp  # can't forget the imports

t = symbols('t')
x = Function('x')(t)
y = Function('y')(t)
z = Function('z')(t)
# these symbols needed to be defined before use
dx, dy, dz = symbols('dx, dy, dz')

b, c, d, e, f, g = symbols('b c d e f g')
q, h, i, p, l, m = symbols('q h i p l m')
a, w, v, u, o, n = symbols('a w v u o n')

eqx = b*x.diff(t) + c*y.diff(t) + d*z.diff(t) + e*z + f*y + g*x
eqy = q*x.diff(t) + h*y.diff(t) + i*z.diff(t) + p*z + l*y + m*x
eqz = a*x.diff(t) + w*y.diff(t) + v*z.diff(t) + u*z + o*y + n*x

eqx = eqx.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})
eqy = eqy.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})
eqz = eqz.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})

# Here is where the main changes start
# To keep things simple, pass all variables to all lambda functions
# The first six are the state variables
# The remainder are parameters
free = [x, y, z, dx, dy, dz, b, c, d, e, f, g, q, h, i, p, l, m, a, w, v, u, o, n]

# You might have forgotten these equations. They relate 
# the derivatives of the first three state variables, `x`, `y`, `z`)
# to the derivative variables `dx`, `dy`, `dz`
dXl = sp.lambdify(free, dx , 'numpy')
dYl = sp.lambdify(free, dy , 'numpy')
dZl = sp.lambdify(free, dz , 'numpy')

# I've prepended `dd` to the names to indicate that these are the
# expressions for the second time derivatives
ddXl = sp.lambdify(free, eqx , 'numpy')
ddYl = sp.lambdify(free, eqy , 'numpy')
ddZl = sp.lambdify(free, eqz , 'numpy')

# The "right hand side" function
# Accepts the state variables passed by `solve_ivp` as a vector
# Accepts the parameters individidually
# Evaluates and returns the derivatives of the state variables
# as defined above
def df(t, state, *params):
    return [dXl(*state, *params), dYl(*state, *params), dZl(*state, *params),
            ddXl(*state, *params), ddYl(*state, *params), ddZl(*state, *params)]

# Here I've set values to zero so that we have
# three independent spring-mass-damper systems,
# but you can set them as you wish to couple the
# systems.
c = d = q = i = a = w = 0
f = e = m = p = n = o = 0

# Negative to get oscillatory solution
b = h = v = -0.1 
g = l = u = -1

params = b, c, d, e, f, g, q, h, i, p, l, m, a, w, v, u, o, n
x = y = z = 1
dx = dy = dz = 0
Xo = [x, y, z, dx, dy, dz]

t_eval = np.linspace(0, 10, 50)
sol = solve_ivp(df, (0, 10), Xo, t_eval=t_eval, args=params)

t = sol.t
x, y, z, dx, dy, dz = sol.y
plt.plot(t, y)