CSC336-Assignment 3 Solved

CSC336S                                                                                  Assignment 3                     

Please write your family and given names and underline your family name on the front page of your paper.

All questions must be written in latex, and compiled to a single pdf. Plots (7 for Q1 and 2 for Q3) must in embedded in latex. All code and output of Q1, Q3 and Q4 should be embedded into latex. Code and output should be embedded with fixed-width fonts, e.g. Courier. Font size of all fonts must be 12. Linespacing set to 1.1.


Thus, the code and plots will be available within latex/pdf, as well as separately.

1. This question involves the study of the outbreak of an epidemic in a population and its eradication with vaccinations. The study of spreading of infectious diseases is mainly done by numerical Ordinary Differential Equations (ODEs), which are not the topic of CSC336 (but of CSC436). However, numerical ODE solvers, essentially convert a system of ODEs to a sequence of nonlinear systems. Each nonlinear system of the sequence relates unknown quantities at a time point (or step), in terms of known quantities at the previous time point, while the initial state is given (known). In this way, the values of various quantities at several time points are computed and studied.

Assume there is a population of N people, which remains constant over time. Assume that each person belongs to one of the five groups: susceptible, exposed (to become infected after incubation), infected/ious, recovered/immune and vaccinated/immune. The infected/ious persons transmit the virus at a given rate β to the susceptible and make them exposed, the exposed become infected/ious after a certain incubation period α−1 (also given), the infected/ious recover at a given rate γ, and stay immune thereafter. Furthermore, the susceptible get vaccinated at a given rate ρ, and become immune immediately after vaccination and stay immune thereafter. The population of N persons is replenished at a given death and birth rate µ, and all newborns are susceptible.

In an instance of this problem, at each time point i, the following 5 × 5 nonlinear system needs to be solved:

                    x1(i) = x1(i−1) + hµN − hµx1(i) − h β x1(i) x3(i) − hρx1(i)                                                                                                                                                             (1a)

N

                          (i)           (i−1)                   (i)                  (i)                 β (i)       (i)

                                x2 = x2 − hα x2 − hµx2 + hx1 x3                                                                                                                                                                                                    (1b)

N

x3(i) = x3(i−1) − hγ x3(i) − hµx3(i) + hα x2(i)   (1c) x4(i) = x4(i−1) + hγ x3(i) − hµx4(i) (1d) x5(i) = x5(i−1) + hρx1(i) − hµx5(i). (1e)

  In the above, x(i) = [x1(i), x2(i), x3(i), x4(i), x5(i)]T represents the vector of numbers of susceptible, exposed, infected/ious, recovered/resistant/immune and vaccinated/immune persons, respectively, at time point i, and these are the five unknowns that need to be computed at the ith time point (step), assuming the quantities of the previous time point i − 1 are already computed, and the initial state, x(0) = [x1(0), x2(0), x3(0), x4(0), x5(0)]T is given. While in practice x(ji) are integers, we treat them as continuous variables. Also, h is a small stepsize in time, that, in this instance, is given and fixed. Also note that, although some simplifications can be applied to this system, we just treat it as a general 5 × 5 nonlinear system. Finally, note that it is not required that you understand the details of how this nonlinear system arises. You take the system as granted.

 

(a)Write the system (1a)-(1e) in standard form f (x(i)) = 0. Indicate all components of f . Write (by hand, i.e.

in latex) the Jacobian matrix for this system. Indicate all entries (some in terms of x(ji))).

(b)Consider Newton’s method for this system. One iteration of Newton’s includes computation of the form solve J(x(i,k−1))s(i,k−1) =− f (x(i,k−1)) compute x(i,k) = x(i,k−1) + s(i,k−1)

