Dr Vivi Andasari http://people.bu.edu/andasari/courses/numericalpython/python.

html

Lecture 8: Euler’s Methods

Forward Euler

The formula for the forward Euler method is given by equation (8.2) in the lecture note
for week 8, as
yi+1 = yi + f (xi , yi )h . (1)
where f (xi , yi ) is the differential equation evaluated at xi and yi .

The implementation of this equation in Matlab or Python is quite straightforward, because

the calculation of yi+1 requires known quantities that have been previously calculated,
which is yi . Because of the straightforward calculation, the forward Euler method is also
referred to as the explicit method.

The first step is to calculate dy/dx = f (xi , yi ) at xi and yi . Evaluating differential equations
dy/dx is done by using feval function. In Matlab, the feval function is already provided
as a built-in function. However, in Python we have to manually define or create this
function. For example, we name the function “feval” and return the evaluated function
using the Python’s built-in command eval, as follows:
def feval ( funcName , * args ) :
return eval ( funcName ) ( * args )

We have briefly learned in lecture 7 on how to pass arguments using the asterisk symbol.
See file “withArgs firstOrderMethods.py”.

To test this function/method, create a .py file (name the file, for example: “example feval.py”).
See an example of this file in lecture 7.
import numpy as np
import math

def feval ( funcName , * args ) :

return eval ( funcName ) ( * args )

a = feval ( ’ np . cos ’ , np . pi )

1
Dr Vivi Andasari http://people.bu.edu/andasari/courses/numericalpython/python.html

print ( " Cos ( pi ) = " , a )

b = feval ( ’ math . log ’ , 1000 , 10 )

print ( " Log10 ( 1000 ) = " , round ( b ) )

def myFunc ( x ) :
return x ** 2

cin = 9 . 0
c = feval ( ’ myFunc ’ , cin )
print ( " If x = " , str ( cin ) , " then x ** 2 = " , c )

Running the file will result in the following outputs:

Cos ( pi ) = - 1 . 0
Log10 ( 1000 ) = 3
If x = 9 . 0 then x ** 2 = 81 . 0

The next step is to solve for yi+1 by directly implementing the above equation (1). We
create a function to implement the forward method. The function takes in several input
arguments to be used in calculating yi+1 . The input arguments consist of the ODE(s) to
be solved (we refer to it as func), initial condition(s) of the ODE(s) (yinit), the range of
interval where the ODE(s) to be evaluated (x range), and the step size (h). The function
definition with the input arguments mentioned above can be written as follows,
def forwardEuler ( func , yinit , x_range , h ) :

We have to make the function as general as possible, and that includes the ability to solve
a system of ODEs. The function needs to know how many ODEs to be solved and we
do this by calculating the length of yinit using the command len. The interval x range
also needs to be divided into smaller subintervals whose length is fixed and defined by the
step size h. The number of subintervals, along with the number of ODEs to be solved, all
will be used in identifying the number of iterations, hence, they must be of integer values.
After the function definition, we add the following lines:
numOfODEs = len ( yinit )
sub_intervals = int (( x_range [ - 1 ] - x_range [ 0 ] ) / h )

Next, we initialize variables x and y and fill them with the values from inputs passed in
the function definition,
x = x_range [ 0 ]
y = yinit

We also need to initialize other variables, xsol and ysol, which will be used as containers
to save the results of calculation of x and y after each iteration. The first element of each
of these containers is filled with the initial values of x and y, respectively,

2
Dr Vivi Andasari http://people.bu.edu/andasari/courses/numericalpython/python.html

xsol = [ x ]
ysol = [ y [ 0 ] ]

After initializations of all variables required in the calculation are done, the next thing do
is to calculate the derivative f (xi , yi ). As has been explained above, the derivative can be
evaluated using the feval function as follows,
yprime = feval ( func , x , y )

and apply the derivative in to the equation (1),

for j in range ( numOfODEs ) :
y [ j ] = y [ j ] + h * yprime [ j ]

Then we increase the position of x and store the new value in xsol,
x += h
xsol . append ( x )

The new value of y (which is the yi+1 ) is also stored in ysol. If there are multiple ODEs,
m
the value of each yi+1 (where m is the number of ODEs) that has been solved is stored one
after the other, such as [y11 , y12 , y13 , · · · , y21 , y22 , y23 , · · · ]. Likewise is the arrangement in ysol,
hence our algorithm needs to access each value of y and uses the following commands,
for r in range ( len ( y ) ) :
ysol . append ( y [ r ] )

and this completes one cycle of iteration. The calculation starts over again from the initial
loop, and so on. Our function returns the values of calculated x and y stored in xsol and
ysol.

