1    Outline

There are two parts to this assignment:

1. Write some code to implement numerical methods to solve ordinary dif- ferential equations (ODEs) using the methods covered in Eng Maths 1.

2. Write a program that can solve ODEs related to a class of electrical cir- cuits.

This document describes only the first part and a separate document describes the second part. To pass this assessment it is enough to complete the first part but to get a higher mark you should attempt the second part which is more open ended.

2    Submission

You should use git and github to keep your work in a git repository (repo) as described in the course notes. Your submission should be a zip file containing the git repo.  Call your repository odesolve and then the submission should be a file odesolve.zip.  Inside the zip file will be the folder odesolve and inside  there will be all of your code files and also the  .git folder.   The zip file MUST be a zip file and NOT some other format (i.e.  NOT 7z or RAR or anything like that). The zipfile should be submitted to Blackboard.

We need to confirm your usage of git so it is essential that the .git directory is in the odesolve folder. It should be possible to download the zip file, extract it and cd into the odesolve directory and then run git commands like git status or git log to see the commits that you have made with git.  It should look like this:

#  Download  the   z i p   f i l e

#  E x t r a c t   the   z i p   f i l e

$  cd   o d e s o l v e

$   g i t   s t a t u s

On  branch  main

n o t h i n g   to  commit ,   working   d i r e c t o r y   c l e a n

IT  IS  ESSENTIAL  THAT  WE  CAN  USE  GIT  TO  SEE  THE

COMMITS THAT YOU HAVE MADE.

Inside the repo there should be the following:

·  A README.md file that explains what is in the repo.

·  A file odesolve.py that contains the main code that you will write for this

part of the assignment.

·  The file test odesolve.py the is provided with this specification on Black-

board. assignment.

·  A Jupyter notebook that shows how to import and use the functions that

are defined in odesolve.py.  This should be used to show any plots of e.g.

solutions to ODEs.

·  Any files that are part of the circuits assignment.

3    Problem summary

The goal here is to make some functions for solving  ODEs using the Euler

method and the 4th order Runge-Kutta method (RK4).  Test code is provided

in the file test example.py that you should include in your repository.  A file

odesolve.py is provided that contains the functions that you should write but

the functions are missing the function body like this:

#  odesolve . py

#

#  Author :  <i n s e r t   name>

#  Date :       <i n s e r t   date>

#  Description :  <i n s e r t   d e s c r i p t i o n >

#

# You  should   f i l l   out   the   code   f o r   the   f u n c t i o n s   below   so   that   they   pass   the #  t e s t s   in   t e s t  o d e s o l v e . py

def   e u l e r ( f ,   x ,   t ,   h ) :

””” Perform  one   step   of   the   Euler  method”””

pass

def   rk4 ( f ,   x ,   t ,   h ) :

””” Perform  one   step   of   the  RK4  method”””

pass

def   s o l v e t o ( f ,   x1 ,   t1 ,   t2 ,  hmax ,   method=e u l e r ) :

”””Use  many  s t e p s   of  method  to   get   from  x1 , t1   to  x2 , t2 ”””

p a s s

d e f   o d e s o l v e ( f ,  X0 ,   t ,  hmax ,   method=e u l e r ) :

”””Compute  the   s o l u t i o n   a t   d i f f e r e n t   v a l u e s   o f   t ””” p a s s

The test script is provided as test odesolve.py.  To run the tests just run that file. When you first run the tests you should see that they fail:

/ $  python   t e s t  o d e s o l v e . py

Checking   t e s t  e u l e r    . . .

FAILED

0   t e s t s   p a s s e d   out   o f  6   t e s t s

t e s t  e u l e r   f a i l e d   with   e x c e p t i o n :

Traceback   ( most   r e c e n t   c a l l   l a s t ) :

F i l e   ” t e s t  o d e s o l v e . py ” ,   l i n e   1 3 5 ,   i n  <module>

r a i s e   exc   from  None

F i l e   ” t e s t  o d e s o l v e . py ” ,   l i n e   1 1 8 ,   i n  <module>

fu n c ( )

F i l e   ” t e s t  o d e s o l v e . py ” ,   l i n e   1 7 ,   i n   t e s t  e u l e r

a s s e r t   e u l e r ( f ,   0 ,   0 ,   1 )  ==  0

A s s e r t i o n E r r o r

Your goal is to make those tests pass one by one.  We can fix the first test by adding the code for the Euler method:

d e f   e u l e r ( f ,   x ,   t ,   h ) :

