MATH64071 Introduction to Uncertainty Quantification
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
MATH64071
Introduction to Uncertainty Quantification
Coursework 2021/22
1. Monte Carlo The Van der Pol oscillator was originally proposed by the Dutch electrical engineer and physicist Balthasar Van der Pol found stable oscillations in electrical circuits employing vacuum tubes1. A variant, the Duffing-Van der Pol oscillator, is given by the second order ODE for x(t):
= αx + βx˙ - x3 - x2 x˙ , x(0) = x˙ (0) = 0.1,
where α, β > 0.
(a) Firstly, show that this system can be rewritten as the following system of two
x˙ = y(t), (1)
y˙ = f (x(t), y(t)), x(0) = y(0) = 0.1,
where you must determine the function f (x, y). Derive the forward Euler ap- proximation of the resulting IVP (1), and implement it on a computer. Use this implementation to approximate the solution of this ODE with α = 1, β = 1, up to time T = 20, for the time steps ∆t = 10 − 1 , 10 −2 , 10 −3 , 10 −4 , and plot the re- sults. In the next question, we will be using this approximation within a Monte Carlo method for estimating the value of E(x(20)) and E(y(20)). Using your plots as an indicator for the degree of accuracy of the Euler approximation for this system, pick a value of ∆t that would be suitable for such a Monte Carlo estimator. Justify your answer by taking into consideration both the error and computational cost.
1 Reference: Wikipedia page on the Van der Pol oscillator
Simon L. Cotter: MATH64071. Coursework 2
(b) Suppose that α = 1, but that we are unsure of the value of β. We model this
uncertainty through a probability distribution,
β ~ U([0.75, 1.25])
Implement a Monte Carlo algorithm in order to estimate E(x(20)) and E(y(20)) = E(x˙ (20)), using the Euler approximation from part (a), and the chosen value of ∆t from part (a). Make appropriate plots to visualise the joint and marginal dis- tributions of the random variables x(20) and y(20) estimated using the Monte Carlo samples. You may use the m-file “hist2d.m” to create your visualisa- tions of the joint distribution, which can be found in the coursework folder on Blackboard2. How many samples would you need to have a standard error less than 10 −3 for your estimate of E(x(20))? Likewise for E(y(20))? Justify your answers.
[5 marks]
2. Numerical approximation of SDEs. An ODE is often not an appropriate type of model for a system which is inherently stochastic. We consider the stochastic version of the Duffing-Van der Pol oscillator:
dXt = Yt dt
dYt = /αXt + βYt - Xt(3) - Xt(2)Yt、dt + σ1Xt dBt(1) + σ2 dBt(2), Xt=0 = Yt=0 = 0.1,
where Bt(i) for i e {1, 2} are independent standard Brownian motions. In this example we will consider the case α = β = σ 1 = σ2 = 1.
(a) Derive the Euler-Maruyama method for the approximation of paths of the so-
lution of this SDE. Implement this method in any computer language you wish. Pick a time step ∆t and justify your choice. Simulate and plot 5 approximate realisations of this SDE, up to time T = 20.
[5 marks]
(b) Implement a Monte Carlo sampler to estimate the values of E(XT ) and E(YT ) at
T = 20. Make appropriate plots to visualise the marginal and joint distributions for the random variables XT and YT at T = 20. How many samples would you need to have a standard error less than 10 −3 for the estimators of each of these expectations? Justify your answer.
[5 marks] [TOTAL: 20 marks]
2 Usage: hist2d(X,Y,N,‘pdf’), where X and Y are arrays containing your samples, and N is the number of bins in each direction (I recommend 100, but make sure you have enough samples to make the results relatively smooth).
2021-12-07