The full implementation of the forward method can be done as follows (see “ex1 forwardEuler.py”),
def forwardEuler ( func , yinit , x_range , h ) :
’’’
This method performs the forward Euler method steps .
’’’
numOfODEs = len ( yinit ) # Number of ODEs
sub_intervals = int (( x_range [ - 1 ] - x_range [ 0 ] ) / h ) # Number of sub -
intervals

x = x_range [ 0 ] # Initialize variables x

y = yinit # Initialize variables y

# Containers for list of solutions

xsol = [ x ]
ysol = [ y [ 0 ] ]

3
Dr Vivi Andasari http://people.bu.edu/andasari/courses/numericalpython/python.html

for i in range ( sub_intervals ) :

yprime = feval ( func , x , y ) # Evaluate dy / dx

for j in range ( numOfODEs ) :

y [ j ] = y [ j ] + h * yprime [ j ] # Eq . ( 8 . 2 )

x + = h # Increase the x - step

xsol . append ( x ) # Save it in the container

for r in range ( len ( y ) ) :

ysol . append ( y [ r ] ) # Save all new y ’s

return [ xsol , ysol ]

When using Numpy, there are little changes such as the following (see “ex1 forwardEuler Numpy.py”):

• To create an empty list/array for container to store solutions,

xsol = np . empty ( 0 )
ysol = np . empty ( 0 )

• To fill the first element of the container,

xsol = np . append ( xsol , x )
ysol = np . append ( ysol , y )

• To save the solution x in xsol,

xsol = np . append ( xsol , x )

and ysol,
ysol = np . append ( ysol , y [ r ] )

Backward Euler

The formula for the backward Euler method is given by equation (8.19) in the lecture note,

yi+1 = yi + f (xi+1 , yi+1 )h . (2)

In this formulation, we first need to calculate the derivative at i + 1, i.e.,

0
yi+1 = f (xi+1 , yi+1 ) , (3)

4
Dr Vivi Andasari http://people.bu.edu/andasari/courses/numericalpython/python.html

which is a bit tricky because at present (at i) we obviously don’t have the future value (at
i + 1) of yi+1 . Hence, the backward method is referred to as being implicit. To solve this
problem, we do calculation by hand using an example equation,

dy
= 3(1 + x) − y . (4)
dx

Let’s take i = 0, and the equation (2) becomes

y1 = y0 + f (x1 , y1 )h = y0 + y10 h , (5)

where, using the example equation (4), the derivative of y1 is

y10 = f (x1 , y1 )
= 3(1 + x1 ) − y1 (6)

Substituting y1 of the equation (5) into equation (6) gives us

y10 = 3(1 + x1 ) − (y0 + y10 h) .

Equating all terms of y10 on the left hand side, we get

y10 + hy10 = 3(1 + x1 ) − y0
y10 (1 + h) = 3(1 + x1 ) − y0
 
0 1
y1 = 3(1 + x1 ) − y0 , (7)
1+h
0
from which, we can formulate yi+1 as
0 1
yi+1 = f (xi+1 , yi ) . (8)
1+h

0
In Python, the implementation of yi+1 can be represented by
yprime = feval ( func , x +h , y )
yp = mult ( yprime , 1 / ( 1 + h ) )

0
where, the formulation of yp is equal to yi+1 . Here we use a function called mult to perform
element-wise multiplication operation of yprime (which is a vector or Python’s list of single
value) and (1/(1+h)) (which is a scalar). The function mult is defined as follows,
def mult ( vector , scalar ) :
newvector = [ 0 ] * len ( vector )
for i in range ( len ( vector ) ) :
newvector [ i ] = vector [ i ] * scalar
return newvector

5
Dr Vivi Andasari http://people.bu.edu/andasari/courses/numericalpython/python.html

Without defining and using the mult function, direct multiplication (using the asterisk) of
list and scalar will result in repetition of the whole list, such as
>>> a = [ 1 . 0 ]
>>> b = 3 * a
>>> b
[ 1 .0 , 1 .0 , 1 . 0 ]

However, when using Numpy, the product of vector/array/list and scalar is straight for-
ward, hence there is no need to involve the mult function,
yprime = feval ( func , x +h , y )
yp = yprime * ( 1 / ( 1 + h ) )

or, even directly with one line of assignment:

yprime = feval ( func , x +h , y ) / ( 1 + h )

The full backward method can be implemented in a function as follows,

def backwardEuler ( func , yinit , x_range , h ) :
numOfODEs = len ( yinit )
sub_intervals = int (( x_range [ - 1 ] - x_range [ 0 ] ) / h )

