Python Forum

Hello all! I'm writing a simple Runge-Kutta solver and seem to be running into issues.

import numpy as np
import matplotlib.pyplot as plt


y0=1

"The amplitude function derivative, a complex exponential"
def amp(x,s):
    return .4j*x

"An explicit form of the same derivative"
def expamp(x,s):
    return .4j*np.exp(.4j*s)

"A simple Runge-Kutta solving algorithm"
def rk4(f,  t):
    y=np.zeros(len(t), dtype=complex) "I feel like the problem originates here"
    y[0]=y0 "Initial Value"
    h=t[1]-t[0] "Iteration time step"
    for i in range(0, len(t)-1):
        K1=(f(y[i], t[i]))
        K2=(f(y[i]+K1/2, t[i]+h/2))
        K3=(f(y[i]+K2/2, t[i]+h/2))
        K4=f(y[i]+K3, t[i]+h)
        y[i+1]=y[i]+h*(K1+K2+K3+K4)/6
        
    return y

"The Analytical Solution"    
def anal(t):
    return np.exp(.4j*t)

time=np.linspace(0,100,1000)
sol=rk4(amp, time)
exp=rk4(expamp, time)
an=anal(time)
plt.plot(time, an, 'b-', label='analytical')
plt.plot(time, exp, 'g.-', label='explicit RK solved')
plt.plot(time, sol, 'k--', label='RK Solved implicit')
plt.ylabel('probability')
plt.xlabel('time')
plt.legend()

Then I plotted the implicit, explicit, and analytical solutions and while the explicit and analytical solutions both show full oscillatory behavior, the implicitly solved one shows something similar to a damped harmonic oscillator(waves decaying eventually down to zero) (I can't figure out how to put the image on here, my school computer won't let me transfer images to get a URL)
I'm sure there's something in some call somewhere that's telling it to do this, but I'm having trouble finding where. Feel free to just copy/paste. I think you need pyplot and numpy the way I've written it. Can anyone shed any light on this for me?

Thank you all so much for your help!

The explicit solution is very far from the analytical solution. Why is that so?
Replacing line 25 with

        y[i+1]=y[i]+h*(K1+2*(K2+K3)+K4)/6

gives a much better explicit solution, following the formula given by Wikipedia.

In fact the other formulas are also wrong. I get good results with both methods and

        K1=(f(y[i], t[i]))
        K2=(f(y[i]+h*K1/2, t[i]+h/2))
        K3=(f(y[i]+h*K2/2, t[i]+h/2))
        K4=f(y[i]+h*K3, t[i]+h)
        y[i+1]=y[i]+h*(K1+2*(K2+K3)+K4)/6