r e t u r n  x  +  f ( x ,   t ) * h

Now if you rerun the tests you should see:

/ $  python   t e s t  o d e s o l v e . py

Checking   t e s t  e u l e r    . . .

PASSED

Checking   t e s t  s o l v e t o    . . .

FAILED

1   t e s t s   p a s s e d   out   o f  6   t e s t s

t e s t  s o l v e t o   f a i l e d   with   e x c e p t i o n :

Traceback   ( most   r e c e n t   c a l l   l a s t ) :

F i l e   ” t e s t  o d e s o l v e . py ” ,   l i n e   1 3 5 ,   i n  <module>

r a i s e   exc   from  None

F i l e   ” t e s t  o d e s o l v e . py ” ,   l i n e   1 1 8 ,   i n  <module>

fu n c ( )

F i l e   ” t e s t  o d e s o l v e . py ” ,   l i n e   3 1 ,   i n   t e s t  s o l v e t o

a s s e r t   s o l v e t o ( f ,   0 ,   0 ,   1 ,   1 )  ==  0

A s s e r t i o n E r r o r

Here the important part is 1  tests  passed out of 6  tests which shows that the first test has passed but the other 5 have still failed. The steps below describe

how to add the rest of the code that you should write yourself in order to make all of the tests pass.

4    Euler method

We are concerned here with solving ordinary differential equations (ODEs).  A first order ODE in standard form looks like

= f (t. t)． dt

For example the ODE might be

dt

in which case the function f is defined by f (t. t) = t.

(1)

(2)

The Euler method is

described in the Eng Maths 1 notes (and in many other places online) and is defined by the rule

tn\1  = tn + f (tn . tn ) h.                                      (3)

where h is the stepsize. In order to make a numerical solver we want to have a function called euler that can be used like this:

#  e u l e r c h e c k . py

from   o d e s o l v e   import   e u l e r

#  D e f i n e s   the  RHS  o f   our  ODE

d e f   f ( x ,   t ) :

r e t u r n  x

#  I n i t i a l   c o n d i t i o n s

t 0  =  1

x0  =  1

h  =  0 . 5

t 1  =  t 0  +  h

x1  =  e u l e r ( f ,   x0 ,   t0 ,   h )

p r i n t ( ’ S o l v i n g   the  ODE  dx/ dt  =  x ’ )

p r i n t ( ’ I n i t i a l   c o n d i t i o n  x  = ’ ,  x0 ,    ’ when  t   = ’ ,   t 0 )

p r i n t ( ’ One  s t e p   o f   the   E u l e r  method  with  a   s t e p s i z e   o f  h  = ’ ,  h ) p r i n t ( ’ Gives  an   e s t i m a t e   o f  x  = ’ ,  x1 ,    ’ when  t   = ’ ,   t 1 )

If the  euler  function is working correctly then the output from running the above code should look like:

$  python   e u l e r c h e c k . py

S o l v i n g   the  ODE  dx/ dt  =  x

I n i t i a l   c o n d i t i o n  x  =  1  when  t  =  1

One  s t e p   o f   the   Euler  method  with  a   s t e p s i z e   o f  h  =  0 . 5 Gives  an   e s t i m a t e   o f  x  =  1 . 5  when  t  =  1 . 5

The code needed to make the euler function work correctly is shown above but to make the other tests pass you will need to write the rest of the code yourself.

5    RK4

The rk4 function should implement the 4th order Runge Kutta method which is more accurate than the Euler method. One step of RK4 looks like:

k1  = f (tn . tn )

k2  = f 4tn + k1  . tn + 、

k2  = f 4tn + k2  . tn + 、

k3  = f (tn + k3 h. tn + h)

h

tn\1  = tn + (k1 + 2k2 + 2k3 + k4 )

Adjusting the code above to use the rk4 function should output:

/ $  python   t . py

S o l v i n g   the  ODE  dx/ dt  =  x

I n i t i a l   c o n d i t i o n  x  =  1  when  t  =  1

One  s t e p   o f   the  RK4  method  with  a   s t e p s i z e   o f  h  =  0 . 5 Gives  an   e s t i m a t e   o f  x  =  1 . 6484375  when  t  =  1 . 5

(4)

(5)

(6)

(7)

(8)

6    solveto

Once you have the RK4 method working you should make the function solveto . This function should be used to make a number of steps with either the Euler or RK4 method. For example, if we have an ODE with initial condition t = 1 when t = 0 and we want to find the value of t when t = 1 using a maximum

step size of 0．1 then we could use

