The backward Euler's method

The examples in the previous post suggested the importance of step-size ‘h’ for numerical integration, and how improper choices of ‘h’ may lead to a divergent solution. The Euler’s method presented in the former post, as indicated earlier, marched from time tn to tn+1 using the knowledge of the solution at time tn. This method is also known as the Forward Euler, since it marches forward in time.

The forward Euler at time tn uses the derivative at tn to march to tn+1. Instead, if the derivative is obtained at time tn+1, the method becomes the Backward Euler. Mathematically, it is represented by the following equation:



Since the function y(x) appears on both sides of the equation, one can no longer simply march forward between time steps. Rather, one needs to perform an iterative root-finding algorithm to determine y(x) at each time step

A simplified algorithm for the Backward-Euler can be summarized as follows:

1)    % Define step size (h), initial function y(0), initial time t0, final time tf, eps
2)    nosteps = (tf – t0)/h
3)    for i ← 1 to nosteps

- Calculate derivative function f[t0,y[i]]
- Iterate until y[i] – (y[i-1] + h*f[x(i), y(i)]) < eps
- t0 ← t0 + h

4)    % Display results

‘eps’ is a user-specified tolerance value for the root-solver. The solver iterates until the equation at each time step is less than ‘eps’. The value that satisfies this condition is the numerical solution at that time. The iteration can be performed using any root-finding algorithm like the Bisection or the Newton-Raphson method. A modified form of the subroutine ‘rtbis’ adapted from Numerical Recipes, (The art of scientific computing, Press et al, 2^nd ed.) was used for this exercise.

The Backward Euler’s method was applied to the following problem (Atkinson et al, Numerical Solution of Ordinary Differential Equations, Eg. 5.8).

 

 The solution is shown in fig. 1.

 

Euler’s method part 2- Does it work for all cases?

Let us now consider applying the Euler’s method to the problem described in Curtiss & Hirschfelder (Integration of stiff equations, Proc. of the Natl. Acad. Of Sciences, 1952, Vol. 38, pp 235-243). This is the first order DE described by the equation

 


The results of the numerical integration are shown in fig. (1) for various step sizes ‘h’, for the initial value y(0)= 1.

 

 Fig. (1): Solution of a quadratic equation for various 'h' values

Clearly the numerical solution changes for various step sizes. So which of the solution curves in fig. (1) is the “correct” solution? The simplistic marching scheme presented by the Euler’s method is clearly not “stable” for all values of step sizes. In other words, the Euler’s method does not converge to a solution for any value of ‘h’. A minimum value of ‘h’ is thus required for the numerical solutions to converge. Note on fig. (1) that the solutions tend to converge as ‘h’ decreases.


The quest for a numerical integration procedure for an ODE is thus to come up with one that is stable for any and every value of ‘h’. The forward-Euler method clearly does not appear to meet this criteria for certain types of DE’s.


To illustrate the point, the solution of the equation

 

 from Ernst Hairer & Christian Lubich, (paper titled Numerical Solution of Ordinary Differential Equations, reference unknown) for various values of ‘h’ is shown in fig. (2).

 

Fig (2): Solution for y'(x) = -50(y-cosx)

These two examples clearly illustrate the need for numerical integration techniques that can converge to a solution independent of the step size ‘h’.

The Euler's Method

Given a differential equation of the form  , a curious mind (the kind of mind that has nothing better to do in life) may wonder how one can go about solving such a DE to produce a variety of colorful numerical results. A simple closed-form analytical solution would obviously be ideal to make things simple. However, let us assume for the moment that an analytical solution either does not exist, or cannot be determined. In most real world cases, this is a perfectly valid assumption. So, how does one go about numerically solving the above differential equation?


Perhaps the simplest of numerical integration techniques is the Euler’s method. The premise of the Euler method is based on the definition of the slope of a curve at a given point from basic calculus.




 

Rearranging the above eqn., one can obtain


 

The newly rearranged equation above represents a method for marching from stage ‘x’ to ‘x+h’, given a value for y’(x) at x=0. In other words, given f(0), one can march along from y(0) with incrementing values of h to obtain a numerical solution for the 1^st order ODE. Thus,


Step 0: f(0) = initial value
Step 1: f(step=1) = f(step=0) + hf’(step=0)
Step 2: f(step=2) = f(step=1) + hf’(step=1)
...

The recursive statements are usually represented in the following format:

 

The algorithm for Euler's method of numerical integration can be summarized as follows:


1) % Define step size (h), initial function y(0), initial time t0, final time tf
2) nosteps = (tf – t0)/h
3) for i ← 1 to nosteps

- Calculate derivative function f[t0,y[i]]
- y[i] ← y[i-1] + h*f[x,y]
- t0 ← t0 + h

4) % Display results

To illustrate the procedure, the Euler’s method was applied to the following initial value problem (Butcher, Numerical methods for ordinary differential equations, 201(a)) with a step size of 0.001.

 

 The result of the integrated function is plotted in fig. (1).