Note that i is the time point index, while k is the index for Newton’s iteration, the former remaining constant through all Newton iterations (all k).

  Write a matlab or equivalent script that, given N = 38e6, α−1 = 6, β = 0. 25, γ = 0. 06, µ = 0. 01/365, ρ = 300µ, h = 1/2, x1(0) = N − x2(0) − x3(0) − x4(0), x2(0) = 20e3, x3(0) = 30e3, x4(0) = 850e3, and x5(0) = 0, computes x(1) by Newton’s method with tolerance 10−6. Use maxit = 10, and as stopping criterion the infinity norm of the residual vector. The Jacobian must be solved using backslash. At each Newton iteration, output x(1,k) (the five values computed), the sum of the five values and the infinity norm of the residual. See script testnl.m (in the course website) for some template. Comment on the number of Newton iterations and on how the residual behaves as the iterations proceed.

- 2 -

(c)Consider now a simulation time period from 0 to tend = 150 (can be viewed as a period of 150 days), divided

tend in nstep =         periods of stepsize (length) h = 1/24 (or dt), and that, at each time point i, giv en the state at time point h

    i − 1, a nonlinear system of the form (1a)-(1e) is solved by Newton’s method with tolerance 10−6. Write a script that, given the same parameters as in (b), except that ρ = 0 (and h = 1/24), computes x(i), for i = 1,..., nstep. Do NOT output and do NOT sav e the results of each Newton iteration. Save, but do NOT output, the number of Newton iterations at each timestep. Save, but do NOT output, the results (x(i)) from each time point. They will be used later for plotting. Output the initial values x(0) and their sum, the computed final x(nstep) (five values at time tend), and their sum, as well as the maximum number of each type of persons among all time points (steps). Output also which day the maximum number of infected persons occurs.

Plot x1(i), x2(i), x3(i), x4(i) and x5(i) versus ih, i = 0,..., nstep (index of time point scaled by h), in one plot, using solid (red), dashed (green), dashed-dotted (blue), dotted (black), and plain dot (cyan, ’c.’) lines. Add proper legend, labels for axes, etc. Use axis tight; after the plot.

In another plot, plot the number of Newton’s iteration at each timestep, versus the index of the timestep. Add labels for axes. Use axis([1 nstep 0.9 2.1]); after the plot. See script nlflu.m (in the course website) for some template.

(d)Do another simulation for time period from 0 to tend = 365, with h = 1/24, and otherwise the same parameters as in (b), i.e. with ρ = 300µ. Output the same quantities and do the same plots as in (c), with the same axis specifications.

(e)Do another simulation for time period from 0 to tend = 365, with h = 1/24, and with the same parameters as in (b) (i.e. with ρ = 300µ), except that β = 0. 25/4, as if measures are taken to reduce contacts by a factor of 4. Output the same quantities and do the same plots as in (c), with the same axis specifications.

In addition, in a another plot, plot x2(i), x3(i), versus ih, i = 0,..., nstep, in one plot, using dashed (green) and dashed-dotted (blue) lines. Add proper legend, labels for axes, etc. Use axis tight; after the plot.

Place the first plot from (c) and the first from (d) side-by-side. Place the first and last plots from (e) side-by-side. Place the three plots of number of iterations side-by-side in 1 × 3 format. Comment on all results of (c), (d) and (e), including output and plots.

Notes: Do not use any symbolic calculation of any sort. The derivatives should be derived by hand and hard-coded into the Jacobian matrix.

2.         [To be done by hand/latex, no coding.] Let f (x) = x sin(x) − 1.

(a)Determine graphically the number of roots of f (x) = 0 in the interval [0, π]. (Hint: Rewrite equation f (x) = 0 to make graphing by hand easy.)

(b)Show mathematically that f has a unique root in [  ,  ].

(c)Consider the fixed-point iteration scheme x(k+1) = g(x(k)), with g(x) =  . Show mathematically that g has a unique fixed point, say x *, in [  ,  ].

(d)Find a (closed) interval I of length at least  , that contains x *, so that, if x(0) ∈I, the iteration x(k+1) = g(x(k)) converges to x *. Specify I and explain. (You can specify the endpoints of I in terms of π.) (e) [4 points] What is the order of convergence of the iteration scheme in (d)? Justify mathematically.

(f)Show that the iteration x(k+1) = g(x(k)) converges to x *, if x 

