ACM41000: Uncertainty Quantification
Dr James Herterich
Assignment 3
This is a graded assignment. The due date is April 22 at 9pm. Upload your
solutions to Brightspace.
Please upload a single PDF file. This is to include scanned written
solutions and PDF outputs of your R code and figures. Please
combine the PDFs into one.
For full marks, ensure your work is clear and well-commented.
1. Why do zebras have stripes and cheetahs have spots? A Turing pattern is a
concept introduced by English mathematician Alan Turing in a 1952 that
describes how patterns that emerge in nature, such as stripes and spots,
may arise naturally and autonomously from a homogeneous, uniform state.
It is often described in coupled models of reacting and diffusing chemicals.
Consider such a model below for the dynamics of X1(x, t) and X2(x, t)
∂X1
∂t
= D1
∂2X1
∂x2
+ 1 +X21X2 − 4X1
∂X2
∂t
= D2
∂2X2
∂x2
+ 3X1 −X21X2
(a) Solve the coupled system numerically in R using the ReacTran package.
Solve in x ∈ [0, 1] and t ∈ [0, 10]. Take the diffusion coefficients to be
D1 = D2 = 0.02,
boundary conditions to be
X1(0, t) = X1(1, t) = 1 X2(0, t) = X2(1, t) = 3
and initial conditions
X1(x, 0) = 1 + sin(2pix) X2(x, 0) = 3
The initial setup is provided to you in an accompanying file on Brightspace.
(b) Derive the sensitivity equations for X1 and X2 with respect to the
diffusivities D1 and D2:
s11 =
∂X1
∂D1
s21 =
∂X2
∂D1
s12 =
∂X1
∂D2
s22 =
∂X2
∂D2
(c) Solve the sensitivity equations taking initial conditions of all sij to be
zero and boundary conditions on sij also to be zero.
(d) Briefly summarise which variables are most sensitive to which param-
eters.
Page 1 of 2
2. Based on the research article,
Goulet, Dave. “Modeling, simulating, and parameter fitting of
biochemical kinetic experiments.” SIAM Review 58.2 (2016):
331-353
you will recreate calculations using methods learned during the module.
Complete the following tasks:
(a) Give a brief introduction to the model described in (3.10), describing
solution variables, model parameters, and initial conditions.
(b) Write a program in R to solve the model equations (3.10) numerically.
Note that the initial conditions are part of the right-hand side of the
ODEs requiring some addional programming thought on your end.
Use k+ = 200 and k− = 2 as parameter values.
The model equations form a stiff system, meaning instead
of calling ode() to solve the system of differential equations,
you must call lsode(). Otherwise your machine may freeze during
computations!!!
(c) Fit the data in Table 2 to the solution of equations (3.10a-c) to find the
best fit k+ and k− according to a Levenberg–Marquardt least-squares
error algorithm. Note a simpler direct version of the sum of square
errors is appropriate here:
(xout-xdata)2 + (yout-ydata)2 + (zout-zdata)2
where you define the above based on your program. In order to solve
the system for two parameters (k+, k−), take the time interval t ∈
[1/30, 20] and use the provided data at t = 1/30 for initial conditions.
(d) Briefly discuss your implementation and results.
(e) Derive the sensitivity equations for the model (3.10), based on param-
eters k+ and k−, using the procedure shown in the lectures. Compare
with (B.1a-1f) in the article. Use either your derived sensitivity equa-
tions, and/or (B.1a-1f) quoted in the article, to discuss and justify any
difference in the best-fit parameters you obtain from that quoted in
the article for the Table 2 dataset.
Note: You do not need to replicate the exact results using the
particular methods as stated in the article. You are using the
model and dataset provided, but performing calculations with
methods you have learned. Some of the algorithms used in the
article are slighlty different than what we use.
Page 2 of 2