#python #scipy #ode #integrate
#python #scipy #ода #интегрировать
Вопрос:
Я экспериментирую с решателем RK45 / RK23 из модуля python scipy. Используя его для решения простых ODE, он не дает мне правильных результатов. Когда я вручную кодирую для 4-го порядка Runge Kutta, он работает отлично, так же как и решатель odeint в модуле, но RK23 / RK45 этого не делает. Если бы кто-нибудь мог помочь мне разобраться в проблеме, это была бы полная помощь. До сих пор я реализовал только простые ODE
dydt = -K * (y ^ 2)
Код:
import numpy as np
from scipy.integrate import solve_ivp,RK45,odeint
import matplotlib.pyplot as plt
# function model fun(y,t)
def model(y,t):
k = 0.3
dydt = -k*(y**2)
return dydt
# intial condition
y0 = np.array([0.5])
y = np.zeros(100)
t = np.linspace(0,20)
t_span = [0,20]
#RK45 implementation
yr = solve_ivp(fun=model,t_span=t_span,y0=y,t_eval=t,method=RK45)
##odeint solver
yy = odeint(func=model,y0=y0,t=t)
##manual implementation
t1 = 0
h = 0.05
y = np.zeros(21)
y[0]=y0;
i=0
k=0.3
##Runge Kutta 4th order implementation
while (t1<1.):
m1 = -k*(y[i]**2)
y1 = y[i] m1*h/2
m2 = -k*(y1**2)
y2 = y1 m2*h/2
m3 = -k*(y2**2)
y3 = y2 m3*h/2
m4 = -k*(y3**2)
i=i 1
y[i] = y[i-1] (m1 2*m2 2*m3 m4)/6
t1 = t1 h
#plotting
t2 = np.linspace(0,20,num=21)
plt.plot(t2,y,'r-',label='RK4')
plt.plot(t,yy,'b--',label='odeint')
#plt.plot(t2,yr.y[0],'g:',label='RK45')
plt.xlabel('time')
plt.ylabel('y(t)')
plt.legend()
plt.show()
Вывод: (без отображения результата RK45)
введите описание изображения здесь
Вывод: (отображается только график RK45)
введите описание изображения здесь
Я не могу выяснить, где я допускаю ошибку
Ответ №1:
Хорошо, я нашел решение. RK45 требует, чтобы определение функции было похоже на fun(t, y), а odeint требует, чтобы оно было func(y, t), из-за чего предоставление им одной и той же функции приведет к разным результатам.
Комментарии:
1. К вашему СВЕДЕНИЮ: если вы хотите использовать одно и то же определение функции с обоими
odeint
иsolve_ivp
, определите его в формеfunc(t, y)
и добавьте параметрtfirst=True
вodeint
вызов функции.2. Спасибо за профессиональный совет :))