How to compare Euler's method to a second-Order Runge-Kutta Method at the same stepsize?

Question

I have two algorithms for a numerical differential equation problem, one called Euler's method and one called a second-order Runge Kutta(RK2) . Essentially Euler's method and RK2 approximate a solution to differential equations. The only difference is that they use different formulas (Euler's uses a first order derivative of Taylor's series whilst the RK2 is a second order derivative of a Taylor series).

I am trying to correct some code that I have written in order for it to return the following solution,

However my code is not returning the correct values when it comes to the RK2 but returns the following,

Note that in my solution I refer to the step size h as dt. I have provided the code that I used to create this below followed by a numeric example of the second-order Runge Kutta method that is working numerically. I am interested in showing that the convergence is quicker with the RK2 than Euler's method.

import matplotlib.pyplot as plt
import numpy as np
from math import exp  # exponential function

dy = lambda x, y: x * y
f = lambda x: exp(x ** 2 / 2)  # analytical solution function
x_final = 2

# analytical solution
x_a = np.arange(0, x_final, 0.01)
y_a = np.zeros(len(x_a))
for i in range(len(x_a)):
    y_a[i] = f(x_a[i])
plt.plot(x_a, y_a, label="analytical")

# Container for step sizes dt /dt
dt = 0.5

x = 0
y = 1
print("dt = " + str(dt))
print("x \t\t y (Euler) \t y (analytical)")
print("%f \t %f \t %f" % (x, y, f(x)))

n = int((x_final - x) / dt)

x_n = np.zeros(n + 1)
y_n = np.zeros(n + 1)
x_n[0] = x
y_n[0] = y

#Plot Euler's method
for i in range(n):
    y += dy(x, y) * dt
    x += dt
    print("%f \t %f \t %f" % (x, y, f(x)))
    x_n[i + 1] = x
    y_n[i + 1] = y

plt.plot(x_n, y_n, "x-", label="Euler dt=" + str(dt))

###################33
# Runge-Kutta's method 2'nd order (RK2)
x = 0
y = 1
print("dt = " + str(dt))
print("x \t\t y (rk2) \t y (analytical)")
print("%f \t %f \t %f" % (x, y, f(x)))

n = int((x_final - x) / dt)

x_n = np.zeros(n + 1)
y_n = np.zeros(n + 1)
x_n[0] = x
y_n[0] = y

# Plot the RK2
for z in range(n):
    K1 = dt*dy(x,y) # Step 1
    K2 = dt*dy(x+dt/2,y+K1/K2) # Step 2
    y += K2 # Step 3
    x += dt
    print("%f \t %f \t %f" % (x, y, f(x)))
    x_n[i + 1] = x
    y_n[i + 1] = y

plt.plot(x_n, y_n, "x-", label="RK2 dt=" + str(dt))

plt.title("Euler's method ")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

This is the numeric code that works for the RK2. I have put step 1,2,and 3 on the changes from the the Euler method that allows us to create the RK2.

from math import exp
dy = lambda x,y: x*y
f = lambda x: exp(x**2/2)
x = 0
xn = 2
y = 1
dt = 0.5
n = int((xn)/dt)
print ('x \t\t y (RK2) \t y (analytical)')
print ('%f \t %f \t %f'% (x,y,f(x)))
# main loop
for i in range(n):
    K1 = dt*dy(x, y) # step 1
    K2 = dt*dy(x + dt/2, y + K1/2) # step 2
    y += K2  # step 3
    x += dt
    print ('%f \t %f \t %f'% (x,y,f(x)))

Sorry If I have asked this question badly. I am new to python, so I am still trying to understand how to go about solving these types of problems. The summary of my question is how can I fix the above function in order to calculate the correct estimate and then plot it on the graph using matplotlib from python.

Answer 1

There are two small details:

1. In the RK-2 loop you're using z, but to store the values you use i
2. In the initial code, you were using y+K1/K2 when you update K2 and that's wrong. I see you fix it in the second code.

So, the fixed code is:

import matplotlib.pyplot as plt
import numpy as np
from math import exp  # exponential function

dy = lambda x, y: x * y
f = lambda x: exp(x ** 2 / 2)  # analytical solution function
x_final = 2

# analytical solution
x_a = np.arange(0, x_final, 0.01)
y_a = np.zeros(len(x_a))
for i in range(len(x_a)):
    y_a[i] = f(x_a[i])

plt.plot(x_a, y_a, label="analytical")

# Container for step sizes dt /dt
dt = 0.5

x = 0
y = 1
print("dt = " + str(dt))
print("x \t\t y (Euler) \t y (analytical)")
print("%f \t %f \t %f" % (x, y, f(x)))

n = int((x_final - x) / dt)

x_e = np.zeros(n + 1)
y_e = np.zeros(n + 1)
x_e[0] = x
y_e[0] = y

#Plot Euler's method
for i in range(n):
    y += dy(x, y) * dt
    x += dt
    print("%f \t %f \t %f" % (x, y, f(x)))
    x_e[i + 1] = x
    y_e[i + 1] = y

plt.plot(x_e, y_e, "x-", label="Euler dt=" + str(dt))

###################33
# Runge-Kutta's method 2'nd order (RK2)
x = 0
y = 1
print("dt = " + str(dt))
print("x \t\t y (rk2) \t y (analytical)")
print("%f \t %f \t %f" % (x, y, f(x)))

n = int((x_final - x) / dt)

x_r = np.zeros(n + 1)
y_r = np.zeros(n + 1)
x_r[0] = x
y_r[0] = y

# Plot the RK2
for i in range(n):
    K1 = dt*dy(x,y) # Step 1
    K2 = dt*dy(x+dt/2,y+K1/2) # Step 2
    y += K2 # Step 3
    x += dt
    print("%f \t %f \t %f" % (x, y, f(x)))
    x_r[i + 1] = x
    y_r[i + 1] = y

plt.plot(x_r, y_r, "s-", label="RK2 dt=" + str(dt))

plt.title("numerical differential equation problem")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()