如何在 python 中求解二阶 ODE？两个二阶微分中的每一个都有两个变量？

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根据函数returns的导数演化。它只会在第一个时间步具有初始值。

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

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()

并给出这个输出：

如何在 python 中求解二阶 ODE？两个二阶微分中的每一个都有两个变量？

How are 2nd order ODEs solved in python? With two variables in each of two second order differentials?

python

algorithm

numerical-methods

differential-equations