Использование решения scipy solve_ivp для решения нелинейного движения маятника - PullRequest
Новая
       35

Использование решения scipy solve_ivp для решения нелинейного движения маятника

1 голос
/ 15 мая 2019

Я все еще пытаюсь понять, как решает execute_ivp против odeint, но как только я понял, что-то произошло.

Я пытаюсь найти движение нелинейного маятника. С odeint все работает как шарм, но на solve_ivp происходит что-то странное: 

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

g = 9.81
l = 0.1


def f(t, r):
    omega = r[0]
    theta = r[1]
    return np.array([-g / l * np.sin(theta), omega])


time = np.linspace(0, 10, 1000)
init_r = [0, np.radians(179)]

results = solve_ivp(f, (0, 10), init_r, method="RK45", t_eval=time) #??????
cenas = odeint(f, init_r, time, tfirst=True)


fig = plt.figure()
ax1 = fig.add_subplot(111)

ax1.plot(results.t, results.y[1])
ax1.plot(time, cenas[:, 1])

plt.show()

enter image description here

Чего мне не хватает?

1 Ответ

1 голос
/ 15 мая 2019

Это численная проблема. Относительные и абсолютные допуски по умолчанию solve_ivp составляют 1e-3 и 1e-6 соответственно. Для многих проблем эти значения слишком низкие. Относительный допуск по умолчанию для odeint составляет 1,49e-8.

Если вы добавите аргумент rtol=1e-8 к вызову solve_ivp, графики согласуются: 

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

g = 9.81
l = 0.1


def f(t, r):
    omega = r[0]
    theta = r[1]
    return np.array([-g / l * np.sin(theta), omega])


time = np.linspace(0, 10, 1000)
init_r = [0, np.radians(179)]

results = solve_ivp(f, (0, 10), init_r, method='RK45', t_eval=time, rtol=1e-8)
cenas = odeint(f, init_r, time, tfirst=True)


fig = plt.figure()
ax1 = fig.add_subplot(111)

ax1.plot(results.t, results.y[1])
ax1.plot(time, cenas[:, 1])

plt.show()

Участок:

plot

...