EMAT10006: Further Computer Programming resit assignment
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
EMAT10006: Further Computer Programming resit assignment
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
dx
dt
For example the ODE might be
dt = x, (2)
in which case the function f is defined by f (x, t) = x. The Euler method is described in the Eng Maths 1 notes (and in many other places online) and is defined by the rule
xn+1 = xn + f (xn , 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 (xn , tn )
k2 = f ╱xn + k1 , tn + 、
k2 = f ╱xn + k2 , tn + 、
k3 = f (xn + k3 h, tn + h)
h
xn+1 = xn + (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
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 x = 1 when t = 0 and we want to find the value of x 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 x when t = 1 (the values of x 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
dx
dt
dy
dt
with initial conditions x(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
0 |
2 |
4 |
6 |
8 |
10 |
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
dx
dt
with initial condition x(0) = 1. The true solution is x = et and its value at t = 1 is x(1) = e1 = e = 2.71 . . . . Use your solveto function with both the Euler and RK4 methods to estimate x(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
Include both the code for the plot and the plot itself in your Jupyter notebook.
10 Markscheme
In this part of the assignment the markscheme is really just looking at the code that you write:
● Does the code work?
● Does it pass the tests?
● Is it written reasonably with good use of loops or functions or anything else?
● Does the code have reasonable comments and/or docstrings to explain what it does?
● Is the code nicely formatted with reasonable variable names?
The Jupyter notebook can be used to show any plots you make. You do not need to write a long report just enough to show a few plots of your code working. Feel free to make extra plots showing the solutions of other differential equations. Primarily the markscheme here is concerned with the code though so you do not need to spend a long time on writing.
2022-08-08