为什么 Python RK23 求解器会爆炸并给出不切实际的结果?
Why does Python RK23 Solver blow up and give unrealistic results?
我正在试验来自 python scipy 模块的 RK45/RK23 求解器。用它来求解简单的 ODE,但没有给出正确的结果。当我为 Runge Kutta 4 阶手动编码时,它工作得很好,模块中的 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 结果)
enter image description here
输出:(仅显示 RK45 图)
enter image description here
我找不到我哪里出错了
好的,我已经找到解决办法了。 RK45 要求函数定义类似于 fun(t,y),而 odeint 要求它类似于 func(y,t),因为给它们相同的函数会导致不同的结果。
我正在试验来自 python scipy 模块的 RK45/RK23 求解器。用它来求解简单的 ODE,但没有给出正确的结果。当我为 Runge Kutta 4 阶手动编码时,它工作得很好,模块中的 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 结果)
enter image description here
输出:(仅显示 RK45 图)
enter image description here
我找不到我哪里出错了
好的,我已经找到解决办法了。 RK45 要求函数定义类似于 fun(t,y),而 odeint 要求它类似于 func(y,t),因为给它们相同的函数会导致不同的结果。