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