Как правильно реализовать scipy.integrate.Radau? - PullRequest
Новая
       84

Как правильно реализовать scipy.integrate.Radau?

0 голосов
/ 31 марта 2020

(я буду признателен даже за ссылку на пример, пока я не нашел ни одного.)

Я пытаюсь использовать Радау из scipy.integrate для решения второго дифференциальное уравнение порядка. Сейчас я пробую простой пример, чтобы понять, как он работает (пока безуспешно).

Допустим, у меня есть следующее уравнение:

d ^ 2y / dt ^ 2 = C,

, что означает, что y = (C t ^ 2) / 2 + B t.

Скажем, например, y (1) = 1 и C = 2. Допустим, я хочу найти значение y для t = 10.

Это мой код: 

from scipy.integrate import Radau
import numpy as np

C = 2.0
y = np.zeros(2)

def fun(t, y):
  #y[0] = C*t
  y[1] = C
  return y

t0 = 1.0
y0 = np.array([1., 1.])
t_bound = 10.0

eq = Radau(fun, t0, y0, t_bound)
print(eq.n)

while(True):
  print(eq.t)
  print(eq.y)
  print(eq.status)
  if eq.status == 'finished':
    break
  eq.step()

Выводы неверны , (Если я раскомментирую эту строку с комментариями в определении fun, это также даст неверный ответ. Но я думаю, что мне даже не нужно было говорить решателю, верно? Обычно я не знаю этого значения.)

Я думаю, что моя самая большая проблема в том, что я не совсем уверен, что должно быть передано как fun. Документация говорит, что это должна быть правая часть системы, поэтому я подумал, что первый вывод должен быть в y [0], второй вывод в y [1] et c.

Что я делаю неправильно ? Как это должно быть реализовано?

1 Ответ

2 голосов
/ 31 марта 2020

Атм вы решаете 

y0' = y0,  y0(1)=1
y1' = 2,   y1(1)=1

, который имеет решение y0(t)=exp(t-1) и y1(t)=2*t-1, что, конечно, не то, что вы хотите. Требуется система первого заказа 

y0' = y1
y1' = C

, чтобы вам потребовалось 

def fun(t,y): return [y[1], C]

Тогда решение будет y1(t)=C*t+B=2*t-1 и y0(t)=0.5*C*t^2+B*t+A=t^2-t+1, и интеграция завершится правильно с eq.y = [91. 19.].

...