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

This would probably fit better on the Computational Science SE — Tyberius, Oct 25 '20 at 22:10
`y[i,:]=y[i-1,:]+(h/6)*(k1+2*k3+2*k3+k4)` needs to be part of the loop. — Peter Meisrimel, Oct 26 '20 at 09:56

score 1 · Answer 1 · answered Oct 26 '20 at 12:42

I assume the indentation is correct in the code that you run, as the numerical solution would look completely wrong else. So my guess is that your exact solution is wrong, or the exact solution is to a different equation than the numerical code.

With the characteristic polynomial q^2+4q+2=(q+2)^2-2=0 the exact solution is

y(t)=exp(-2t)*(A*cosh(sqrt(2)*t)+B*sinh(sqrt(2)*t))

By the initial conditions, A=1 and -2*A+sqrt(2)*B=3

Correcting the indentation of the RK4 time loop and adding the exact solution gives the plot

Python evaluating a second order ODE with RK4

1 Answers1