Description
1. Lorenz equations and chaos. (20 marks)
You are asked to implement the Ralston 2-stage Runge-Kutta method for y
0
(x) = f(x, y(x)),
given by
k1 = f(xn, yn),
k2 = f
xn +
2
3
h, yn +
2
3
h k1
,
yn+1 = yn + h
1
4
k1 +
3
4
k2
,
1
for the numerical solution of ODE systems
U
0
(x) = F(x, U(x)),
and apply it to simulate chaotic solutions of the nonlinear Lorenz equation system.
Download the files driver1.m, FESystem.m and EqSystemTest.m from the class website.
These files provide details about the required input and output variables for the routines
and an example implementation of the Euler method (which is also called the Forward Euler
method, or FE), and will allow you to test your implementation.
Study the implementation of the system Forward Euler method that is provided in Matlab
m-file FESystem.m, with function heading
function U=FESystem(F,a,b,U0,steps).
Here [a,b] is the interval where the approximate solution is sought, and steps is the number
of subintervals. For the n-dimensional system U
0
(x) = F(x, U), U is an n × (steps + 1)
result matrix with the discrete approximations, and U0 is a column vector with the initial
condition. The RHS function F can be specified in a separate m-file (see driver1.m and
EqSystemTest.m).
(a) Implement the 2-stage Ralston method in a Matlab m-file that you name
RalstonSystem.m, with function heading
function U=RalstonSystem(F,a,b,U0,steps).
Test your implementation with driver1.m, and include the plot it generates
in your assignment package to be submitted on Crowdmark. What do you
observe, which of FE or Ralston gives the most accurate solution? Explain
why. Also, include a copy of your RalstonSystem.m in your assignment package.
(b) Compare the order of accuracy of the two methods applied to the test problem specified
in EqSystemTest.m. (Extend driver1.m to do this.) Define the error
ei = kUi − Vik2
where Ui
is the exact solution vector and Vi the numerical approximation at discrete location xi
. For each of the FE and Ralston methods and for steps = 16, 32, 64, 128, 256,
make a table of the error ei at the endpoint of the interval (i.e., i is the index
of the endpoint of the interval). What can you conclude about the global
order of accuracy for each of the methods? Include a copy of your extended
driver1.m in your assignment package.
(c) Consider the Lorenz equations, which are given by
x
0
(t) = σ (y(t) − x(t)),
y
0
(t) = r x(t) − y(t) − x(t)z(t),
z
0
(t) = x(t)y(t) − b z(t),
where σ, b and r are fixed parameters ∈ R (assumed to be positive) and t is time.
These equations are derived from fluid mechanics, and are related to weather prediction.
The equations are nonlinear (e.g., there are products x(t)z(t) and x(t)y(t)) and there
are no closed-form solutions, so numerical solutions are required.
2
– Implement the RHS function F of the Lorenz system in a Matlab m-file that you
name EqLorenz.m, with function heading
function F=EqLorenz(t,U),
where U(1) = x, U(2) = y, and U(3) = z and x, y, z are functions of t. Use
parameters σ = 10, b = 8/3, r = 28. Include a copy of your code in your assignment
package to be submitted on Crowdmark.
– Write a Matlab pogram driverLorenz.m that uses the Ralston method (or the
system Forward Euler method if you did not get Ralston to work) to calculate
and plot solutions for the Lorenz system. (You can adapt driver1.m.) For initial
condition (x0, y0, z0) = (1, 0, 0) and t ∈ [0, 50], produce a plot of x(t) and a plot of
the solution trajectory (x(t), z(t)) in the xz-plane. A good value for the time step
is ∆t = 0.002. (With larger time steps, the method may become unstable, or the
accuracy may be too low.) Include a copy of your code driverLorenz.m and the
two plots (generated by your code) in your assignment package.
– Hmm, the solution you generated looks very strange. Just like Lorenz in the 1960s,
you have re-discovered (mathematical) “Chaos”! (Lorenz also used the FE method
when he first observed Chaos numerically.) You can look up “Chaos theory” on
Wikipedia and have a look at some of the properties of this phenomenon as described there. Your plot (x(t), z(t)) in the xz-plane reflects a “strange attractor”.
The solution remains in a bounded region of space, but it does not converge to an
equilibrium point or cycle, and it appears to fill up the space locally. The attractor
has a fractal structure: it fills up the local 3D (x, y, z)−space more than a plane
(dimension 2), but not entirely (dimension 3); its Haussdorf dimension has been
measured to be about 2.06.
One of the properties of chaos is sensitivity to initial conditions: small differences in
initial conditions (such as those due to rounding errors in numerical computation)
yield widely diverging outcomes (for a solution that remains bounded), rendering
long-term prediction problematic. This has been summarized as: “Chaos: When the
present determines the future, but the approximate present does not approximately
determine the future.” Following Lorenz’ 1972 paper with title “Predictability: Does
the Flap of a Butterfly’s Wings in Brazil set off a Tornado in Texas?”, sensitivity
to initial conditions is now popularly known as the “butterfly effect”.
You are asked to illustrate the property of “sensitivity to initial conditions”: change the x0 initial condition to x0 = 1.0001, and make a plot that
overlays the values of x(t) with x(t) for the case x0 = 1. (Use “hold on” and diffferent colours or linestyles.) Submit this plot as part of your assignment package, with
a brief discussion (at most one paragraph) on how the plot illustrates sensitivity to
initial conditions.
2. Adaptive RK45 method. (20 marks)
Download the matlab m-files adaptiveRK.m, driver2.m and FBrusselator.m from the class
website. You are asked to implement the adaptive RK45 method that was described in class,
and apply it to a nonlinear ODE system called the ‘Brusselator’ ODE system
U
0
(t) = F(t, U(t)),
given by
3
u
0
1
(t) = 1 + u1(t)
2u2(t) − 4u1(t),
u
0
2
(t) = 3u1(t) − u1(t)
2u2(t).
This system describes autocatalytic reactions in chemistry. The Brusselator IVP you will
solve is given by
U
0
(t) = F(t, U(t)),
u1(0) = 1.5 and u2(0) = 3,
t ∈ [0, 20],
and the numerical solution you will reproduce looks like this:
(a) The adaptive RK45 method uses a Runge-Kutta base method, method A, with global
accuracy of order 4, and a secondary Runge-Kutta method, method B, with global
accuracy of order 5. For efficiency, methods A and B each use 6 stages where the
intermediate RK slopes Kj , j = 1, . . . , 6 computed in time step n are shared between
the two methods. Using the slopes Kj , j = 1, . . . , 6 computed in step n with step length
hn, methods A and B update the solution at time tn+1 from the solution at time tn using
a final slope that is a linear combination of the slopes Kj , j = 1, . . . , 6, as follows:
U
(A)
n+1 = Un + hn
X
6
j=1
w
(A)
j Kj
,
and
U
(B)
n+1 = Un + hn
X
6
j=1
w
(B)
j Kj
.
4
The coefficients w
(A)
j
for method A are chosen such that the global order for method
A is 4, and the coefficients w
(B)
j
for method B are chosen to produce an approximation
with global order 5. The slopes Kj , j = 1, . . . , 6 in step n are computed by
Kj = F
tn + hnαj , Un + hn
X
j−1
i=1
βijKi
!! ,
where the coefficients βij and αj are shared between methods A and B. The specific coefficients βij , αj , w
(A)
j
and w
(B)
j we will use can be found in file adaptiveRK.m. All details
about the input and output arguments required are also given in the file adaptiveRK.m,
and the script driver2.m provides the input parameters that will allow you to reproduce
the adaptive simulation of the so-called ‘Brusselator’ ODE system shown above.
You are asked to implement the adaptive RK45 code into file adaptiveRK.m, and submit
the file to the LEARN drop box on the course webpage. Also, include a copy of the code
in your assignment package to be submitted on Crowdmark.
(b) Run your adaptive RK45 code applied to the Brusselator IVP with initial timestep 0.1
and local truncation error tolerances δ = 10−2
, δ = 10−4 and δ = 10−6
.
– For all three values of δ, make plots of U as a function of t.
– For δ = 10−6
, make a phase-space plot, i.e., plot (u1(t), u2(t)) in the (u1, u2)-plane.
How does the solution behave for large t?
– For δ = 10−6
, how many function evaluations (computations of a Kj ) does the
algorithm perform? What is the smallest step taken by the adaptive method?
– One way to compare the efficiency of the adaptive method with a fixed timestep
method, is to count how many function evaluations would be needed in a fixed step
method with the the same steplength as the smallest steplength in the adaptive
simulation. (This would guarantee that local errors would not be much larger in
the fixed step method than in the adaptive simulation.) Make this comparison for
δ = 10−6
, and conclude whether an adaptive or a fixed step calculation is most
efficient for this IVP.
As always, include all answers, plots and code in your Crowdmark submission.
3. Rootfinding methods. (20 marks)
(a) Implement each of the following four root finding algorithms discussed in class in MATLAB in a seperate m-file. You should use the function header given for each below.
(a1) Bisection method.
function root = bisection(f, a, b, maxit)
% Usage: bisection(f, a, b, maxit)
% This method uses bisection to find roots of
% f in the interval [a,b] with maximum number of steps maxit.
(a2) Fixed point iteration.
function root = fixedpoint(g, x0, maxit)
% Usage: fixedpoint(g, x0, maxit)
% This method uses fixed point iteration to
% find fixed points of g with maximum number of steps maxit
% using initial estimate x0.
5
(a3) Newton’s method.
function root = newton(f, fprime, x0, maxit)
% Usage: newton(f, fprime, x0, maxit)
% This method uses Newton’s method to
% find roots of f (with derivative function fprime)
% with maximum number of steps maxit
% using initial estimate x0.
(a4) Secant method.
function root = secant(f, x0, x1, maxit)
% Usage: secant(f, x0, x1, maxit)
% This method uses the secant method to
% find roots of f with maximum number of
% steps maxit using initial estimates x0 and x1.
Write an m-file driver script, question3.m, in which you use the four functions to find
the root of
f(x) = x
2 − exp(−x)/2.
The script should produce a table of the approximate solution for successive steps of
each of the four methods, and a table of the absolute errors for all steps. Use format
long. Use enough steps to reach full accuracy for all methods. You can use the inline
function to construct f, f prime and g in the driver script (or you can use separate mfiles). For the fixed point method, use g(x) = x − f(x). For bisection, use [a, b] = [0, 2]
as the initial interval. Use initial guess x0 = 1 for the fixed point method, and x0 = 2
for Newton. For the Secant method, use x0 = 2 and x1 = 0.
Investigate the convergence order and convergence constants for the Newton, Secant and
fixed point methods. Make tables of the error ratios ck for each step for each of these
methods, and verify that ck asymptotically tends to a constant. Verify that this is the
theoretically predicted constant for the Newton and fixed point methods.
Submit your driver and function m-files to the LEARN drop box on the course webpage.
Also, include a copy of the code in your assignment package submitted on Crowdmark.
Your driver should produce the tables requested. Also provide printouts of your tables
in your Crowdmark submission.
You will next apply your methods to a few more examples.
(b) Consider f(x) = (x − π)|x − π|. This function has one root, located at x = π. Plot
this function using the ezplot command on the interval [2, 4]. Ensure that your plot
has appropriate title and axis labels. Using the initial estimate x0 = 2 apply Newton’s
method and the Secant method (with x1 = 4) to find the root at x = π. What do you
notice about the converge orders of each method? Explain why this behaviour does not
contradict what you know about convergence of these methods.
(c) Consider
f(x) = −
exp(3x)
9
−
x
2
2
+
1
9
.
Plot f(x) and g(x) = x − f(x) using ezplot and provide the resulting plot. Using fixed
point iteration on g(x), create a table of the results you get with initial condition x0 = 1.
Why are we unable to achieve convergence? Suggest a different g(x) and x0 that gives
convergence with fixed point iteration.
6
(d) Consider f(x) = −|x|
1/3
. Plot the function using ezplot and provide the resulting plot.
Using Newton’s method, create a table of the results you get with initial condition x0 = 1.
Explain graphically why we are unable to achieve convergence. Show analytically that
Newton’s method diverges for any initial guess. Does this contradict what you know
about convergence of Newton’s method?
(e) Consider f(x) = x/(x
2 + 1). Plot the function using ezplot and provide the resulting
plot. Using Newton’s method, create a table of the results you get with the initial
condition x0 = 0.75. What do you observe? Does this contradict what you know about
convergence of Newton’s method? How can you fix this?