The distinct approximatons in the weighted average for my rk4 implementation are all the same

Question

Okay so I'm trying to write an implementation of the 4th order runge kutta method for the numerical approximation of a differential equation just as part of my maths course but also to learn some programming but the problem is that it uses these coefficients inside each step of the approximation. So we start out with an x value and y value and it wants to find the y value for a later x value and i give it a step size, usually 0.1, and it moves up by 0.1 in x and gives me a new approximation for the y value, and at each of these steps it's performing 4 approximations which I call k1,k2,k3 and k4. So then it takes a weighted average of these 4 approximations and gives me the final approximation which should be quite accurate. The problem is that for each step I'm getting k1=k2=k3, but k4 is different. This doesn't seem right. I'm testing it on a function, feeding it f(x,y)=2x-3y+1, using a step size h, and my k's are as follows:

k1=f(x,y)
k2=f(x+.5h,y+.5hk1)
k3=f(x+.5h,y+.5hk2)
k4=f(x+h,y+hk3)

So just examining this algebraically to have k1=k2 I end up with 9y=1-6x, but the initial values I'm feeding it are (x0,y0)=(1,5),and 45 obviously !=1-6

So this all seems wrong. I'm not getting the right answer according to the textbook, and these k's really shouldn't be the same. So anyway here's my code. The print statements in between the k's are just working as a test to see that the k's are actually changing as we cycle through n.

import numpy

def rk4(x0,y0,xf,h,f):
    y=[]
    x=numpy.linspace(x0,xf,(xf-x0)/h+1)
    y.insert(0,y0)
    for n in range(len(x)-1):
        k1=f(x[n],y[n])
        print k1
        k2=f(x[n]+(1/2)*h,y[n]+(1/2)*h*k1)
        print k2
        k3=f(x[n]+(1/2)*h,y[n]+(1/2)*h*k2)
        print k3
        k4=f(x[n]+h,y[n]+h*k3)
        print k4
        y.insert(n+1,y[n]+(h/6)*(k1+2*k2+2*k3+k4))
        print y[n]+(h/6)*(k1+2*k2+2*k3+k4)

    print x
    print y

def twoxminusthreeyplusone(x,y):
    return 2*x-3*y+1

rk4(1,5,1.5,0.1,twoxminusthreeyplusone)

So then I get this output

/usr/bin/python2.7 /home/t/PycharmProjects/untitled/chunk.py
-12.0
-12.0
-12.0
-8.2
3.86333333333
-8.39
-8.39
-8.39
-5.673
3.06961666667
-5.80885
-5.80885
-5.80885
-3.866195
2.52110925
-3.96332775
-3.96332775
-3.96332775
-2.574329425
2.14792644708
-2.64377934125
-2.64377934125
-2.64377934125
-1.65064553888
1.900100743
[ 1.   1.1  1.2  1.3  1.4  1.5]
[5, 3.8633333333333333, 3.0696166666666667, 2.5211092500000003, 2.1479264470833335, 1.9001007429979169]

Process finished with exit code 0

So I don't get it. The k's being the same seems like the main problem. I've tried it a few different ways, like before instead of just having the k's defined as functions of x[n],y[n] I had each k start as an empty list like y and populate it after each step of the for loop cycling through n using .insert(n,x) or maybe .insert(n-1,x) depending on how I defined each step but that seemed unnecessarily complicated and I think actually ended up with the same problem

Hi Tobias, welcome to SO! You provide to much information unrelated to python. You should consider, that a) not everyone is a mathematician, or b) interested in an underlying mathematical problem. Therefore you should simplify your example and explanation. — Don Question, Oct 14 '15 at 22:04
@DonQuestion Thanks, sorry. I don't really know how to explain the problem without the extraneous maths stuff, but thankfully it's solved now anyway. — Tobias Nash, Oct 14 '15 at 22:28

Lutz Lehmann · Accepted Answer · 2015-10-14T22:23:13.497

0

You have a typical problem at an usually unexpected point.

Replace all your factors (1/2) by 0.5 (or (1/2)*h by h/2, but avoid division wherever possible). The integer division results in 0 and thus in the observed behavior.

But note that with python 3 your code would work, in turn, the integer behavior has to be forced.

edited Oct 14 '15 at 22:23

answered Oct 14 '15 at 22:17

Lutz Lehmann

25,219
2
22
51

YERRRRRRSSS I now have a working rk4 implementation. Thanks dude. Typical maths guy avoiding decimals, I guess. – Tobias Nash Oct 14 '15 at 22:22
But in your short rk4 description, you used floating point just fine? Note that 0.5 as a power of 2 is exact in (binary) floating point. – Lutz Lehmann Oct 14 '15 at 22:24
Yeah I know, it's just an instinct without any actual justification, and in the description I think I used it to avoid ambiguity or having to write an extra set of brackets. – Tobias Nash Oct 14 '15 at 22:26

The distinct approximatons in the weighted average for my rk4 implementation are all the same

1 Answers1