d e f   f ( x ,   t ) :

r e t u r n  x

x0  =  1

t0  =  0

hmax  =  0 . 1

t1  =  1

x1  =  s o l v e t o ( f ,   x0 ,   t0 ,   t1 ,  hmax)

The solveto function will perform many Euler steps to go from t0 to t1 in steps of size 0.1. There will be an Euler step from t = 0 to t = 0．1, then from t = 0．1 to t = 0．2 and so on until t =  1.   The solveto function will then return an estimate for the value of t when t = 1 (the values of t at the intermediate steps are not returned).  It should be possible to use solveto with either the Euler method or the RK4 method like this:

from   o d e s o l v e   import   e u l e r ,   rk4 ,   s o l v e t o

x1  =  s o l v e t o ( f ,   x0 ,   t0 ,   t1 ,  hmax ,   e u l e r )  #  Uses   E u l e r  method x1  =  s o l v e t o ( f ,   x0 ,   t0 ,   t1 ,  hmax ,   rk4 )      #  Uses  RK4  method

x1  =  s o l v e t o ( f ,   x0 ,   t0 ,   t1 ,  hmax)                  #  Uses   Eu l e r  method  by   d e f a u l t

If the stepsize does not evenly divide the interval from t0 to t1 then solveto must make a final smaller step. For example to from t = 0 to t = 1 in steps of size t = 0．3 then there would be 3 steps of size 0．3 and a final step of size 0．1. It is important to get this right for the accuracy of the method.

7    Arrays

Once you have solveto working you should check that your functions can all work with numpy arrays. This is needed to be able to handle systems of ODEs. For example, suppose we want to solve the system

dt

dt

dy

dt

with initial conditions t(0) = 1 and y(0) = 0.  Then we should be able to do that with the following code:

import  numpy  as  np

from   o d e s o l v e   import   s o l v e t o

#  D e f i n e s   the  RHS  o f   our  ODE

d e f   f (X,   t ) :

x ,   y  = X

dxdt  =  y

dydt  = -x

r e t u r n  np . a r r a y ( [ dxdt ,   dydt ] )

#  I n i t i a l   c o n d i t i o n s

t0  =  0

x0  =  1

y0  =  0

X0  =  np . a r r a y ( [ x0 ,   y0 ] )

h  =  0 . 5

t 1  =  1

X1  =  s o l v e t o ( f ,  X0 ,   t0 ,   t1 ,   h )

p r i n t (X1)

This should output

[   0 . 7 5   - 1 .     ]

8    odesolve

Finally you should make a function odesolve. The purpose of this function is to give many different values of the solution for 北 suitable for a plot. For example:

import  numpy  as  np

import   m a t p l o t l i b . p y p l o t   as   p l t

from   o d e s o l v e   import   o d e s o l v e

#  D e f i n e s   the  RHS  o f   our  ODE

d e f   f (X,   t ) :

x ,   y  = X

dxdt  =  y

dydt  = -x

r e t u r n  np . a r r a y ( [ dxdt ,   dydt ] )

#  I n i t i a l   c o n d i t i o n s

x0  =  1

y0  =  0

X0  =  np . a r r a y ( [ x0 ,   y0 ] )

h  =  0 . 0 1                                             # max  s t e p   s i z e

t  =  np . l i n s p a c e ( 0 ,   1 0 ,   100)  #  t i m e s   t o   p l o t   the   s o l u t i o n Xt  =  o d e s o l v e ( f ,  X0 ,   t ,   h )

#  Use   . T  t o   t r a n s p o s e   the   a r r a y

p l t . p l o t ( t ,   Xt . T)

p l t . s a v e f i g ( ’ shm . pdf ’ )

p l t . show ( )

This should give a plot like the one shown below:

1.0

0.5

0.0

0.5

1.0

Include code and a plot in your notebook.  The odesolve function receives an array of times t and an initial condition X0 which specifies the value of the variables at time t [0] .  It should use the solveto function repeatedly in a loop to make an 2D array representing the solution.

9    Errors

The final step is to consider the errors of the numerical methods for different step sizes. Consider the following ODE

dt

dt

with initial condition t(0) = 1.  The true solution is t = et  and its value at t = 1 is t(1) = e1  = e = 2．71．．．. Use your solveto function with both the Euler and RK4 methods to estimate t(1).  Comparing with the true value e you can compute the errors of the two methods. Produce a plot like this one:

Comparison of errors for Euler and RK4

10  5                                 10  4                                 10  3                                 10  2                                 10  1

h