x = x_range [ 0 ]
y = yinit

xsol = [ x ]
ysol = [ y [ 0 ] ]

for i in range ( sub_intervals ) :

yprime = feval ( func , x +h , y )

yp = mult ( yprime , ( 1 / ( 1 + h ) ) )

for j in range ( numOfODEs ) :

y [ j ] = y [ j ] + h * yp [ j ]

x += h
xsol . append ( x )

for r in range ( len ( y ) ) :

ysol . append ( y [ r ] )

return [ xsol , ysol ]

6
Dr Vivi Andasari http://people.bu.edu/andasari/courses/numericalpython/python.html

Heun’s Method

This method is formulated as in equation (8.15) in the lecture note,

h h 0
yi+1 = yi + f (xi , yi ) + f (xi+1 , yi+1 ), (9)
2 2
where
0
yi+1 = yi + f (xi , yi )h . (10)

To break it down, in our Python code the derivative f (xi , yi ) is represented by,
y0prime = feval ( func , x , y )

while the second term on the right hand side of equation (10) or multiplication of f (xi , yi )
and h is represented by
k1 = mult ( y0prime , h )

the predictor equation (10) is coded as (element-wise addition of arrays)

ypredictor = [ u + v for u , v in zip (y , k1 ) ]

or if using Numpy, the element-wise addition operation can be done directly:

ypredictor = y + k1

0
the derivative with predictor equation f (xi+1 , yi+1 ) is represented by,
y1prime = feval ( func , x +h , ypredictor )

and the equation for Heun’s method in (9) is represented by

for j in range ( numOfODEs ) :
y [ j ] = y [ j ] + ( h / 2 ) * y0prime [ j ] + ( h / 2 ) * y1prime [ j ]

The full function for implementing Heun’s method is given by

def HeunsMethod ( func , yinit , x_range , h ) :
numOfODEs = len ( yinit )
sub_intervals = int (( x_range [ - 1 ] - x_range [ 0 ] ) / h )

x = x_range [ 0 ]
y = yinit

xsol = [ x ]
ysol = [ y [ 0 ] ]

7
Dr Vivi Andasari http://people.bu.edu/andasari/courses/numericalpython/python.html

for i in range ( sub_intervals ) :

y0prime = feval ( func , x , y )

k1 = mult ( y0prime , h )

ypredictor = [ u + v for u , v in zip (y , k1 ) ]

y1prime = feval ( func , x +h , ypredictor )

for j in range ( numOfODEs ) :

y [ j ] = y [ j ] + ( h / 2 ) * y0prime [ j ] + ( h / 2 ) * y1prime [ j ]

x += h
xsol . append ( x )

for r in range ( len ( y ) ) :

ysol . append ( y [ r ] )

return [ xsol , ysol ]

Midpoint Method

The midpoint method is formulated as in equation (8.16) of the lecture note,

yi+1 = yi + hf (xi+ 1 , yi+ 1 ) , (11)

2 2

where
h
yi+ 1 = yi + f (xi , yi ) . (12)
2 2

Python implementation of the midpoint method is almost similar to the implementation

of Heun’s method, except for the second term on the right hand side of equation (12) or
h
2
f (xi , yi ) is coded as
k1 = mult ( y0prime , h / 2 )

and the derivative in equation (11) or f (xi+ 1 , yi+ 1 ) is represented by

2 2

y1prime = feval ( func , x + h /2 , ypredictor )

and the midpoint equation (11),

for j in range ( numOfODEs ) :
y [ j ] = y [ j ] + h * y1prime [ j ]

8
Dr Vivi Andasari http://people.bu.edu/andasari/courses/numericalpython/python.html

The full function for midpoint method is given by

def midpoint ( func , yinit , x_range , h ) :
numOfODEs = len ( yinit )
sub_intervals = int (( x_range [ - 1 ] - x_range [ 0 ] ) / h )

x = x_range [ 0 ]
y = yinit

xsol = [ x ]
ysol = [ y [ 0 ] ]

for i in range ( sub_intervals ) :

y0prime = feval ( func , x , y )

k1 = mult ( y0prime , h / 2 )

ypredictor = [ u + v for u , v in zip (y , k1 ) ]

y1prime = feval ( func , x + h /2 , ypredictor )

for j in range ( numOfODEs ) :

y [ j ] = y [ j ] + h * y1prime [ j ]

x += h
xsol . append ( x )

for r in range ( len ( y ) ) :

ysol . append ( y [ r ] )

return [ xsol , ysol ]