Computational Physics I Computational Astrophysics I PH32005 PH32010, Semester 2 201920
Module leaders:
Prof. Ulrich Zachariae Dr Scott Gregory
Fridays 25pm, Durajczik Lab, Harris Bldg.
Email: u.zachariaedundee.ac.uk sgregory001dundee.ac.uk
Email ONLY for general questions about lecture content etc, NOT for help with assignments!
Class Week 1, 24 Jan 2, 31 Jan
3, 7 Feb 4, 14 Feb 5, 21 Feb 6, 28 Feb 7,6Mar
8, 13 Mar
9, 20 Mar 10, 27 Mar 11, 3 Apr
Overview
Exercise Title
Solution of ODEs with predictor corrector methods
Solution of 1D ODEs with finite difference methods
Solution of 1D PDEs with an implicit FD method
Planetary Motion
No class
Magnetospheric models root finding and field line tracing
Solution of Laplaces equation with the GaussSeidel overrelaxation method
The MonteCarlo method and random walks Simulation of a 2D ideal gas
Molecular Dynamics Simulation
Class test
Module Principles
about 45 min 15 introductory lecture
4 class exercises
45 minutes to solve questions, then break 15 min
up to 1 hour to finish questions, revise comment code and save
programming platform to be used is MATLAB
Fri 25pm
Module principles Class Exercises
3 Questions to be solved during workshop hours, can also be completed afterwards
Workshop exercises to be submitted to Blackboard by Wednesday 6pm of the following week.
MATLAB code to be submitted as .m file or files it is necessary to comment code to explain what you did !
Graphs or other graphical output of the code and answers to specific questions to be submitted in a pdf as figures with captions, tables or text, respectively.
The more you tell me what you do in each step the better! Plain code without explanations and comments will lead to low grades even if it works.
Adapted from WIKI:
MATLAB
MATLAB matrix laboratory is a numerical computing environment and programming language. MATLAB allows matrix manipulations, plotting of functions and data, implementation of algorithms, creation of user interfaces, and interfacing with programs written in other languages, including C, C, Java, Fortran and Python.
MATLAB is intended primarily for numerical computing, and has around one million users across industry and academia. MATLAB users come from various backgrounds of engineering, science, and economics.
MATLAB is widely used in academic and research institutions as well as industrial enterprises.
MATLABR2008aforLinux User: Baosheng english Wikipedia
MATLAB
MATLAB
MATLAB code: Simple, many builtin functions available esp. for science engineering Matrix manipulation, data processing graphing particularly easy
Reading List:
Attendance Sheets
M M Woolfson, G J Pert: An introduction to computer simulation Oxford University Press
P Harrison: Computational methods in Physics, Chemistry and Biology, Wiley
Klein and Gudunov: Introductory computational physics Cambridge University Press
M M Woolfson and M S Woolfson: Mathematics for physics Oxford University Press
Class 1 Solution of ODEs with predictor corrector methods
In many situations in computational physics we are faced with finding a solution to a generalised firstorder Ordinary Differential Equation ODE of the form:
Boundary condition:
Recap ODE: 1 independent variable; First order: first derivative of any variable
Unless the function fx,y ft,y is relatively simple, it may not be possible to integrate Eqn. 1 directly analytical solution, and we need to seek numerical solutions.
In this laboratory we will consider several numerical approaches and compare their computational complexity and underlying accuracy.
Eqn. 1
Solution 1 The Euler method
The simplest method of solving a first order ODE is the Euler method. The principle of
this approach is illustrated in the Figure. For the Euler method we consider a short section of the solution curve yx near the point x0,y0 here: yt near t0,y0.
We find the next point on the curve
yt by moving a small distance h along The taxis and finding the new y value by approximating the curve as a straight line with a gradient given by the original differential equation at the starting
point.
San Diego State University
Solution 1 The Euler method
This may be generalised for the mth step taken along the t or x axis to:
Eqn. 2a
San Diego State University
Eqn. 2b
Solution 1 The Euler method
It is intuitively clear
that unless h is very small, the solution obtained
by the Euler method will very quickly deviate from the true solution curve yt.
The error in approximating over a fixed interval h is Oh.
San Diego State University
Solution 2 The Euler predictor
corrector method Heuns method
A significant improvement in the accuracy of the Euler method is provided by the EPC methodology, whereby one step in the integration now uses an estimate of the slope of the chord connecting the interval endpoints.
Wikipedia
The slope of this chord is estimated as the average of the slopes of the curves at the two ends, so that this approach requires two function evaluations per integration step compared to the single evaluation of the Euler method.
The accuracy of the method is improved and is now Oh2.
Solution 2 The Euler predictor
corrector method Heuns method
To find the slope of the curve at the two ends, we use an Euler predictor step:
with
… and then a corrector
step based on the
trapezoidal rule: Wikipedia
… Oh2 so by halving the interval h the error is reduced by a factor of 4 compared to the Euler method.
Eqn. 3a
Eqn. 3b
Solution 3 The RungeKutta method
The RungeKutta method is a very popular method for a simple yet accurate numerical solution to ODEs. It requires 4 function evaluations per step. These involve obtaining estimates of the increment in y at the initial point, midpoint, and endpoint locations of the interval. For the point xn, yn, the four estimates are calculated as:
1 and 4 are essentially Euler estimates of the y increment based upon the slope at the initial and end points respectively. 2 and 3 are both estimates of the y increment based on the slope at x xn h2 .
Eqns. 4ad
Solution 3 The RungeKutta method
Then, a weighted average is taken in the ratio 1:2:2:1 cf. Simpsons rule in order to estimate ym at the next x location:
Eqn. 5
The error of the RK method is Oh4.
Solution 4 The AdamsMoulton method
The previous methods 13 are linear singlestep methods, i.e. information from previous approximation steps is discarded. Multistep methods try to improve efficiency by using the information gained in previous steps. AM is one of a family of more accurate and sophisticated approximation techniques. It requires values from 4 previous approximation steps.
The error is similar to RungeKutta, but AM only requires two evaluations per step. AM is a linear multistep method, i.e. a linear combination of previous points and derivatives is used.
Solution 4 The AdamsMoulton method
Let us assume that the required 4 starting values are obtained from RK: x0,y0…x3,y3
A first predictor stage is then implemented by fitting a cubic Lagrange polynomial to these 4 points, which is then used to obtain an estimate YAM for y4 at x x4.
With the predicted value YAM for y4, a good estimate may subsequently be made for fx4,YAM. This polynomial may then finally be used to obtain a better estimate of y4.
Eqn. 6
Module principles Class Exercises
MATLAB code to be submitted as .m file or files extensively commented code required to explain what you did !
Graphs or other graphical output of the code and answers to specific questions to be submitted in a pdf as figures with captions, tables or text, respectively.
The more you tell me what you do in each step the better! Plain code without explanations and comments will lead to low grades even if it works.
Class test on 03 Apr
Class Exercises Overview
EXAMPLE:
Implementation of Euler ODE initial value problem by a MATLAB function taking the variables f function, xinit starting value for x, yinit starting value for y and n number of increments in evaluation interval, and returning vectors of x and y:
function x,yeulerbackwardf,xinit,yinit,xfinal,n
end
Euler approximation for ODE initial value problem Euler backward method
File prepared by DJT 21012017
Calculation of h from xinit, xfinal, and n hxfinalxinitn;
Initialization of x and y as column vectors
xxinit zeros1,n;
yyinit zeros1,n;
Calculation of x and y by Euler method, Eqn X in lecture 1 for i1:n
xi1xih; ynewyihfxi,yi; yi1yihfxi1,ynew;
end
plot x,y,dr;
The function plots an x,y graph of the independent and dependent variables x, y with red diamond symbols, no line between function points plot x,y,dr.
Class Exercises
The aim of this class is to solve a differential equation that has the format of Eqn 1. using the Euler, predictorcorrector and RK methods. The equation to be solved has fx,y y
so that:
The analytical solution is accordingly:
Eqn. 7
The numerical output may thus be compared directly with the known exact values. The equation is to be solved for the interval 0 x 1 . The evaluation step size h is to be set to h 0.1 . The point x0, y0 is to be taken as 0, 1 .
1
where the x and y values are stored in appropriate vectors and can be referenced by an appropriate running variable. Plot the numerical solution as individual symbols on an xygraph and compare this with the analytical solution on the same graph dashed curve. Also print the numerical and analytical solutions side by side with x in your final script e.g. in a table. Save your code as an individual .m file.
2 Solve Eqn 7 in MATLAB using the predictorcorrector method. Write code to implement Eqns 3a and 3b, and put x and y values into vectors. Plot and tabulate x and y as in 1.
3 Solve Eqn 7 in MATLAB using the RungeKutta method. Write code to implement Eqns 4ad and 5, and put x and y values into vectors. Plot and tabulate x and y as in 1.
Solve Eqn 7 in MATLAB using the Euler method. Write code to implement Eqns 2a and 2b
Class Exercises
4. Investigate the relative accuracy of the Euler, predictorcorrector and RK methods as the evaluation step is both decreased and increased by a factor of 2. Include a MATLAB plot of the relative inaccuracieserror magnitudes for each method and each step size.
Workshop exercises and homework to be submitted to Blackboard by Wednesday 6pm of the following week.
Submit code as .m files other output graphs, tables, text to be submitted separately as .pdf file!
Include your name in the filename please.