Search code examples
pythonnumerical-methodsrunge-kutta

Runge–Kutta fourth-order method

I'm implementing Runge–Kutta fourth-order method for system of two equations.

enter image description here

h is number of segments, so T/h is step.

def cauchy(f1, f2, x10, x20, T, h):
    x1 = [x10]
    x2 = [x20]

    for i in range(1, h):
        k11 = f1((i-1)*T/h, x1[-1], x2[-1])
        k12 = f2((i-1)*T/h, x1[-1], x2[-1])
        k21 = f1((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k11, x2[-1] + T/h/2*k12)
        k22 = f2((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k11, x2[-1] + T/h/2*k12)
        k31 = f1((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k21, x2[-1] + T/h/2*k22)
        k32 = f2((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k21, x2[-1] + T/h/2*k22)
        k41 = f1((i-1)*T/h + T/h, x1[-1] + T/h*k31, x2[-1] + T/h*k32)
        k42 = f2((i-1)*T/h + T/h, x1[-1] + T/h*k31, x2[-1] + T/h*k32)

        x1.append(x1[-1] + T/h/6*(k11 + 2*k21 + 2*k31 + k41))
        x2.append(x2[-1] + T/h/6*(k12 + 2*k22 + 2*k32 + k42))

    return x1, x2

Then I'm testing it on this system:

enter image description here

def f1(t, x1, x2):
    return x2

def f2(t, x1, x2):
    return -x1

def true_x1(t):
    return np.cos(t) + np.sin(t)

def true_x2(t):
    return np.cos(t) - np.sin(t)

It seems to be working fine (I also tested it with different initial values and different functions: all works fine):

x10 = 1
x20 = 1
T = 1
h = 10

x1, x2 = cauchy(f1, f2, x10, x20, T, h)
t = np.linspace(0, T, h)

plt.xlabel('t')
plt.ylabel('x1')
plt.plot(t, true_x1(t), "blue", label="true_x1")
plt.plot(t, x1, "red", label="approximation_x1")
plt.legend(bbox_to_anchor=(0.97, 0.27))
plt.show()

plt.xlabel('t')
plt.ylabel('x2')
plt.plot(t, true_x2(t), "blue", label="true_x2")
plt.plot(t, x2, "red", label="approximation_x2")
plt.legend(bbox_to_anchor=(0.97, 0.97))
plt.show()

enter image description here enter image description here

Then I want to check if the error is on the order of O(step^4), so I reduce step and compute error like this:

enter image description here

step = []
x1_error = []
x2_error = []
for segm in reversed(range(10, 1000)):
    x1, x2 = cauchy(f1, f2, x10, x20, T, segm)
    t = np.linspace(0, T, segm)
    step.append(1/segm)
    x1_error.append(np.linalg.norm(x1 - true_x1(t), np.inf))
    x2_error.append(np.linalg.norm(x2 - true_x2(t), np.inf))

And I get this:

plt.plot(step, x1_error, label="x1_error")
plt.plot(step, x2_error, label="x2_error")
plt.legend()

enter image description here

So, the error is linear from step. This is really strange, because it is supposed to be on the order of O(step^4). Can anyone tell me what I did wrong?

Solution 
  • for i in range(1, h):

    This will iterate from 1 to h-1. As the last step is missing, the difference from x[h-1] at time T-T/h to the exact solution at time T is O(T/h).

    Thus use 

    for i in range(1,h+1):

    for h steps from i-1 to i or 

    for i in range(h):

    for h steps from i to i+1.

    Also, np.linspace(0,1,4) will produce 4 equally spaced numbers where the first is 0 and the last is 1, resulting in 

    array([ 0.        ,  0.33333333,  0.66666667,  1.        ])

    which is probably not what you were expecting. Thus with the above correction use

    t = np.linspace(0, T, segm+1)

    to use the same time points in both computations.

    It would be easier to follow your code if you would use the letters in their usual meaning, where h or dt is the step size and N is the number of steps. Then define before the loop h=T/N or dt=T/N to avoid the repeated use of T/N in the function calls.