CSCI 2227 - HW 7

Date: Fall 2023

Due Date: Tuesday Dec 5, 11:59pm EST

Update 11/30/2023: On the bonus problem, please demo your  how_much_vaccination_to_save()   function by printing its result for  x=.5, .7, .9 . In the first uploaded version I had  x = .5, .4, .3 instead (which doesn’t really make sense since the answer for   .4 and   .3 is 0).

Update 11/29/2023: The first uploaded version had a typo in the Written Exercise 1: a  was missing in the differential equation’s denominator, which threw off the true value.

How to submit

Submit your work on Gradescope

(https://bostoncollege.instructure.com/courses/1647277/external_tools/127906)as an IPython notebook (  .ipynb file). I suggest you just use Google Colab

(https://colab.research.google.com/)(see Homework 1 assignment for more details).

Honor code

At the top of your notebook, please include a text cell and paste in the following honor code statement:

I expect this to be a true statement for everyone. If any part of this is not true, then please explain the situation in this text cell.

Written exercises

Formatting for Written Exercises

The written exercises typically involve some writing, and some coding, so use a combination of    Text cells and Code cells. Please preface each of your solutions with a Text cell containing e.g. # Written Exercise 1  . The leading  #<space> is important, since it creates a bold heading.

To type math, use Latex commands in Text cells between dollar signs, such as $1+x+x^2 = \sum_{n=0}^{2} x^n$ . Here is a good cheat sheet

(https://joshua.smcvt.edu/undergradmath/undergradmath.pdf), here is a good Latex Equation   Editor     (                                                                                                        )https://www.tutorialspoint.com/latex_equation_editor.htm . Ask on Piazza if you need more help with Latex.

Written Exercise 1 “ Leaps and Bounds” (5 points)

Here is an IVP we want to solve:      and            , and we want to estimate  .

Part (a) (1 point)

Solve this IVP using one step of Euler’s Method. You can do the math by hand or write a short    Python script. Don’t bother with an all-purpose  euler() function yet: you’re only doing one step.

Part (b) (2 points)

Solve this IVP using one step of the RK2 Method (with ). You may do the math by hand or write a short script similar to part (a).

Part (c) (2 points)

We’re now going to analyze the error term when doing  steps of Euler’s Method. The true value of  is, in fact,  , because the true solution turns out to be  .

Write a Python script to solve the above IVP using  steps of Euler’s Method, where

N=1,2,4,8,…,406 (you can and should use an all-purpose  euler() function this time). Make an aligned table with three columns: , your estimate, and the absolute error.

How fast are the errors going to zero? Answer by computing the ratio of successive errors, and

explain (in a print statement or Text cell) what this means for the order of the error term in terms of N. You don’t have to create a table showing the ratio of successive errors, but do mention them in your explanation.

Written Exercise 2 “Cooling Coffee” (5 points)

In Chapter 1, we introduced Newton’s Law of Cooling:

T(0)=TO

This differential equation models the temperature  of an object immersed in an environment of constant temperature  .

Part (a) (2 points)

We will use the constants  ,  (you can use 0.693), and  (temperatures in degrees Fahrenheit). These values are chosen to represent a cooling cup of coffee/tea in a

room-temperature office (the choice of  was chosen to make the math easier). In this problem, you will estimate the temperature after 5 minutes.

Use the composite RK2 method to estimate  using  steps. When you write up RK2, use     a＝0-5                                                                             [0,5.Makeagraphofyourestimatedsolutioncurveovertheinterval . Make sure your graph has an x-axis going from  to  . Add labels to the x-axis and y-axis.

Feel free to beautify your plot further (not required).

Part (b) (1 point)

How quickly does theory predict the error in the composite RK2 method will go to  in terms of the number of steps ? Give a short theoretical argument, similar to the one we gave for Euler’s

Method in class.

Part (c) (2 points)

Verify your answer in part (b) holds in this Newton’s Law of Cooling example. To do so, use the

composite RK2 method to estimate  using  steps, where  .    Make a table of these estimates, and add a column showing the absolute error (for this you’ll have to figure out the true solution to Newton’s Law of Cooling, which is on pg. 145 of the textbook or     you can derive it yourself with separation of variables).

How fast are the errors going to zero? Answer by computing the ratio of successive errors, and

explain (in a print statement or Text cell) what this means for the order of the error term in terms of N. You don’t have to create a table showing the ratio of successive errors, but do mention them in your explanation.

Written Exercise 3“Euler Slick”(5 points)

Consider the following IVP:      and      and we want to know  . This looks like a perfectly harmless IVP, so let’s solve it!

Part (a) (1 point)

Do 1 step of Euler’s Method by hand to estimate  .

Part (b) (1 point)

The general solution to      is in fact    . What is the correct value of  based on the initial condition     ? What do you conclude is the error in your estimate from part (a)?

Part (c) (1 point)

Can this error be reduced by doing more than 1 step of Euler’s Method? To answer, try doing 999 steps of Euler’s Method over the interval  and graph the resulting solution curve with

Python. Are we getting closer to the true solution at ?

Part (d) (2 points)

Try to explain what went wrong here. After all, the IVP looked perfectly harmless!

Hints:

·  What happens to the true solution curve if you change the initial condition by 0.000001?

  Recall that at one point I said in class: “After you do one step of Euler’s Method, you’re not on  the true solution anymore.” Unpack what this means, considering the previous bullet point, and you’ll find out what is going wrong here.

Coding Exercises

Coding Formatting

Type your code into cells in your notebook. Preface your solution with a Text cell containing  # Coding Exercise 1 . Unless otherwise indicated, you don’t need to explain your answers.

Coding Exercise 1“Hundred Steps”(5 points)

Recall the Lotka-Volterra Predator Prey model we simulated in class:

R(t)＝   0.7R(t)-0.005R(t)F(t)

F(t)＝-0.2F(t)+0.001R(t)F(t)

The initial condition is  and  .

Make a graph of just the fox population over the interval   [0,59] using 100 steps of vectorized

Euler’s Method. You can use the vectorized Euler’s Method from the class slides, or you can write your own.

Add to this graph the solution produced by a vectorized RK2 method, 100 steps. You’ll need to write your own vectorized RK2 method (use ).

Also add to this graph the solution produced by the vectorized RK4 method, 100 steps. You can use the Vectorized RK4 Method from the textbook, or you can write your own.

Make sure your graph has the following features:

·  Different colors or patterns for the three different curves.

·  A legend indicating which curve is which.

  Informative x-axis and y-axis labels.

Coding Exercise 2“Controlling the Curve”(8 points)

Part (a) (2 points)

In this problem, you will simulate the following 4x4 differential equation over the interval te [0,15o (use 9999 steps of Euler’s Method throughout).

This is the usual SIR model with one additional group added: the “Vaccinated” population (they get to go straight from susceptible to immunity, without having to get infected first).

For initial condition, please start with  1/8000 in the infected group,  1-1/8000 in susceptible, and the other groups  0 .

We’re going to be making a lot of plots where we tweak the parameters ,  , and  , so write a function:

which plots the SIRV model using the current values of the variables  b, k, p . You should be able to reproduce the plot we made in class with:

We’restarting with , which is the regular SIR model (no vaccines yet).

Part (b) (2 points)

During the COVID-19 pandemic, health officials aimed to “flatten the curve”: the idea was that via  social distancing or other measures, we could lower the strain on health care systems by reducing the maximum amount of people that were ever simultaneously infected.

Suppose that, with no social distancing,  b=1/2 and  k=1/3 are realistic parameter choices. If     everyone starts social distancing, which of these parameters would this likely change? Would it increase or decrease?

Try tweaking the relevant parameter by  0.1 and make a new plot of the results. Comparing the new plot to the part (a) plot, did it have the desired effect of “reducing the maximum amount of   people that are ever simultaneously infected”?

To answer this quantitatively, go back and add some print statements to  simulate_SIRV() which output:

·  The maximum amount (proportion) that were ever infected simultaneously (this is the maximum height of “the curve” we want to flatten).

·  The proportion that were never infected (including the vaccinated), by the end of our simulation.

These stats can be computed in one line each.

Part (c) (2 point)

Let us call “Parameter 1” your answer to part (b) and “Parameter 2” the other

parameter(whichever is the one among  b and  k you didn’t answer in part (b)). From a public

health standpoint, do we want Parameter 2 to increase or decrease? Give one plausible example  of an intervention/innovation that you think could produce this desired effect. Put your answer in a comment or a text cell.

Put all the parameters back at their part (a) values, then try tweaking this Parameter 2 by  0.05  in the direction you determined is favorable for public health, and make a plot.

Which of the printed stats from part (b) are impacted (compared to part (a))? Give a brief

explanation (in comment or Text cell) of why changing this parameter would have this impact.

Part (d) (2 points)

Vaccination is controlled by parameter  p . Leaving all the other constants as they were in part (a), what value of this parameter is needed so that about 50% of the population never suffered any

infection at all? You can estimate this by guess and check, try to get it within 3 decimal places of 50%.

(You will notice that if even just a few people get vaccinated, it generally has a massive positive impact on everyone, including the non-vaccinated).

Part (e) (bonus problem, worth 1 point of extra credit)

If you’re up for it, write a function called  how_much_vaccination_to_save() which takes a single

parameter  x (between  and ) and returnshow large the vaccination parameter has to be so  that proportion  x of the population never suffers any infection at all (vaccinated individuals also count as never suffering any infection). Your return value will only be an estimate, not the exact  value.

Notice that this can be done using the nonlinear solvers from Chapter 5 (did your guess-and-   check in part (d) remind you of anything?). But it’s also not that easy to set up the function  f() that you want to be  0 , so that’s why it’s a bonus problem.

Demo your  how_much_vaccination_to_save() function by printing its result for  x=.5, .7, .9 .

Each exercise will receive a score out of 5 points (unless otherwise indicated).

·  5 means perfect/exemplary understanding.

  4 means mainly good understanding, maybe a small mistake or two.

  2-3 means an understanding that is in the right direction but not correct.

·  0-1 means mainly incorrect, not close to correct.

Remember, there is a -25% late penalty per day of lateness.