COSC2500/COSC7500—Semester 2, 2024
Exercises—due 3:00 pm Friday 25th October
Required
The grade of 1–5 awarded for these exercises is based on the required exercises. If all of the required exercises are completed correctly, a grade of 5 will be obtained.
R4.1 Do either part (a), (b), or (c) below. (If you want to do more than one of these, see the additional exercises!) These are Q1 in Sauer computer prob- lems 8.1, 8.2, and 8.3.
a) From Sauer 8.1: Solve ut = 2uxx for 0 ≤ x ≤ 1, 0 ≤ t ≤ 1, for the sets of boundary conditions
i. u (x, 0) = 2 cosh x for 0 ≤ x ≤ 1
u (0, t) = 2 exp(2t) for 0 ≤ t ≤ 1
u (1, t) = (exp(2) + 1) exp(2t − 1) for 0 ≤ t ≤ 1
(Solution is exp(2t + x) + exp(2t − x))
ii. u (x, 0) = exp x for 0 ≤ x ≤ 1
u (0, t) = exp(2t) for 0 ≤ t ≤ 1
u (1, t) = exp(2t + 1) for 0 ≤ t ≤ 1
(Solution is exp(2t + x))
using the forward difference method for step sizes h = 0.1 and k = 0.002. Plot the approximate solution (the mesh command might be useful). What happens if you use k > 0.003? Compare with the exact solutions.
HINT: You can use Program 8.1 (heatfd.m) from Sauer.
b) From Sauer 8.2: Solve the following initial–boundary value problems using the finite difference method with h = 0.05 and k = h/c. Plot the solutions.
i. utt = 16u xx
u (x, 0) = sin(πx) for 0 ≤ x ≤ 1
ut(x, 0) = 0 for 0 ≤ x ≤ 1
u (0, t) = 0 for 0 ≤ t ≤ 1
u (1, t) = 0 for 0 ≤ t ≤ 1
(Solution is u (x, t) = sin(πx) cos (4πt))
ii. utt = 4u xx
u (x, 0) = exp(−x) for 0 ≤ x ≤ 1
ut(x, 0) = −2 exp(−x) for 0 ≤ x ≤ 1
u (0, t) = exp(−2t) for 0 ≤ t ≤ 1
u (1, t) = exp(−1 − 2t) for 0 ≤ t ≤ 1
(Solution is u (x, t) = exp(−x − 2t))
c) From Sauer 8.3: Solve the Laplace equation for the following boundary conditions using the finite difference method with h = k = 0.1. Plot the solutions.
i. u (x, 0) = sin(πx) for 0 ≤ x ≤ 1
u (x, 1) = exp(−π) sin(πx) for 0 ≤ x ≤ 1
u (0, y) = 0 for 0 ≤ y ≤ 1
u (1, y) = 0 for 0 ≤ y ≤ 1
(Solution is u (x, y) = exp(−πy) sin(πx))
ii. u (x, 0) = 0 for 0 ≤ x ≤ 1
u (x, 1) = 0 for 0 ≤ x ≤ 1
u (0, y) = 0 for 0 ≤ y ≤ 1
u (1, y) = sinh π sin(πy) for 0 ≤ y ≤ 1
(Solution is u (x, y) = sinh(πx) sin(πy))
HINT: You can use Program 8.5 (poisson.m) from Sauer.
R4.2 Use a Monte Carlo method to find the area of a circle, the volume of a sphere, and so on, for higher-dimensional “circles” and “spheres” . Make a statistical estimate of the accuracy of your result, and compare with known results.
R4.3 From computer problems 9.3 in Sauer: In a biased random walk, the proba- bility of going up one unit is 0 < p < 1, and the probability of going down one unit is q = 1 − p. Design a Monte Carlo simulation with n = 10000 to find the probability that the biased random walk with p = 0.7 on the intervals
a) [−2, 5]
b) [−5, 3]
c) [−8, 3]
reaches the top. Calculate the error by comparing with the correct answer
(q/p)
a+b − 1/(q/p)
b − 1
for p ≠ q and interval [−b, a].
Additional
Attempts at these exercises can earn additional marks, but will not count towards the grade of 1–5 for the exercises. Completing all of these exercises does not mean that 6 marks will be obtained—the marks depend on the quality of the answers. It is possible to earn all 6 marks without completing all of these additional exercises.
A4.1 Do the other part (i.e., (a) or (b)) for R6.2.
A4.2 Find a PDE solver. Discuss the code, if source code is available. Test, and discuss.
A4.3 From computer problems 9.4 in Sauer: Use the Euler-Maruyama method to solve
dy = Btdt + (9y2 )1/3 dBt
over the interval t = [0, 1] with initial condition y(0) = 0. Note that Bn = Σ0 ∆Bi. Use step sizes of h = 0.1, 0.01, 0.001, and use 5000 simulations for each step size to find the mean value of y(1) and the error. How does the error vary with step size?
A4.4 Implement and test a random number generator.
A4.5 Generate a distribution of random numbers other than uniform. or normal. What kind of problem would this distribution be useful for?
A4.6 Find or implement code for high-dimensional optimisation. Test, and dis- cuss.
A4.7 Implement Monte Carlo integration for circles, spheres, etc., using C or some other compiled language. Compare the performance of your code with your Matlab version. Investigate the effect of using single precision instead of double precision.
A4.8 Write a popular science article describing a problem and the use of Monte Carlo methods to solve it.
Programming hints
Note that rand() and randn() can return a vector or matrix of random numbers. Instead of getting one random number at a time in a loop, you can get them all at once!
• rand(N) andrandn(N) will return N × N matrices of random numbers.
• rand(M,N) andrandn(M,N) will return M × N matrices of random num- bers.
R4.2 You need to demonstrate a sufficient understanding of the question. This means that you will need to do the integration for multiple dimensions. Going from a 1-sphere (asphere in 1 dimension) to a 10-sphere is sufficient (That is: 1, 2, 3, 4, 5, 6, 7, 8, 9 and 10).
If you are using rand() in MATLAB or similar, one way to do this is to sam- ple points inside a box. This method is finding the ratio of In/N, where In is how many random points lie inside the sphere, and N ishow many points used in the simulation. It is important to note that this is a ratio of the “box” it covers; this means you will need to multiply the ratio by the volume of the box.
Algorithm R4.3 (For a unit sphere).
Draw X from U(0,1)
Find the distance the point is from the centre . Count how many points are inside the sphere .
Calculate ratio = inside/N .
Calculate volume = ratio*boxVolume
R4.3 There are several ways to go about this question. When the walk reaches the top or the bottom of the interval (that is, if at any point the value of the walk equals a or −b) you break the calculation immediately and go to the next calculation.
The stating point for the walk is 0.
For some, it maybe helpful to think of it as a Discrete Markov Chain. Each state is finite, and once it reaches the boundary, it can keep going back to boundary with a probability of 1.
Algorithm R4.3 Easy algorithm, ignoring Murphy’s Law.
Initialse ReachedTheTop
Do N times
Initialise W
* Random walk process *
Stop
* Random Walk process *
while not at the boundary
Draw u from U
if u > q
go up
otherwise
go down
if W is at a
Increment ReachedTheTop
stop
if W is at -b
stop
There is a more advanced algorithm. This is only useful if you are going to vectorise; if you are not going to vectorise, do the one above.
Initialise W
Do M times
Draw dW from sign ( U - q )
W = W + dW
if Wi = a , set Wi to be inf
if Wi = -b , set Wi to be - inf
Stop
Sum up W equal to inf
Question: Why does this algorithm work? Is it mathematically correct? What does W(W==a) do?
A4.4 Bt is a partial sum, and is related to every random number drawn before (for that walk).
Algorithm R4.4 Non-vectorised.
Calculate the number of steps .
Do N times
*Euler -Maruyama step*
Stop
Find the mean .
*Euler -Maruyama step*
Do number of steps times Draw Z from N(0,1) . Calculate dB and Bt .
Apply Euler -Maruyama scheme
Stop
And of course, a more advanced vectorised version. (Thomas wrote these hints, and he vectorises a lot.)
Draw Z from N(0,1)
Calculate Bt and initialise Y
Calculate the number of steps
Do number of steps times
Apply Euler Maruyama - scheme
(using element -wise operations)
Stop
Find the mean .