# Python evaluating a second order ODE with RK4

Pasted below is my python code. It is a 4th order runge kutta that evaluates the 2nd order ode: y'' +4y'+2y=0 with initial conditions y(0)=1, y'(0)=3.

I need help fixing it. When I run my code, my analytical solution does not match my numerical solution, my professor said they should be the same. I have tried editing this a bunch and cannot seem to figure out what's wrong. If anyone could review my code and let me know if there is something wrong that would be great. Thank you!

import numpy as np
import matplotlib.pyplot as plt

def ode(y):
return np.array([y[1],(-2*y[0]-4*y[1])])

tStart=0

tEnd=5

h=.1

t=np.arange(tStart,tEnd+h,h)

y=np.zeros((len(t),2))

y[0,:]=[1,3]

for i in range(1,len(t)):
k1=ode(y[i-1,:])
k2=ode(y[i-1,:]+h*k1/2)
k3=ode(y[i-1,:]+h*k2/2)
k4=ode(y[i-1,:]+h*k3)

y[i,:]=y[i-1,:]+(h/6)*(k1+2*k3+2*k3+k4)

plt.plot(t,y[:,0])
plt.plot(t,1-t)
plt.grid()
plt.gca().legend(('y(t)',"y'(t)"))
plt.show()
$$$$

• Your update is not in the for-loop body
– VoB
Commented Oct 26, 2020 at 8:50
• One suggestion: to debug a relatively complex algorithm like RK4 it is usually a good idea to simplify the algorithm to something very basic, like RK1, and see if that works. Commented Oct 27, 2020 at 5:38
• Cross-post to stackoverflow.com/questions/64529228/… with my answer guessing that the exact solution that is compared against is not for the given IVP. A freshly computed exact solution matches the numerical solution (with corrected indentation) perfectly. Commented Oct 27, 2020 at 20:20

Well, your problem is the numerical solution to an initial value problem. First, we must transform the second-order ODE into a 2x2 system of first-order ODE. Note that despite the transformation, the problem remains well posed, as we have two ODE with two initial conditions. First, I will use the method developed in scipy, scipy.integrate.odeint.

# Author: Carlos eduardo da Silva Lima
# Solving EDO initial value problem (IVP) via scipy and 4Order Runge-Kutta
# scipy

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

# Initial conditions
t_initial = 0.0
t_final = 10.0
y0  = 1.0
u0  = 3.0
tol = 1e-8
N   = 10000

# Enter the definition of the set of ordinary differential equations
def ode(s,t):
y = s[0]; u = s[1]
ode_1 = u
ode_2 = -2*y-4*u
return np.array([ode_1,ode_2])

# Resolution of the initial value problem (IVP) via scipy.inetgrate.odeint
t = np.linspace(t_initial,t_final,N)
s0 = np.array([y0,u0])
sol = odeint(ode,s0,t,rtol=tol)

# Result applied to y and u arrays
y = sol[:,0]
u = sol[:,1]

# Graphics
plt.style.use('dark_background')
plt.figure(figsize=(7,7))
plt.xlabel(r'$$t(s)$$')
plt.ylabel(r'$$y(t)$$ and $$u(t)$$')
plt.title(r'$$\frac{d^{2}y(x)}{dt^{2}}+ + 4\frac{dy(x)}{dt} + 2y(x) = 0$$ with $$y(t_{0} = 0) = 1$$ and $$\frac{dy(0)}{dt} = 3$$')
plt.plot(t,y,'b-o',t,u,'r-o')
plt.grid()


Graphic (Odeint)

Now we will apply the 4th Order Runge-Kutta algorithm.

# Author: Carlos eduardo da Silva Lima
# Solving EDO initial value problem (IVP) via scipy and 4Order Runge-Kutta
# 4Order Runge-Kutta

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

# Initial conditions
t_initial = 0.0
t_final = 50.0
y0  = 1.0
u0  = 3.0
N   = 10000
h   = 1e-3 # Step

# Enter the definition of the set of ordinary differential equations
def ode(t,y,u):
ode_1 = u
ode_2 = -2*y-4*u
return np.array([ode_1,ode_2])

# RK4
t = np.empty(N)
y = np.empty(N); u = np.empty(N)

t[0] = t_initial
y[0] = y0; u[0] = u0

for i in range(0,N-1,1):

k11 = h*ode(t[i],y[i],u[i])[0]
k12 = h*ode(t[i],y[i],u[i])[1]

k21 = h*ode(t[i]+(h/2),y[i]+(k11/2),u[i]+(k12/2))[0]
k22 = h*ode(t[i]+(h/2),y[i]+(k11/2),u[i]+(k12/2))[1]

k31 = h*ode(t[i]+(h/2),y[i]+(k21/2),u[i]+(k22/2))[0]
k32 = h*ode(t[i]+(h/2),y[i]+(k21/2),u[i]+(k22/2))[1]

k41 = h*ode(t[i]+h,y[i]+k31,u[i]+k32)[0]
k42 = h*ode(t[i]+h,y[i]+k31,u[i]+k32)[1]

y[i+1] = y[i] + ((k11+2*k21+2*k31+k41)/6)
u[i+1] = u[i] + ((k12+2*k22+2*k32+k42)/6)
t[i+1] = t[i] + h

# Graphics
plt.style.use('dark_background')
plt.figure(figsize=(7,7))
plt.xlabel(r'$$t(s)$$')
plt.ylabel(r'$$y(t)$$ and $$u(t)$$')
plt.title(r'$$\frac{d^{2}y(x)}{dt^{2}}+ + 4\frac{dy(x)}{dt} + 2y(x) = 0$$ with $$y(t_{0} = 0) = 1$$ and $$\frac{dy(0)}{dt} = 3$$')
plt.plot(t,y,'b-o',t,u,'r-o')
plt.grid()
`

Graphic (Runge-Kutta 4Order)

Some results and comparisons between odeint and RK4

• This answer is not substantially better than the question. You have also an indentation error in the first code, the OP code already contains the correct transformation to a first-order system, and it uses the better vector implementation of RK4 that is less error-prone than the component-wise implementation. So the information gained is that, in addition to the exact solution from the SO answer also the numerical solution with an established method confirms the result of the OP implementation. Commented Dec 19, 2022 at 11:51