MATH 3001 Computational applied maths 2019/20
MATH3001: Computational Applied Maths
Contents
1 Introduction
1
Supervisors:
Adrian Barker (A.J.Barker@leeds.ac.uk) Stephen Griffiths (S.D.Griffiths@leeds.ac.uk) Daniel Read (D.J.Read@leeds.ac.uk) Rob Sturman (R.J.Sturman@leeds.ac.uk)
October 22, 2019
1.1 Writingyourreport ………………………………. 2 1.2 TypesettingusingLATEX…………………………….. 2 1.3 Programming …………………………………. 3
2 The Rayleigh–Ritz method 4
3 Fast Fourier transforms 5
4 How can we accurately evolve planetary orbits using geometric integrators? 8
5 Solving the heat equation 10
6 Using random numbers to integrate 13
7 Determining the period of a pendulum 16
1 Introduction
Most mathematical problems that arise in academic research or industry are not solvable without the use of computers. This project will explore the application of computational methods to solve various prob- lems in applied mathematics, including solving partial differential equations and evaluating integrals by various means.
This project will consist of 3 independent investigations, chosen from those described in sections 2–7. Each investigation has well-defined tasks and the opportunity for independent exploration (i.e., for you to set your own tasks). These tasks fall into three main categories:
1. computational – requiring computer programs to be written, or used;
2. analytical – requiring mathematical calculations, on paper;
3. research – requiring you to read parts of books and research papers, and summarise the results.
To complete all the tasks in detail in all of your 3 chosen investigations would probably constitute an M.Sc. dissertation worth 60 credits (or more). So, you need not complete all the tasks, and those that you do complete could be investigated to different degrees. All that is required is that your final report is of the right length, with each of your 3 investigations presented as a coherent set of results.
1
MATH 3001 Computational applied maths 2019/20
1.1 Writing your report
To quote from the module handbook, ‘the length of the project report should be between 20–25 pages (including the contents page and bibliography but excluding the title page and appendices)’. It will be natural to structure the report into three main sections (i.e., one for each investigation). Each of these sections should contain a short introduction (describing the background to the investigation and giving an overview of what is included), individual subsections for each part of your results, and a conclusion (summarising what has been achieved and what more could be done). You will probably also want to include a table of contents and a list of references to cover the entire report. Note that giving and appropriately citing additional references is an important part of the projects: a good report would have more than 10 cited works (books, research paper, websites). Finally, all python codes that have been used to generate results given in your report should be included as appendices, so that the assessors can re-run the codes and isolate errors, if necessary.
It would be natural to aim for each investigation to be written up over 7–9 pages. However, it is fine to go into more depth in some investigations than others. For example, your three investigations might appear in the report as 9+9+6 pages. We would not be happy with a split of 11+11+2 pages!
It is hard to write a good scientific report. All assertions and results need to be supported by either an explicit proof, tables of data, graphs, or a citation. Note that Figures (and Tables) must be referred to in the text by their number – so one would write ‘Figure 2’, rather than ‘the Figure at the top of the next page’, for example. Each Figure (and Table) must also have a caption, which typically starts with a single sentence describing the contents, followed by more detailed information (giving parameters, or explaining line colours, etc.). Even when all the components are in place, you will need to read and rewrite your report several times to obtain a polished product. If this is not done, it is almost inevitable that you will not score well for your report, however interesting is the mathematical content.
1.2 Typesetting using LATEX
As you may know, LATEX is a typesetting package designed for writing mathematical reports, and it is recommended that you use it for this project. You can access the software in different ways:
1. use LATEX on the University machines – the best version is called TeXworks.
2. Download LATEX freely from the web. Nice versions are TeXworks, or TeXshop (for Mac).
3. Use a web interface such as www.overleaf.com. This avoids the hassle of installing software, and you can work on your files from anywhere (including from tablets and smartphones).
There are many LATEX tutorials along with other documentation on the web. For example, you may find it helpful to consult
http://www.maths.leeds.ac.uk/latex/
http://www.maths.leeds.ac.uk/latex/texworks23.pdf
When including your python codes in your report appendix, it will be best to use the LATEX listings package. This can directly import code from your python files using commands of the form
\lstinputlisting{source_filename.py}
i.e., there is no need to copy and paste to your LATEX file.
2
MATH 3001 Computational applied maths 2019/20
1.3 Programming
You should write your programs in the programming language python, which is used in other modules in years 1 and 2. It is also a versatile language that is increasingly used in a variety of workplaces, so it is a good choice for employability. There is a good online course at
https://www.codecademy.com/learn/learn-python-3
There are various ways to run python on different machines. Sometimes it can be invoked by typing python at a command-line shell; for example, on a Mac, you can start a Terminal window (in Applica- tions), and then just type python. There is also a common interactive environment called IDLE, which is available for most platforms. However, it is easiest to interact with python via one of many nice graph- ical interfaces, such as spyder, which is part of the anaconda package that you can download for free (for Windows, Mac, or UNIX) from https://www.anaconda.com/distribution/ (or googleanaconda python).Itisbesttodownloadpython3ratherthanpython2.
Since anaconda is used for teaching python at Leeds, it is also installed on University machines. You can find it by searching from the Start Menu, or by looking in All Programs / Departmental Software, although it is often tucked away in a further subfolder (e.g., All Programs / Departmental Software / Engineering / Computing, or All Programs / Departmental Software / MAPS / Maths). Again, ensure that you choose python3 rather than python2.
Once you have figured out how to run a python code, save this simple code as test.py:
Run the code; the output should be the number 3, printed somewhere on screen. There is no point going further until you can get this simple code to work!
a=1
b=2 c=a+b print ( c )
3
MATH 3001 Computational applied maths 2019/20
2 The Rayleigh–Ritz method
A circular drum supports oscillating solutions, obeying Bessel’s equation of order zero:
d2y + 1 dy + λy = 0, dx2 xdx
for 0 ≤ x ≤ 1, subject to the boundary conditions y(1) = 0, and y non-singular at x = 0. As discussed in MATH2375 and MATH2340, λ = ω2/c2, where ω is the frequency of oscillation of the drum, and which can take many different values corresponding to the natural harmonics of the instrument, and c is the wave speed. However, the higher harmonics typically decay quickly, so the lower frequencies are of greater interest.
The exact solution to Bessel’s equation is given by
y(x)=J λ1x, 02
1
where λ2 is a zero of the Bessel function J0.
Find the first four zeros of J0(x), and plot the four corresponding solutions y(x) to Bessel’s equation. Note that the library scipy.special contains the functions jn zeros and j0, which should make this straightforward.
To use the Rayleigh-Ritz method to estimate solutions of Bessel’s equation, we express the equa- tion in Sturm-Liouville self-adjoint form, and estimate the eigenvalues. The equivalent Sturm-Liouville problem is
d dy
−dx xdx =λxy,
again for 0 ≤ x ≤ 1, with the same boundary conditions as above. Recall from MATH2650 (or consult a suitable textbook, such as Mathematical Methods for Physicists, by G. Arfken) that we set
1 1
I = xy′2dx, J = xy2dx
00
The extrema of I/J are exactly the eigenvalues λ of the Sturm-Liouville, and since in Bessel’s equation we have a finite lowest eigenvalue λ0, any evaluation of I/J provides an upper bound on λ0.
Here we assume that y(x) will be an even function, and so take y(x)=ax6 +bx4 +cx2 +d
for arbitrary parameters a, b, c, d, with a + b + c + d = 0. This condition leaves three degrees of freedom, but for Bessel’s problem we can also set d = 1, and be concerned only with ratios of the parameters.
Find explicit, analytic formulæ for I and J using sympy.
For a first estimate, take a = 0, and use b as the only remaining independent parameter. Evaluate I/J for a suitable range of values of b, and plot a graph. We seek the value of b which minimises this function, and you may find it by trial and error; or by writing an optimization code; or by using sympy.utilities.lambdify to access scipy.optimize.minimize. Plot the trial function which corresponds to this minimum, and compare, quantitatively to the exact solution.
For a second estimate, we re-introduce parameter a, and consider the two-parameter space of a and b. Draw a contour plot of I/J, and demonstrate that there is a minimum close to that of the first estimate, perhaps by overcoming multi-dimensional difficulties with using lambdify and minimize.
Draw a contour map for I/J in the range −8 ≤ a ≤ 0, 8 ≤ b ≤ 12. Locate a saddle point here, either by trial and error, or by some other method. Draw the graph of the trial function corresponding to the saddle point and comment on its relation to the exact solution to the original equation.
4
MATH 3001 Computational applied maths 2019/20
3 Fast Fourier transforms
In MATH2375 (or similar) you will have been introduced both to Fourier series and Fourier transforms. For a function x(t), the Fourier transform is defined as
which has inverse transform
+∞ −∞
1 +∞
x(t) = 2π
X(ω) =
x(t) exp (−iωt) dt
X(ω) exp (iωt) dω.
−∞
In the above, we assume that we have a function x(t) defined over the continuous domain t ∈ (−∞, ∞). But, for computational problems, it is simply not possible to store or work with infinite sets of data: in reality we would more likely have data obtained at (for example) a finite number of regularly spaced in- tervals in t. So, for example, we might have a set of n values of x, obtained at times t = 0,h,2h,3h,…
xm =x(t=mh), m=0,1,2,3,…,n−1. For such data, we can work with a discrete Fourier transform, defined as:
n−1 mk
Xk=xmexp −2πin , k=0,1,2,3,…,n−1
with inverse transform
m=0
1 n−1 mk
xm = Xkexp 2πi nn
, m=0,1,2,3,…,n−1. These transforms enable us to work with discretely defined sets of data.
It is important to notice that for k = n − j, exp 2πimk = exp −2πimj . Hence, for k > n/2, nn
it is better to think of the transform as representing negative frequencies given by index k′ = (k − n). This makes sense because, by definition, Xk (and in fact xm) are n-periodic. The frequency obtained from k = n/2 (for n even) is called the Nyquist frequency, and is the largest meaningful frequency appropriate for the discretely sampled data.
On the face of it, the time taken for a computer to perform the discrete Fourier transform appears to scale with number of points n as O(n2), since we need to calculate Xk for each of n values of k, and each of these requires us to sum over n values of m. Potentially this might take a long time if n is large. Fortunately a nice algorithm to perform this more quickly exists, the Fast Fourier Transform (FFT) [1]. This takes a time which scales as O(n log n) to perform the Fourier transform (i.e. much faster if n is large). An implementation of the FFT is available within the numpy Python package [2], and we can simply use it for calculations.
The Fast Fourier Transform is a hugely important algorithm, and is used especially in signal pro- cessing. But we can also use it as the basis for a spectral method for solution of PDEs. We will look at both applications in this project. Here are some suggested tasks, some of which are open ended.
1. First, make sure you are comfortable with the principles of the discrete Fourier transform and its implementation in numpy. Noting that t = mh, to what angular frequency ω does the kth component Xk correspond? [Hint: look at the quantities inside the exponential of the discrete and continuous Fourier transforms]. Consider taking the limits h → 0 and n → ∞ to convince yourself (at least roughly) that the continuous limit of the discrete transform gives the Fourier transform. You might find reading the first few pages of chapter 12 of [1] useful here, or find an alternative reference.
k=0
5
MATH 3001 Computational applied maths 2019/20
2. Next write yourself a short piece of code to convince yourself that the numpy FFT package is working as you expect it to do. You can find descriptions of the various routines at [2]. Your code should obtain the forward transform of a small set of data (e.g. x=[1,2,-1,4] where the elements of the list could be real numbers, or complex numbers). The code should then take the reverse transform of the result, to obtain the initial set of data. Check that the results of both the forward and reverse transforms are what you expect, by calculating the answer yourself with pen and paper! Investigate the difference between numpy.fft.fft (which takes the transform of a set of complex number data, and numpy.fft.rfft which assumes the input data are real numbers. Pay particular attention to the ordering of the output from the forward transform, and the frequencies that the output corresponds to. Why is there roughly half as much output when the input is real numbers? Note in particular the Nyquist frequency, and the order in which the data is output (see section 12.1 of [1] and [2]).
You could also take a look at the routines found at [2] for computing FFTs of 2D data, again paying close attention to the ordering and shape of the output, especially for the Fourier transform of real data.
3. Now you should be in a position to demonstrate filtering a set of data. Your task is to create a suitable set (or sets) of input data, and then demonstrate the following applications: (i) high pass filtering (in which the low frequency data is greatly reduced), (ii) low pass filtering (in which the high frequency data is greatly reduced) and (iii) removal of noise at a particular frequency from the data (this type of filtering might be important if, for example, data was corrupted by electrical mains noise at 50Hz). In each case, you should construct input data containing a large number of points (where n equals say 512 or 1024), take the Fast Fourier transform of the data, manipulate the transformed data appropriately (e.g. reduce the amplitude of some frequencies) and then transform the filtered data back.
[Extension: I don’t know how easy this will be, but can you do the same for 2D data using the 2D FFT? 2D data could be used to construct a 2D greyscale image, in which case a high-pass filter (for example) should smooth the image out. Can you perform this, and represent the results visually?]
4. It is also possible to use the FFT to compute extremely accurate solutions to PDEs, especially linear PDEs. Consider, for example, the heat equation:
∂f =K∂2f. (3.1)
∂t ∂x2
Wemayconsiderdiscretisingthespacevariablexintoasetofpointsx = mh,m = 0,1,2,3,…,n− 1 and then representing the function f (x, t) in terms of its discrete Fourier transform as:
1 n−1 mk
f(x,t) ≈ fm(t) = Fk(t)exp 2πi nn
,
k=0
1 n−1 xk
= Fk(t)exp 2πi . n hn
k=0
Then, the second derivative w.r.t. x is:
n−1
∂2f 1 4π2k2 xk ∂x2=n −h2n2Fk(t)exp 2πihn ,
k=0
and substitution into Eq. 3.1 allows a series of ODEs to be written for each of the Fourier com- ponents Fk(t), which can be solved easily if the initial value Fk(t = 0) is known (this process is analagous to the separation of variables method learnt in MATH2375). Hence, a solution to Eq. 3.1 for given initial condition can be obtained by the following method:
6
MATH 3001 Computational applied maths 2019/20
• Specify an initial condition fm(t = 0) and perform the FFT to obtain Fk(t = 0). • Evolve the obtained ODEs to obtain Fk(t) at whichever time(s) are desired.
• Perform the inverse FFT to obtain fm(t).
You should write a code to achieve this, for a suitable initial condition (e.g. a sharp spike in f in the middle of your domain which will spread out), displaying the result at various appropriate times t or (better) as an animation.
[Extensions: What happens if you put your sharp spike at the edge of the domain rather than in the middle? And, can you generalize the above to 2D?]
5. In principle, the same methodology can be used for other PDEs, such as the wave equation:
or the damped wave equation:
∂2f =c2∂2f ∂t2 ∂x2
∂2f +μ∂f =c2∂2f. ∂t2 ∂t ∂x2
In each case you would need to supply initial conditions on the rate of change, ∂f , as well as on ∂t
f itself. Can you write code to solve these problems, or a PDE of your own, and explore different initial conditions?
References
[1] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery Numerical Recipes, 2nd Edition, Cam- bridge University Press (1986). See chapters 12 and 13.
[2] https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html
7
MATH 3001 Computational applied maths 2019/20
4 How can we accurately evolve planetary orbits using geometric inte- grators?
Consider the gravitational two-body problem for a planet in orbit about a star of mass M
d2r =F(r)=−GMr, (4.1)
dt2 r3
where G is the gravitational constant, and r = rrˆ is the position vector of the planet relative to the star (and rˆ is the corresponding unit vector). This is approximately valid to describe the orbit of a single planet with a mass that is much smaller than that of the star.
We will investigate several different numerical methods for the solution of these equations to deter- mine which methods might be suitable for long-term orbit integrations. For one example, we require a method that can predict the evolution of the Solar System for millions to billions of years (in principle) without accumulating large errors that would make the method unusable.
We can write Eq. 4.1 as
r ̇=v, v ̇=F(r)=−GMrˆ, (4.2) r2
where v is the orbital velocity. We will consider the following integration schemes with timestep size h, which give the updated position (rn+1) and velocity of the planet (vn+1):
(i) Forward Euler: rn+1 = rn + hvn; vn+1 = vn + hF (rn).
(ii) Modified Euler: rn+1 = rn + hvn; vn+1 = vn + hF (rn+1).
(iii) Leapfrog: r′ = rn + hvn; vn+1 = vn + hF(r′); rn+1 = r′ + hvn+1. 22
(iv)Runge-Kutta(RK4):ifwedefinew=(r,v)T,sothatoursystemisw ̇ =f(w,t),thenthisscheme is given by: wn+1 = wn + 1 (k1 + 2k2 + 2k3 + k4),
6
where k1 = hf(wn,tn), k2 = hf(wn + 1k1,tn + 1h), k3 = hf(wn + 1k2,tn + 1h) and k4 = 22 22
hf(wn +k3,tn +h).
Euler methods are first-order; leapfrog is second-order; Runge-Kutta is fourth order.
(a) Begin by writing a code to follow the motion of the planet orbiting a star of mass M , with semi-
major axis a and eccentricity e (a measure of the ellipticity of the orbit when e ∈ [0, 1)), using method (i). You may assume that the motion is in the x − y plane and that the orbit starts from apocentre (the furthest distance from the planet to the star). The initial conditions are then x = a(1 + e), y = 0, vx =
0,vy =GM(1−e).TheperiodoftheorbitisP =2π a3 .YoumayassumeGM =a=1. a(1+e) GM
Evolve one whole orbit of the planet if e = 0.5, using 1000 force evaluations per orbit (i.e. h = P/1000). Plot its orbit on the x − y plane, over-plotting the analytical prediction for the shape of the orbit given by
x=a(cosE+e), y=bsinE (4.3)
√
whereb=a 1−e2 andE∈[0,2π].Markthepositionofthestarinyourplot.
Now modify your code so that it can use each of the remaining three integration schemes listed above, keeping the same number of force evaluations per orbit (which means h = P/250 for RK4), and produce a similar plot for each of these. What do you notice?
(b)TheenergyandangularmomentumaregivenbyE=1(v2+v2)−GM andL=xv −yv. 2xyryx
Show analytically that these are conserved for a Keplerian orbit (you may use sympy if you wish).
We will use each of the methods (i)–(iv) to follow the motion for N orbital periods, so that the end time is NP. For e = 0.5, 300 force evaluations per orbit and N = 102, plot the fractional error in the energy |E − E0 |/|E0 |, where E0 is the initial energy, as a function of time for all 4 integration schemes on a log-log scale (using different colours or symbols for each method). Then produce another plot showing the fractional error in the angular momentum |L − L0|/|L0|, where L0 is the initial angular momentum. Do any of these schemes conserve energy and/or angular momentum? Which method works best for
one orbit? Which method appears to work best for long-term orbit integrations? 8
MATH 3001 Computational applied maths 2019/20
Also, as in (a), plot your orbits on the x − y plane, showing all orbits from t = 0 to t = NP (this requires storing the positions for many or all time-steps). Does the deviation from the analytical orbit solution in each case agree with what you expect from considering the numerical evolution of the energy and angular momentum?
Now produce similar plots of the fractional error in the energy and angular momentum, and plots showing the orbit in the x − y plane for each integrator for e = 0.9, for N = 102 orbital periods and 1000 force evaluations per orbit. How does this differ from e = 0.5?
We will now explore in more detail some of the reasons why some methods are more suitable than others for long-term orbit integrations.
(c) We will investigate if any of these methods are time-reversible. Starting from the configuration in (b) for e = 0.5, integrate forward in time for N = 10 orbits, reverse all velocities, integrate backward in time for N = 10 orbits, then reverse the velocities again. Then compute the magnitude of the error between the initial and final positions and velocities. Do this for all 4 schemes. Which schemes are time-reversible (allowing for small round-off errors < 10−10) and which are not? Time-reversibility of a numerical scheme is related to the absence of numerical dissipation, so it might be considered a desirable property.
(d) Eqs. 4.1 can be derived from a Hamiltonian, so we refer to this system as a Hamiltonian system. This means that Liouville’s theorem applies, which states that the flow in phase space generated by a dynamical system (Eqs. 4.1) governed by a Hamiltonian conserves phase-space volume. The phase- space is the 4-dimensional space containing points (r, v). By considering analytically the numerical schemes (i) and (ii) (you may use sympy if you wish), show that one of them conserves phase space volume, but the other does not. To do this, you should consider the Jacobian determinant of the mapping from (rn, vn) to (rn+1, vn+1). (What about (iii)? sympy might be helpful here.)
A geometric integration algorithm is a numerical integration algorithm that exactly preserves some geometric property of the original set of differential equations. Their motivation is that preserving the phase-space (r, v) geometry of the flow (e.g. phase-space volume) determined by the real system is more important than minimising the one-step error, at least for long-term orbit integrations i.e. that higher- order is not necessarily better. Modified Euler and leapfrog are geometric methods but Runge-Kutta and forward Euler are not.
(e) Consider the N-body problem for N > 2 and consider investigating similar issues as those discussed above. You might also wish to explore some even more accurate numerical schemes e.g. Wisdom-Holman integrators (Wisdom & Holman, 1991, Astronomical Journal, 102, 1528) or IAS15 (https://arxiv.org/abs/1409.4779).
Some initial references
Goldstein, Poole & Safko, Classical Mechanics (3rd edition), Pearson/Addison Wesley, 2002. Murray & Dermott, Solar System Dynamics, Cambridge University Press, 1999.
9
MATH 3001 Computational applied maths 2019/20
5 Solving the heat equation
The conduction of heat down a solid metal bar of length L may be described by the one-dimensional diffusion equation (commonly referred to as the heat equation)
∂θ =K∂2θ, (5.1) ∂t ∂x2
where θ(x, t) is the temperature averaged over the cross-section at a distance x along the bar and K is the thermal diffusivity (a positive constant). This equation arises in many different problems in nature. It is appropriate here if there is negligible heat flux through the sides of the bar.
The aim of this investigation is to compute the solutions of this equation numerically using a finite- difference method, and to explore its performance by comparison with analytical solutions. The simplest case is where both ends of the bar are held at a fixed temperature, which can be taken to be zero without loss of generality i.e. θ(0, t) = θ(L, t) = 0 (why?).
(i) Derive the analytical solution for the evolution of an initial temperature distribution along the bar θ(x,0) = f(x),
∞ nπx nπ2 θ(x,t)= Ansin L exp −K L t ,
n=1
whereAn = 2 Lf(x)sinnπxdx.Whatisthesolutionwhenf(x)=sinπx?
(ii) We will now solve this problem numerically in a domain of unit length (L = 1) using a simple finite-difference method. You may also assume K = 1 (which can be thought of as re-scaling our time variable). We will represent the solution in the domain x ∈ [0, 1] using a discrete set of N points xi (or N − 1 intervals), such that the interval length is δx = 1/(N − 1). We also approximate the time derivative using a first-order forward difference scheme:
(5.2)
L0LL
∂θ(x, t) = θ(x, t + δt) − θ(x, t) + O(δt), ∂t δt
and the second-order space derivative by the central difference scheme (at the current time)
∂2θ(x,t) = θ(x+δx,t)−2θ(x,t)+θ(x−δx,t) +O((δx)2). ∂x2 (δx)2
This gives the numerical scheme
θm+1=θm+Kδtθm −2θm+θm , i i (δx)2 i+1 i i−1
where θm is an approximation to θ(iδx, mδt). Define the Courant number, C = Kδt . i (δx)2
(5.3)
(5.4)
(5.5)
Write a program to implement this numerical scheme with initial condition θ(x, 0) = f (x) = sin (πx). For the case N = 5 and C = 1 , produce a plot that compares the analytical and numerical solutions,
3
then tabulate the values of both of these and the magnitude of the error at each point, at t = 0.25 and 0.75.
We may define the root-mean-square error E = 1 |θ(x, t) − θan(x, t)|2dx, where θan(x, t) is 0
the analytical solution from (i). This can be used to measure the accuracy of the numerical scheme
for different values of N and C. Run your code for several different values of N and C (assuming
C < 1 ) until t = 0.5, and produce separate plots showing how E varies with N and C . Are your results 2
consistent with the theoretical order of accuracy of the scheme in space and time?
(iii) Run your code with N = 10 and C = 2 until t = 1. Plot the numerical and analytical solutions
3
as a function of x. What happens in this case?
10
MATH 3001 Computational applied maths 2019/20
In order to explain this strange behaviour it is helpful to perform a Von-Neumann stability analysis. This allows us to determine the stability of the numerical scheme i.e. whether or not the error rapidly grows with time to dominate the solution. Since Eq. 5.5 is linear and has constant coefficients, we may seek independent solutions (where we have now used an index j instead of i to avoid confusion with the imaginary number i)
θ(jδx, mδt) = ξmeikjδx, (5.6)
where k is a real spatial wavenumber and ξ(k) is an amplitude which depends on k, which is raised to the power m in this expression. Substitute this into Eq. 5.5 and divide both sides by θ(jδx,mδt) to obtain an equation for ξ = θ(jδx, (m + 1)δt)/θ(jδx, mδt), which can be thought of an as amplification factor. Hence, show that
2 212
|ξ| = 1−4Csin (2kδx) . (5.7)
A numerical scheme is said to be stable if |ξ|2 ≤ 1 for all values of kδx (why?). Determine a condition
on C for this numerical scheme to be stable. Can you now explain the strange behaviour that you
obtained when C = 2 ? 3
(iv) We will now implement a more stable and accurate numerical time-integration scheme that is ideally suited to solving the heat equation: an implicit Crank-Nicolson scheme, which approximates the time-derivative using
θm+1 − θm K ∂2θm+1 ∂2θm
i i= 2 + 2 , (5.8)
δt 2 ∂x i ∂x i
where the spatial derivatives are approximated using Eq. 5.4. This is an implicit scheme, since the
solution at time level m + 1 depends on the solution at the current time level (m + 1) in addition to the
previous one (m), so we must solve a matrix equation to find θm+1. By contrast, the method in (ii) is i
an explicit scheme, in which the the solution at time level m + 1 only depends on the solution at the previous time level (m).
For N = 5, we can express the values of the interior points (those not at the boundaries, for which θ0m =θ4m =0)forθattimelevelmasthevector
θ 1m
θm=θ2m , (5.9)
θ3m
and we can express the spatial derivatives in Eq. 5.4 using a derivative matrix, which for the case N = 5
can be defined by
0 1 −2
(Why is this a 3x3 rather than a 5x5 matrix?) Hence, we can express Eq. 5.8 in the form
Aθm+1 = b(θm),
−2 1 0
D2=1 1−21.
(5.10)
(5.11)
(δx)2
where A is a matrix (to be determined) and b(θm) is a vector that depends on the values at the previous time level. This can be solved to find θm+1 using the routine numpy.linalg.solve, and you may find it helpful to use the identity matrix numpy.eye to construct A.
Implement this numerical scheme. By computing E for various values of C (with N = 5), show that this scheme gives second-order convergence in time (i.e. that E ∝ C2), and is therefore more accurate than the scheme used in (ii).
11
MATH 3001 Computational applied maths 2019/20
GeneraliseyourcodetocopewiththecaseN =10,andrunitwithN =10andC = 2 untilt=1. 3
Plot the numerical and analytical solutions as a function of x. How does this behave in comparison with the scheme in (ii) for these parameters?
(v) In this example we have solved a linear PDE which has a simple analytical solution. But this is not typically the case for the nonlinear PDEs that often arise in reality, where we usually have to resort numerical methods such as those described here. Try any of the numerical schemes above on a PDE (or several PDEs) of your choice, which might not have an analytical solution.
Some initial references
G D Smith, Numerical solution of partial differential equations: finite-difference methods (3rd edition), Oxford University Press, 1985.
W. F. Ames, Numerical Methods for Partial Differential Equations (2nd edition), New York : Academic Press, 1977.
12
MATH 3001 Computational applied maths 2019/20
6 Using random numbers to integrate Introduction to Monte Carlo methods
Most computer programming languages have a command which allows the user to generate random, or pseudorandom, numbers. For example, the “Python” language has a set of standard commands for generating random numbers (either integers or real numbers) from a wide range of standard distributions, such as the uniform and normal distributions [1]. Random numbers can be used for statistical sampling (also called “Monte Carlo sampling”). A common application of this is in numerical integration, in which the function to be integrated is sampled at random points within the domain of integration, in order to estimate the value of the integral [2, 3]. For integration in one dimension, other numerical integration methods (such as Simpson’s rule) are usually more efficient than Monte Carlo Methods. But for multidimensional integrals, Monte Carlo methods become competitive, and for very large numbers of dimensions these are the only practical method for numerical integration.
There are a number of Starter exercises suggested in this document, and several open ended extension tasks. You are not necessarily expected to complete all extension tasks!
Starter task 1:
• Therandompackageinpython(seehttps://docs.python.org/2/library/random.html) contains several routines to generate random numbers, e.g. uniform distribution, normal distribu-
tion, exponential distribution. Write a Python code to evaluate the mean and variance of a random
sample drawn from some probability distribution(s). Then adapt the code to investigate different parameters and/or distributions. Check that the results make sense to you. If they do not, it may indicate a bug in the program. Edit the program until you understand what is happening.
Monte Carlo Integration
A good reference for this is chapter 7.6 of Ref. [4].
A common task in mathematics is the evaluation of integrals. Although in the problems set in under-
graduate mathematics courses integrals can be performed analytically, this is rarely true in applications. Consequently we need numerical methods to approximate integrals as sums, so that
b a
( b − a ) N
N f(xi). (6.1)
i=1
f(x)dx ≈
Furthermore we can also estimate the error, from the mean squared error as
b N f(x)dx ≈
wif(xi),
a
i=1
where we take a weighted sum of the values of the function at N points xi. In traditional quadrature schemes such as the trapezium and Simpson’s rules [4] the positions of the points are predetermined, e.g. equally spaced points.
However, in Monte Carlo integration the points are chosen randomly from a distribution. Since the mean value of f can be defined as
1b
⟨f⟩ = b − a
we can determine the value of the integral by sampling f to find its mean value. Thus if we draw N
points,x1,...xN fromauniformdistributionintheinterval[a,b]thenwecanestimatetheintegralas,
a
f(x)dx,
σ
(b−a)√ (6.2)
N
13
MATH 3001 Computational applied maths 2019/20
where σ is the standard derivation of f over the interval. This means that the accuracy of the estimate improves as the square root of the number of points sampled.
Starter task 2: Obtain a formula for σ and hence an error estimate. How does this compare with quadrature schemes such as the trapezium rule or Simpson’s rule? Try estimating 1 2 cos2(x) dx using
20
Monte-Carlo sampling. Compare the answer with the trapezium and/or Simpson’s rule.
Multi-dimensional Integration
The above example shows that there are better ways of integrating in one-dimension than Monte-Carlo (MC) methods. However, we shall discover that in multi-dimensional integration MC methods have their advantages.
Suppose we now consider the integration of the function f (x) over a D-dimensional domain x ∈ V. As shown in the vector calculus course this can be computed via nested one-dimensional integrals as
b1 bD
f(x)dV = ··· f(x)dxD ...dx1
V a1 aD
where V is contained in the hypercuboid (a, b) and we define f to be zero for points outside V . If we use a quadrature scheme with n points for each of the integrals we require a total of N = nD function evaluations. Now unless V fits precisely to the boundaries of the domain, f will be discontinuous. This means the error in each quadrature integral will scale as 1/n (since the higher order convergence of the trapezium and Simpson’s rules depend on function having continuous derivatives). Thus the total error will scale as N−1/D. So when D is large the convergence is very slow. Even where we have a rectangular domain over which the function is differentiable the convergence rate we decrease as the number of dimensions increases. This is the “curse of dimensionality”.
Now suppose we use Monte-Carlo sampling. The integral can be estimated as
VN σ
f(x)dV ≈ N f(xi) ± V √ . (6.3)
V i=1 N Hence the error scales as N −1/2 irrespective of the number of dimensions.
Starter task 3:
As an example let us consider calculating the area of the unit circle by integrating the function
1 x21+x2<1
f=
0 otherwise
over the domain V = {(x1, x2) : |x1| < 1, |x2| < 1}. The answer is of course π. Write and execute your own program to estimate the integral using Monte-Carlo sampling. By plotting a graph, examine how the error changes with N. Run repeated calculations for particular values of N and compare the variance with the error estimator. How many samples would you need to calculate π to 3, 5 and 10 decimal places?
Extension task i:
Modify your programme to calculate the volume of a unit sphere and higher dimensional hyper- spheres. The analytical formulas for the answers are given at http://thinkzone.wlonk.com/MathFun/Dimens.htm. How does increasing the dimension modify the convergence?
14
MATH 3001 Computational applied maths 2019/20
Extension task ii:
Examine what happens if you replace f with the function
rα r<1
f=
0 otherwise
where r = x21 + x2. In particular look at what happens for the cases −2 < α < 0 when the integrand is singular at (0, 0) but integral remains convergent. Also see what happens at α = −2 when the integral has a logarithmic divergence.
Extension task iii:
So far we have performed integrations using random numbers drawn from uniform distributions. But, supposing we draw our N points, x1, . . . xN from a probability distribution p(x), then:
( b − a ) N
f(x)p(x)dx ≈ N
• Integration over infinite domains. For example, to perform the integral
∞
f (x) exp(−x/λ)dx
0
we can sample from the exponential probability distribution p(x) = exp(−x/λ)/λ. Can you demonstrate this working in practice?
• Variance reduction. By choosing good probability distributions p(x) which concentrate the ran- dompointsx1,...xN intheregionwheretheintegrandoff(x)p(x)dxislarge,thentheerror in the Monte Carlo integration is reduced. Can you prove this is true, and demonstrate the effect with a suitably chosen example? [Hint: consider the variance in the quantity
( b − a ) N
⟨f⟩ =
This can then be used to evaluate the integral f (x)p(x)dx. This allows:
f(xi)
when the N points, x1, . . . xN a drawn from a probability distribution p(x).
References
[1] https://docs.python.org/2.7/library/random.html
[2] Malvin H. Kalos, Paula A. Whitlock, “Monte Carlo methods” 2nd edition (Wiley-Blackwell,
2008.)
[3] William L. Dunn and J. Kenneth Shultis, “Exploring Monte Carlo methods” (Elsevier/Academic Press, Amsterdam and London, 2012.)
[4] William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery, “Numerical Recipes”, In various editions, the books have been in print since 1986. Cambridge University. You can also look at a limited number of pages online.
f(xi).
i=1
N
i=1
15
MATH 3001 Computational applied maths 2019/20
7 Determining the period of a pendulum Background
A frictionless pendulum supported by a massless rigid rod of length l moves under the influence of a downward gravitational acceleration g. At time τ, the angle from the downward vertical is denoted by θ, and can be shown to satisfy d2θ/dτ2 = −(g/l)sinθ. Introducing a rescaled nondimensional time t = τ/k, for a suitably chosen constant k with units of time, this governing equation can be rewritten as
θ ̈= −sinθ, θ(0) = a, θ ̇(0) = 0, (7.1) where ̇denotes differentiation in t, and the initial conditions correspond to a release from rest at an angle
a (which satisifies 0 < a < π, without loss of generality). • Task 1: Derive the results stated above.
Despite the apparent simplicity of the system (7.1), there is no (known) formula for the period T of the resulting oscillation in terms of the release angle a. Here we will construct the dependence of T upon a using numerical computations. Our numerical computations should recover three known results:
• when a ≪ 1, there is a small-amplitude approximation for the equation of motion (7.1), which can be solved exactly (see below). In this limit, there is a remarkable property: the period of oscillation T is always 2π, regardless of the release angle a, i.e.,
T ≈2πwhena≪1.
• Some particular values of T are known; for example,
(Γ(1/4))2 π T = √π ≈7.4163whena= 2,
where Γ is the gamma function, which can be evaluated using scipy.special.gamma.
• It can be argued on physical grounds that
T →∞asa→π.
The small-amplitude approximation
When a ≪ 1, it follows that |θ| ≪ 1 for all t, so sinθ ≈ θ, and (7.1) becomes θ ̈= −θ, θ(0) = a, θ ̇(0) = 0.
• Task 2: show that the standard solution method for solving the linear ODE (7.5) leads to (7.2).
(7.2)
(7.3)
(7.4)
(7.5)
The period can also be derived using an alternative integral method. Multiplying (7.5) by θ ̇ and integrating, for the initial downwards swing we obtain θ ̇ = −√a2 − θ2, and integrating over 0 < t <
T/4 it follows that
a dθ
T=4 √2 2, (7.6)
0a−θ
i.e., that T can be expressed in integral form.
• Task 3: flesh out the derivation of (7.6), evaluate T from (7.6) using an appropriate substitution, and hence recover (7.2).
16
MATH 3001 Computational applied maths 2019/20
The extension to arbitrary amplitude
Now return to the unapproximated system (7.1), which cannot be solved in terms of simple functions. • Task 4: Derive the relevant analogue of (7.6), which is
√a dθ
√cosθ−cosa. (7.7)
Unfortunately, this integral cannot be evaluated (exactly) using a simple substitution. We thus resort to a numerical integration (or quadrature) scheme. However, this is not as straightforwards as it seems, since the integrand becomes unbounded as θ → a−.
• Task 5: We will use the trapezoidal rule (TR) and Simpson’s rule (SR) for the numerical integration. Before applying these to (7.7), test these schemes on a simple integral, such as
π/4
I= cosθdθ. (7.8)
0
Divide the interval [0, π/4] into n intervals, each of width ∆θ = (π/4)/n. By successively doubling n (starting at 100, say, and calculating for about 4 other values), calculate a sequence of approximations In for each method. These should converge towards I; quantify the rate of convergence by calculating the errors En = |In − I|. You should find that En ∝ n−p, where p > 0 can be found by noting that log En should be linear in log n. You should recover the theoretically predicted values of p = 2 for TR and p = 4 for SR.
• Task 6: Now return to (7.7) with a = π/2, for which the exact value of T is known (via (7.3)). Divide the integration interval [0, π/2] into n intervals, each of width ∆θ = (π/2)/n. Evaluate (7.7) using a ‘hope for the best approach’: apply either TR or SR, but simply neglect the end-point at θ = a where the integrand becomes unbounded. By successively doubling n (starting at 100, say, and calculating for about 4 other values), calculate a sequence of approximations Tn for each method. From the implied errors En = |Tn − T |, calculate the rate of convergence for each of TR and SR. You should find that the rate of convergence is disappointingly slow in both cases.
• Task 7: A more sophisticated technique is to perform the improper integral associated with the singular endpoint exactly, and to integrate the remainder using a standard numerical method. By expanding cos θ as a Taylor series about θ = a, show that
1 1 √ 3/2 − √cosθ−cosa=(sina)(a−θ)+b a−θ+c(a−θ) +··· as θ→a , (7.9)
for some constants b and c, which can be determined if you wish. Hence show that (7.7) may be usefully rewritten as
T=2 2
0
2a √a 1 1 T =4 sina+2 2 √cosθ−cosa−
explaining why the remaining integral should be well-behaved. Calculate a sequence of approximations for Tn based upon (7.10), again taking a = π/2 for each of TR and SR. By considering the scaling of the corresponding errors En with n, determine the rate of convergence for each of TR and SR. What do you notice about the rates of convergence relative to those from tasks 5 and 6? Is there a way to improve the rate of convergence for SR towards its theoretical maximum?
• Task 8: Given a sequence of approximations Tn with a known rate of convergence p, we can write Tn = T + k , when n ≫ 1,
np 17
0 (sin a)(a − θ)
dθ, (7.10)
MATH 3001 Computational applied maths 2019/20
for some constant k. So we can estimate T by making a linear fit of Tn against n−p. By doing this at a = π/2, you should find that you can then estimate T to a high degree of accuracy (meaning within 10−10 for SR).
• Task 9: Use one of your methods (e.g., TR with singularity removed) to calculate T (a) as a varies over 0 < a < π. You should produce (i) a table giving T when a = π/6, π/3, π/2, 2π/3, 5π/6, (ii) a graph showing T (with a finer set of a), (iii) a corresponding graph showing T − 2π. To what accuracy (i.e., how many decimal places) do you think your results are given?
Extensions
• How much ‘better’ is the method based upon (7.10) relative to that based upon (7.7)? One way to measure ‘better’ is through efficiency, or speed of computation. How long does each scheme take to calculate T to a specified accuracy (e.g., 8 decimal places)? How does the computation time depend upon the specified accuracy? For example, does a calculation correct to 10 decimal places take twice as long as one to 5 decimal places?
• The weakly nonlinear regime (i.e., a small but not vanishingly small) can be analysed using per- turbation theory, as is taught in level 3 (or 5) Mathematical Methods. You could try either or all of the approaches below, each of which lead to the approximation
a2
T≈2π 1+16 . (7.11)
– Writeθ=aθˆ1(t)+a2θˆ2(t)+···,useaTaylorseriesforsinθin(7.1),gathertogetherterms of like order in a, and solve the resulting differential equation and boundary conditions. By calculating terms up to θˆ3(t) and then finding the first zero crossing in θ, which occurs at t = T/4, derive (7.11).
– Use the method of multiple scales. Introduce a slow time scale t ̃ = a2t, and write θ = ˆ ̃3ˆ ̃ ˆ
aθ1(t, t) + a θ3(t, t) + · · · . Eliminating the secularity for θ3, you will find a leading-order ˆ ̃ ̃
solution θ1 = cos(t/16) cos t + sin(t/16) sin t, from which (7.11) follows.
– UseTaylorseriesexpansionsforcosθandcosain(7.6),andevaluatetheintegralsthatresult
at O(1) and O(a2), perhaps using a standard trigonometrical substitution.
If you are determined, higher-order approximations (i.e., the corrections in T of O(a4)) can be
derived using these methods.
• Canyouavoidtheproblemsinevaluatingtheintegral(7.7),bymakingaprofitablesubstitution?In particular, can you obtain an equivalent integral without a singular integrand, so that the numerical evaluation becomes straightforward?
• VariousapproximationsforT(a)arederivedintheresearchpaperliterature.Youcoulddocument some of these, and perhaps add the approximations to your results from Task 9 (e.g., additional columns for the table, or dashed lines for your graphs).
• Investigatecompleteellipticintegralsofthefirstkindtofindoutwheretheexactresult(7.3)comes from. Are there any other values of a for which T can be determined exactly in a similar way?
• Feel free to implement any other methods for calculating T that seem sensible to you.... There are many different approaches.
18