Runge Kutta Algorithm using while loop

Question

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

`t = t + h` floating point accumulation of a fixed step: don't do that — Jean-François Fabre, Dec 14 '18 at 20:55
How do I fix this and do you think that this is the problem? — Ian Ronk, Dec 14 '18 at 20:57
it's not the problem. It's one of the problems. The most blatant one is that `(k1 + 2*k2 + 2*k3 + k4)/6` should be `h*(k1 + 2*k2 + 2*k3 + k4)/6`. — Jean-François Fabre, Dec 14 '18 at 20:58

Jean-François Fabre · Accepted Answer · 2018-12-15T08:24:44.697

3

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

and coefficients too:

(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

edited Dec 15 '18 at 08:24

answered Dec 14 '18 at 21:03

Jean-François Fabre

137,073
23
153
219

I see it now, I used the wikipedia formulas and looked over the h. Thank you very much, will set this as the answer when possible. – Ian Ronk Dec 14 '18 at 21:10
This still gives the value at `t=5.01`, one step too many. Properly ended the value should be `1.24585816213e+18` which is close enough to the exact value `np.exp(5**3/3) = 1.2462449532958323e+18`. – Lutz Lehmann Dec 15 '18 at 01:29
@LutzL thanks A LOT for this. while loops with conditions are one-off error prone, I know that... but did it anyway. Changed to a `for` loop. and got the result right :) – Jean-François Fabre Dec 15 '18 at 08:25
@IanRonk check the update. There was a one-off error so the result was too high – Jean-François Fabre Dec 15 '18 at 08:26

Runge Kutta Algorithm using while loop

1 Answers1