Search code examples
pythonnumerical-methodsrunge-kutta

Runge Kutta Algorithm using while loop


I am currently trying to get the Runge Kutta 4 Integration to work, but it outputs the following:

Runge Kutta value: inf

While it is supposed to give a value in the range of:

Runge Kutta value: 8.476271005220534e+16

I used the following code, but can't seem to output the right approximation.

endtime = 5
h = 0.01

def f(x, z):
    return x*x*z


t = 0
y = 1

while t < endtime:
    k1 = f(t, y)
    k2 = f(t+(h/2), y+(k1/2))
    k3 = f(t+(h/2), y+(k2/2))
    k4 = f(t+h, y+k3)

    y = y + h*(k1 + 2*k2 + 2*k3 + k4)/6

    t = t + h

print("Runge Kutta value: " + str(y))

Is there anyone who knows where I made the mistake


Solution

  • your formulas are wrong because most of the time you're forgetting to multiply by h

    enter image description here

    and coefficients too:

    enter image description here

    (source: Wikipedia)

    so:

    endtime = 5
    h = 0.01
    
    def f(x, z):
        return x*x*z
    
    y = 1
    count = 0
    t = 0
    
    for count in range(int(endtime/h)):
        t = h*count
    
        k1 = f(t, y)
        k2 = f(t+(h/2), y+h*(k1/2))
        k3 = f(t+(h/2), y+h*(k2/2))
        k4 = f(t+h, y+h*k3)
    
        y += h*(k1 + 2*k2 + 2*k3 + k4)/6
    
    print("Runge Kutta value: " + str(y))
    

    Also avoid floating point accumulation by computing the value of t each time, instead of adding the fixed step, and trade the while loop for a for loop.

    not sure it's perfect but I get a finite result now :)

     Runge Kutta value: 1.245858162131811e+18