File:Forward and backward Euler method.png

import numpy as np
import matplotlib.pyplot as plt

a = float(5)   # Maximal value of the domain
n = float(30)  # Number of divisions
dt = a/n       # Size of each subdivision

t1 = np.arange(0, a, 0.01)  # For plotting the analytical solution
t2 = np.arange(0, a, dt)    # For plotting the approximations

# The following function solves y'(t) = -y(t)^2 and y(0) = 1
def y(t):
    return 1/(1+t)

# Our approximations will simply be given as an array of approximate values
forward_euler = [1]   # The initial condition
backward_euler = [1]  # is the same for both approximations

for k in range(len(t2)-1):
    forward_euler.append(forward_euler[-1] - dt*forward_euler[-1]**2)
    backward_euler.append((-1 + np.sqrt(1+4*dt*backward_euler[-1])) / (2*dt))

# Plotting time
plt.plot(t1, y(t1), 'r-', label='Analytical solution')
plt.plot(t2, forward_euler, 'bs', label='Forward Euler method')
plt.plot(t2, backward_euler, 'g^', label='Backward Euler method')
plt.legend()
plt.show()