pythonscipyodenumerical-integration

Why does Python RK23 Solver blow up and give unrealistic results?


I am experimenting with RK45/RK23 solver from python scipy module. Using it to solve simple ODEs and it is not giving me the correct results. When I manually code for Runge Kutta 4th order it works perfectly so does the odeint solver in the module but RK23/RK45 does not. If someone could aid me in figuring out the issue it would be help full. I have so far implemented only simple ODEs

dydt = -K*(y^2)

Code:

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()

Output: (without showing RK45 result)

enter image description here

Output: (with just RK45 plot displayed)

enter image description here

I cant find out where am I making a mistake


Solution

  • Okay I have figured out the solution. RK45 requires function definition to be like fun(t,y) and odeint requires it to be func(y,t) due to which giving them same function will result in different results.