#!/usr/bin/python3
# Comparing numerical errors between Euler and RK4

import numpy as np
import matplotlib.pyplot as plt

def f(x, t): # the RHS for some Cauchy problem
    return x - x**2/2
    
a, b, = 0.0, 10.0        # interval
N = 5000                 # number of points
h = (b - a) / N          # increment


tpoints = np.arange(a, b, h)

# Euler's method:
xpoints = []
x = 1.0                 
for t in tpoints:
    xpoints.append(x)
    x += h * f(x, t)
xpointsE = np.array(xpoints)

#  RK4:    
xpoints = []
x = 1.0
for t in tpoints:
    xpoints.append(x)
    k1 = h * f(x, t)
    k2 = h * f(x + k1/2, t + h/2)
    k3 = h * f(x + k2/2, t + h/2)
    k4 = h * f(x + k3, t + h)
    x += (k1 + 2*k2 + 2*k3 + k4)/6
xpointsRK4 = np.array(xpoints)

exactx = 2./(1. + np.exp(-tpoints))

plt.xlabel("t")
plt.ylabel("Numerical error with Euler")    
plt.plot(tpoints, abs(xpointsE - exactx))

plt.figure()
plt.xlabel("t")
plt.ylabel("Numerical error with Runge-Kutta 4")    
plt.plot(tpoints, abs(xpointsRK4 - exactx))

plt.show()