(Note: If you have already shown this in (d), you can skip this question.)

 (0) = π , and indicate the computed (g) [3 points] Apply by hand one Newton iteration to f (x) = 0, with starting guess x

2 x(1).

3.        [To be done by hand/latex, except (e).]

(a)Use Newton’s Divided Differences, to construct a interpolating polynomial p1(x) of the function f (x) = e−x (i.e. f (x) = exp(−x)), based on the data points x0 =− 1 and x1 = 1. Show the NDD table, and the interpolating polynomial p1(x).

Note: It is fine to leave quantities such as e−1 and e in the presentation of the table entries and polynomial coefficients, as long as you reasonably simplify all terms.

- 3 -

(b)Giv e the polynomial interpolation error formula for this problem. Note that the formula involves an unknown point ξ.

Using the formula, find a (sharp) bound of the (absolute value of the) error for any x∈[−1, 1]. Indicate the interval ξ belongs to, when x∈[−1, 1].

Using the formula, find a (sharp) bound of the (absolute value of the) error for any x∈[1, 2]. Indicate the interval ξ belongs to, when x∈[1, 2].

(c)Use Newton’s Divided Differences, to construct the lowest degree interpolating polynomial p2(x) of the function f (x) = e−x (i.e. f (x) = exp(−x)), based on the data points x0 =− 1, x1 = 0 and x2 = 1. Show the NDD table, and the interpolating polynomial p2(x).

Note: It is fine to leave quantities such as e−1 and e in the presentation of the table entries and polynomial coefficients, as long as you reasonably simplify all terms.

(d)Giv e the polynomial interpolation error formula for this problem. Note that the formula involves an unknown point ξ.

Using the formula, find a (sharp) bound of the (absolute value of the) error for any x∈[−1, 1]. Indicate the interval ξ belongs to, when x∈[−1, 1].

Using the formula, find a (sharp) bound of the (absolute value of the) error for any x∈[1, 2]. Indicate the interval ξ belongs to, when x∈[1, 2].

(e)Giv e the Taylor polynomials t1(x) and t2(x) of degree 1 and 2, respectively, that arise by expansion of f (x) = e−x about x = 0. Write a MATLAB script that uses polyfit to compute the polynomials p1(x) and p2(x), interpolating f (x) at 2 and 3 equidistant data points in [−1, 1], respectively. The script then evaluates the four polynomials, p1, p2, t1, t2, and the function f at 100 equidistant evaluation points in [−1, 1]. Let v1, v2, u1, u2 and v f be the vectors of values of p1, p2, t1, t2 and f , respectively, at 100 equidistant evaluation points in [−1, 1]. In one plot, graph v1, v2, t1 and t2 and v f versus the evaluation points. Use solid, dashed, dotted, dotted-dashed and solid lines, respectively, and appropriate legend and labels. In another plot, graph v f − v1, v f − v2, v f − u1, v f − u2, versus the evaluation points. Use solid, dashed, dotted and dotted-dashed lines, respectively, and appropriate legend and labels. Calculate and output max |v f − v1|, max |v f − v2|, max |v f − u1|, and max |v f − u2|,. Place the plots side-by-side with appro-

                                eval pts                           eval pts                           eval pts                                      eval pts

priate captions/subcaptions. Comment on the results.

4.Consider using three types of interpolation for the function f (x) = exp(−x), namely, polynomial, linear spline and cubic spline interpolation, as in script 

Add the appropriate extra lines to the script to compute the cubic spline interpolant at the points indicated. Then run the script and collect the output. (If you get a warning from polyfit, you are not necessarily doing anything wrong.) Comment on the results, in particular, how the error of each interpolant behaves as the number of data points increases, and whether this behaviour agrees with what expected from theory.

Do not change the format according to which the errors and other quantities are output. You need to explain what the other quantities represent.

General note: Plotting quantity y versus quantity x, means that x is in the x-axis and y is on the y-axis, i.e. what follows "versus" is in the horizontal axis.