Python evaluating a second order ODE with RK4

Question

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

Answer 1

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