Как решаются ODE 2-го порядка в python? С двумя переменными в каждой из двух дифференциалов второго порядка? - PullRequest
Новая
       97

Как решаются ODE 2-го порядка в python? С двумя переменными в каждой из двух дифференциалов второго порядка?

0 голосов
/ 08 января 2019

Мне дали два ODE второго порядка, и меня попросили решить их с помощью odeint в python.

Это уравнения: 

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 и помещаю их в один момент.

Но я вижу, что я что-то действительно неправильно делаю с моей действующей функцией настройки. Я не понимаю, как разделить оды второго порядка.

1 Ответ

0 голосов
/ 08 января 2019

У вас есть несколько ошибок здесь.

Первое: 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()

и выдает следующее:

Output of odeint

...