python中如何解决二阶常微分方程？在两个二阶微分中的每一个中都有两个变量？

Question

我得到了两个二阶 ODE，我被要求用 python 中的 odeint 解决它们。这些是等式：d^x(t)/dt^2 = 10dy(t)/dt + x(t) - (k + 1)(x(t))/z^3d^2y(t)/dt^2 = - 10dy(t)/dt + y(t) - ((k+1)(y(t) + k))/z^3在哪里 z = np.sqrt((y+k)^2+x^2))我已经获得了初始变量 (x, y, dxdt, dydt) 我知道它们的值，但我不会坚持输入它们，所以我不会把它们放在这里。def function(init, time, k):    xt = init[0]    yt = init[1]    z = np.sqrt((init[1]+k)^2+init[0]^2))    dxdt = init[2]    dydt = init[3]    ddxddt = 10*dydt + xt - ((k+1)(xt))/z^3    ddyddt = -10*dxdt + xt - ((k+1)(yt + k))/z^3    return(dxdt, ddxddt, dydt, ddyddt)init = [0.921, 0, 0, 3.0]values = odeint(function, initial, time, args(k,))在此之后，我定义了初始值，并定义了时间 k，并将它们放入 odeint。但是我可以看到我的实际设置功能确实有问题。我不明白如何拆分二阶颂歌。

Answer 1

你这里有几个错误。

第一：z^3不是幂，而是异或运算。在 Python 中，功能是使用**运算符完成的，因此您需要编写z**3.

第二：您错误地命名了函数的参数。代替：

def function(init, time, k):

你应该有

def function(state, time, k):

因为state根据函数返回的导数进化。它只会在第一个时间步中具有初始值。

第三：您的状态解释和状态增量不一致。你写：

xt   = init[0]

yt   = init[1]

dxdt = init[2]

dydt = init[3]

但后来

return dxdt, ddxddt, dydt, ddyddt

这意味着，除其他外，dydt=ddxddt. 你应该写：

xt, yt, dxdt, dydt = state

[....]

return dxdt, dydt, ddxddt, ddyddt

请注意，您必须确保您的初始条件与您订购状态的方式一致。

正确实现的最小工作示例可能如下所示：

import numpy as np

import scipy.integrate

import matplotlib.pyplot as plt

def function(state, time, k):

  xt,yt,dxdt,dydt = state

  z               = np.sqrt((yt+k)**2+xt**2)

  ddxddt          = 10*dxdt + xt - ((k+1)*(xt    ))/z**3

  ddyddt          = -10*dydt + yt - ((k+1)*(yt + k))/z**3

  return dxdt, dydt, ddxddt, ddyddt

init = [

  0.921, #x[0]

  0,     #y[0]

  0,     #x'[0]

  3.0    #y'[0]

]

k = 1

times  = np.linspace(0,1,1000)

values = scipy.integrate.odeint(function, init, times, args=(k,), tfirst=False)

plt.plot(values)

plt.show()

并给出这个输出：