[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.

What to submit:

• The single pdf file 00A3.pdf (with embedded code and output).
• The latex source file of A3.
• Any source code you wrote (for computation, plots, etc).
• Plot files (7 for Q1 and 2 for Q3) in eps and png format.

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 ¹ (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:

^x₁(i) = ^x₁(i1) + hN h^x₁(i) h ^x₁(i) ^x₃(i) h^x₁(i) (1a)

N

(i) (i1) (i) (i) (i) (i)

x₂ = x₂ h x₂ hx₂ + hx₁ x₃ (1b)

N

^x₃(i) = ^x₃(i1) h ^x₃(i) h^x₃(i) + h ^x₂(i) (1c) ^x₄(i) = ^x₄(i1) + h ^x₃(i) h^x₄(i) (1d) ^x₅(i) = ^x₅(i1) + h^x₁(i) h^x₅(i). (1e)

In the above, x⁽ⁱ⁾ = [x₁⁽ⁱ⁾, x₂⁽ⁱ⁾, x₃⁽ⁱ⁾, x₄⁽ⁱ⁾, x₅⁽ⁱ⁾]^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⁽⁰⁾ = [x₁⁽⁰⁾, x₂⁽⁰⁾, x₃⁽⁰⁾, x₄⁽⁰⁾, x₅⁽⁰⁾]^T is given. While in practice x⁽_jⁱ⁾ 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.

• [5 points] Write the system (1a)-(1e) in standard form f (x⁽ⁱ⁾) = 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⁽_jⁱ⁾⁾).

• [8 points] Consider Newtons method for this system. One iteration of Newtons includes computation of the form solve J(x(i,k1))s(i,k1) = f (x(i,k1)) compute _x(i,k) = _x(i,k1) + _s(i,k1)

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

Write a matlab or equivalent script that, given N = 38e6, ¹ = 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 Newtons method with tolerance 10⁶. 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

• [7 points] Consider now a simulation time period from 0 to t_end = 150 (can be viewed as a period of 150 days), divided

t^end 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 Newtons method with tolerance 10⁶. Write a script that, given the same parameters as in (b), except that = 0 (and h = 1/24), computes x⁽ⁱ⁾, 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⁽ⁱ⁾) from each time point. They will be used later for plotting. Output the initial values x⁽⁰⁾ and their sum, the computed final x^(nstep) (five values at time t_end), 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 x₁⁽ⁱ⁾, x₂⁽ⁱ⁾, x₃⁽ⁱ⁾, x₄⁽ⁱ⁾ and x₅⁽ⁱ⁾ 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 Newtons 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.

• [5 points] Do another simulation for time period from 0 to t_end = 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.
• [15 points] Do another simulation for time period from 0 to t_end = 365, with h = 1/24, and with the same parameters as in (b) (i.e. with = 300), except that = 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 x₂⁽ⁱ⁾, x₃⁽ⁱ⁾, 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.

• [1 points] 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.)
• [2 points] Show mathematically that f has a unique root in [ , ].
• [1 points] 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 [ , ].
• [8 points] Find a (closed) interval I of length at least , that contains x *, so that, if x⁽⁰⁾ 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) [6 points] 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.)

⁽⁰⁾ = , 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).]

• [3 points] Use Newtons Divided Differences, to construct a interpolating polynomial p₁(x) of the function f (x) = e^x (i.e. f (x) = exp(x)), based on the data points x₀ = 1 and x₁ = Show the NDD table, and the interpolating polynomial p₁(x).

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

3

• [3 points] 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].

• [4 points] Use Newtons Divided Differences, to construct the lowest degree interpolating polynomial p₂(x) of the function f (x) = e^x (i.e. f (x) = exp(x)), based on the data points x₀ = 1, x₁ = 0 and x₂ = Show the NDD table, and the interpolating polynomial p₂(x).

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

• [5 points] 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].

• [10 points] Giv e the Taylor polynomials t₁(x) and t₂(x) of degree 1 and 2, respectively, that arise by expansion of f (x) = e^x about x = Write a MATLAB script that uses polyfit to compute the polynomials p₁(x) and p₂(x), interpolating f (x) at 2 and 3 equidistant data points in [1, 1], respectively. The script then evaluates the four polynomials, p₁, p₂, t₁, t₂, and the function f at 100 equidistant evaluation points in [1, 1]. Let v₁, v₂, u₁, u₂ and v _f be the vectors of values of p₁, p₂, t₁, t₂ and f , respectively, at 100 equidistant evaluation points in [1, 1]. In one plot, graph v₁, v₂, t₁ and t₂ 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 v₁, v _f v₂, v _f u₁, v _f u₂, versus the evaluation points. Use solid, dashed, dotted and dotted-dashed lines, respectively, and appropriate legend and labels. Calculate and output max |v _f v₁|, max |v _f v₂|, max |v _f u₁|, and max |v _f u₂|,. Place the plots side-by-side with appro-

eval pts eval pts eval pts eval pts

priate captions/subcaptions. Comment on the results.

4. [10 points] Consider using three types of interpolation for the function f (x) = exp(x), namely, polynomial, linear spline and cubic spline interpolation, as in script http://www.cs.toronto.edu/ccc/Courses/336/scriptint.m

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.