Ошибка Python scipy.integrate.odeint для простого моделирования гравитации - PullRequest
Новая
       137

Ошибка Python scipy.integrate.odeint для простого моделирования гравитации

0 голосов
/ 04 марта 2019

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

Проблема в том, что я получаю следующее сообщение об ошибке: 

ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)

А также, что-то явно идет не так -уравнения не интегрируются правильно, а движение неверно.Ниже приведен график движения для начальных условий, которые должны давать круговое движение вокруг начала координат:

image

Это код: 

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

G=1
m=1
def f_grav(y, t):
    x1, x2, v1, v2 = y
    m = t
    dydt = [v1, v2, -x1*G*m/(x1**2+x2**2)**(3/2), -x2*G*m/(x1**2+x2**2)**(3/2)]
    return dydt

t = np.linspace(0, 100, 1001)
init = [0, 1, 1, 0]
ans = odeint(f_grav, init, t)
print(ans)

x = []
y = []
for i in range (100):
    x.append(ans[i][0])
    y.append(ans[i][1])
plt.plot(x, y)
plt.show()

Обратите внимание, что я использовал эту функцию раньше, и написание почти идентичного кода для дифференциального уравнения SHM дает правильные результаты.Изменение чисел в t не помогает.Кто-нибудь знает, почему это так плохо?

1 Ответ

0 голосов
/ 04 марта 2019

Неправильное движение, вероятно, является числовой нестабильностью, в любом случае, из документации odeint:

примечание: для нового кода используйте scipy.integrate.solve_ivp для решения дифференциального уравнения.

solve_ivp на самом деле принимает только границы и определяет количество точек, чтобы метод интегрирования был устойчивым для уравнения.Вы также можете выбрать метод интеграции. 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

G=1
m=1
def f_grav(t, y):
    x1, x2, v1, v2 = y
    m = t
    dydt = [v1, v2, -x1*G*m/(x1**2+x2**2)**(3/2), -x2*G*m/(x1**2+x2**2)**(3/2)]
    return dydt

domain = (0, 100)
init = [0, 1, 1, 0]
ans = solve_ivp(fun=f_grav, t_span=domain, y0=init)

plt.plot(ans['y'][0], ans['y'][1])
plt.show()

, при этом я не получаю никаких предупреждений, и симуляция выглядит лучше (обратите внимание, что функция должна иметь параметры в порядке (t, y)).

...