Backwards Euler Method’s error converges to step size (Python)

computational mathematicseuler's methodordinary differential equationspython

I am trying to programmatically solve the ODE,
$$\displaystyle Y'( t) =-Y( t) +\frac{1}{1+t^{2}} +\tan^{-1}( t)$$
with initial condition $Y(0)=0$. I know the analytical solution is $\displaystyle Y( t) =\tan^{-1}( t)$.

I am using the formula for the backwards Euler Method, $$\displaystyle \boxed{y_{n+1} =y_{n} +hf( t_{n+1} ,\ y_{n+1})}$$
and am using an initial guess value as instructed from a book
$\displaystyle \overset{( 0)}{y_{n+1}} =y_{n+1,\ FE}$ which I obtain from the forward Euler Method, and then iterate over $$\displaystyle \overset{( k)}{y_{n+1}} =y_{n} +hf( t_{n+1} ,\ \overset{( k-1)}{y_{n+1})}$$
until it is sufficient enough to stop.

I have used code in python (found at the end), and have plotted the results obtained using both the forwards and backwards Euler Method for reference, as well as the difference of each method from the 'true result'.
Forwards and Backwards Euler Method Plots

My main question is regarding the difference of the backwards Euler Method's result from the true result obtained from np.arctan(t), and why it converges to the step size $h$.

The code,

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

h = 0.01 # step size
t0 = 0; tf = 10
numel = int(np.floor(((tf-t0)/h) + 1))

t = np.linspace(t0, tf, numel)
yf = np.zeros(numel) # forwards array
ybold = np.zeros(numel) # backwards arrays
ybnew = np.zeros(numel)
ytrue = np.zeros(numel)

yf[0] = 0; ybold[0] = 0; ybnew[0] = 0; # initial condition

for i in range(0, numel-1):
  # Forward Method
  yf[i+1] = yf[i] + h*(-yf[i] + 1/(1+(t[i]**2)) + np.arctan(t[i]))

  # Backward Method
  ybold[i+1] = yf[i+1]
  c = 0
  diff = 1
  while  diff > 0.00001 and c < 5:
    ybnew[i+1] = ybnew[i] + h*(-ybold[i+1] + 1/(1+(t[i+1]**2)) + np.arctan(t[i+1]))
    diff = ybnew[i+1] - ybold[i+1]
    ybold[i+1] = ybnew[i+1]
    c = c + 1

ytrue = np.arctan(t)

I believe I could have made a mistake when coding it, although if I have my mistake eludes me. Or, I am misunderstanding the errors of the Euler Methods.

Any help would be appreciated, thank you.

Best Answer

  while  diff > 0.00001 and c < 5:

will break the loop when diff = -1, so change to absolute values for the error bound

  while  abs(diff) > 0.00001 and c < 5:

Note that for a first-order method an error proportional to the step size falls in the expected behavior.

Related Question