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

Question

我得到了两个二阶 ODE，我被要求用 python 中的 odeint 解决它们。

这些是等式：

d^x(t)/dt^2 = 10dy(t)/dt + x(t) - (k + 1)(x(t))/z^3

d^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 npimport scipy.integrateimport matplotlib.pyplot as pltdef 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, ddyddtinit = [  0.921, #x[0]  0,     #y[0]  0,     #x'[0]  3.0    #y'[0]]k = 1times  = np.linspace(0,1,1000)values = scipy.integrate.odeint(function, init, times, args=(k,), tfirst=False)plt.plot(values)plt.show()并给出这个输出：