These are draft notes from subtitles. Please help to improve them. Thank you!
How can we model the spreading of a contagious disease? The straight forward approach may be to model individual persons as either infectious or susceptible If one infectious and one susceptible person meet, there is a certain chance that the susceptible person becomes infectious too and then the first infectious person as well as the newly infected person further spread the disease. This seems to cry out loud for a spatial model in which we take care of the location of these persons that would be complex. Luckily, there is another way. This is a much simpler way of thinking about such diseases. The SIR Model. We do not look at the spatial distribution at all. We simply say that our population is made up of susceptible, infectious, and recovered persons and we only look at the total number of these persons. Susceptible persons can get infected and then become infectious. The rate of that process is controlled by the number of infectious people. The more infectious people there are, the larger the probability is that a susceptible person meets an infectious person. Hopefully, people are going to recover from the disease, so that is the process that takes us from infectious to recovered. And then the standard version of that model, we assume that the recovered persons are immune and cannot be infected again. So they form a third compartment. This is a much simpler way of dealing with our problem. We forget about the spatial distribution altogether and simply construct a mechanical model if you will in which the population consists of three compartments, susceptible persons, infectious persons, recovered persons depending on the number of infectious persons, susceptible persons become infected and hence infectious and after a while, infectious persons recover and we assume that recovered persons are immune against the disease and cannot be infected again. This is the SIR Model. Susceptible, Infectious, Recovered.
Now, let's turn the SIR model into equations. We're doing that one part after the other--first, the recovery. If I(t) is the number of infectious persons--it's I(t), and R(t) is the number of recovered persons--it's R(t), so which of these equations makes the most sense. The derivative of the number of recovered persons with respect to time equals 5 days times the number of infected persons 1/5 days times the number of infected persons 1/5 days times the number of recovered persons or minus 1/5 days times the number of infected persons.
With this answer, the rate of change of the number of recovered persons would be negative, meaning that that number shrinks. This cannot be true. That number has to grow, so we can rule out this answer. This derivative of a number of persons with respect to time has to be something about persons per days. It cannot be days times persons, so this cannot be true. The remaining two answers are correct in that respect. This answer which need to unbound its exponential growth. This works like compound interest. The rate of change is proportional to the current value. This cannot be true. So we are left with this answer, which make sense. For instance, if we have twice the number of infectious persons, this rate would go up by a factor of 2--twice as many persons will recover per day. And imagine you have 5,000 infectious persons, then you get a recovery rate of 5,000 persons/5 days, 1,000 persons per day that recover. These 5 days seem to be something like a typical duration for the disease. So this is right answer.
Now, let's look at the infection part of the SIR model. We've got a number of susceptible persons S and we've got a number of infectious persons I, and we watch the model the spread of the disease, so which of these equations makes more sense for the infection process. The derivative of the number of infectious persons with respect to time equals 5 * 10⁻⁹ per day times the number of susceptible persons minus 5 * 10⁻⁹ per day times the number of infectious persons 5 * 10⁻⁹ per day the number of susceptible persons divided by the number of infectious persons or 5 * 10⁻⁹ per day and person times the number of infectious persons times the number of susceptible persons.
This answer would simply lead to a decay of the number of infectious persons no matter what the rate of change is a negative factor times the current value. This doesn't make sense With this answer, the rate of change which depend only on the number of susceptible persons. It would not matter to the rate of change of the number of infectious persons whether we started with 1 or 10 or 1 million infectious persons. This doesn't make sense. So we have to have a solution that incorporates both quantities. The number of susceptible persons and the number of infectious persons. So what about that ratio? This would mean that if the number of infectious person is low, this ratio is high. We divide by a small number and the rate would be high which again doesn't make sense. And there is another reason why this doesn't makes sense, the unit would be 1/days times persons/persons. This canceled and unit of the result is 1/days but what we need is persons per day, so this cannot be true for several reasons. This is what is left and it makes sense. Let's look at the units. A number divided by days and persons, times persons squared. What we are left with is person per day, that's what we need and this product has the right behavior. If we increase the number of infectious persons, we increase that rate. We have more encounters. That makes sense.
We already look into the rate of change of the number of infectious persons. It was a constant times the number of infectious persons times the number of susceptible persons. If this is the number of persons per day that move from susceptible to infectious, the same number has to get lost from the susceptible compartment. So the rate of change of the number of susceptible persons has to be minus the same quantity, and then there's the second effect mainly that people recover. We've already looked at that term. This is the increase of the number of recovered persons per time. If this number of persons is added to the recovered compartment each day, we have to subtract it from the infectious compartment as well. We have to subtract it from the infectious compartment because these persons leave the infectious compartment and join the recovered compartment. This completes the equations of the SIR model. Obviously, this is a very rough model. In particular, as the numbers that we are dealing with here are not integer numbers, in real life they should be say 1,000 or 1,000 at 1 or 1,000 at 2 persons, but with these equations, we're dealing with say 1,000.1, 2 , 3, 4, 5 persons. These equations can only make sense for large numbers of persons so that the fractional part can be ignored, but if these numbers are large, the numerical constants that we're using here may for instance be used to model influenza.
Now, that we have the equations for the SIR model, compute the sum of rate of changes. The rate of change of the number of susceptible persons plus the rate of change of the number of infectious persons plus the rate of change of the number of recovered persons, then look at the result. What do you learn from that result. Do you learn from the result that an epidemic will develop or that no persons are lost or that everybody recovers eventually or that the number of infections per time stays constant.
The correct answer is that no persons are lost. Let's look into the reason for that. These are the creations of the SIR Model again. If we form the sum of these three derivatives and see that this term minus coefficient times I times S plus the same coefficient times I times S, cancelled in the sum and the same happens with this term, minus coefficient times I plus the same coefficient times I. So the sum of these derivatives is zero. For simple interpretation of that left-hand side here. If you look at this number, the number of susceptible persons plus the number of infectious persons plus the number of recovered persons form the derivative with respect to time, we get this expression. The rate of change of the sum is the sum of the rates of change. You see, in this case, I'm switching to the d/t notation because it's hard to place the dots on top of the sum. So what we learned is, that the rate of change of the sum is zero. If the rate of change is zero, there is no change and if there is no change, the sum is constant which means nobody is lost. Persons move from susceptible to infectious and from infectious to recovered but never move out of that system. We could see the same equation also if there was immigration, or births and deaths that compensate each other precisely.
If you look at the equations of the SIR model, the meaning of this constant in front of the I here is pretty clear. This number of 5 days controls the infectious period of the disease, but what about the first constant here that controls the rate of infection. Assume that we're looking at a secluded small town with 1,000 inhabitants in total and assume that everyone meets 25 persons per day, and third, assume that the infection probability is 2% per contact, which number should be entered here. Just input that number.
This number would be 0.0005 per day and person. Let's do the computation. Let's look at the given data. What do we know about each infected person? As every other person, each infected person that meets 25 persons a day. All this susceptible persons can be infected, so the question is how many susceptible persons are among these 25? We can compute that by ratio. If there are 1,000 people in total, then the ratio of the number of susceptible persons divided by the total of 1,000 persons times these 25 persons that would be susceptible. Imagine, everybody was susceptible then we would compute 1,000 persons times 25. All 25 would be susceptible. If every second person was susceptible, 500 divided by 1,000 results in one half times 25, every second person would be susceptible. So this is the number of contacts between a single infected person and susceptible persons per day and each of these contacts will lead to infection with a probability of 2%. So now, we can compute how many persons are infected by the single infected person per day versus the number of contact of susceptible persons per day and multiply with probability. This will be the number of new infections per day. We can cancel the persons, these three numbers combine to 0.0005 per day and per infected person, and that's precisely the expression we are needing in the SIR equations. Constant times the number of susceptible persons times the number of infected persons. But keep in mind that we are dealing with average values here, 25 persons per day give or take a probability of 2%, things can occur more often or less often.
Go back to this. Let's classify the types of differential equation that we are seeing here. What's the order--is this is a system, is it explicit, is it linear?
The highest derivative is the first one. We have the first derivative, the rate of change, but we do not have a second derivative of a rate of change. Yes, it's the system. We have several equations in parallel. It's explicit. The highest derivative appears alone on one side and nowhere else, and this maybe a surprise, it's not linear. The problem is this term and this term as well. They are multiplying one component of the solution by another component of the solution. This is something like x(t)².
If we take the standard SIR model and add this term 1,000 persons per day in the equation for the rate of change of susceptible persons, which effect could that be-- vaccination, immigration, births, or deaths?
If we increase of change of susceptible persons, that cannot possibly be vaccination. We also would increase the total number with that if you form S+I+R, the total number. You'll see that the rate of change of the total number would be precisely these 1,000 person a day positive rate of change. It cannot be deaths, so we're left with births and immigration. Births do not make too much sense. The total number of persons would increase, but as the total number of persons increases, the number of births should also be increasing. It should not be constant. What's left is immigration, which makes much sense. 1,000 susceptible persons are added to the population per day.
Again, we start from the SIR equations. Now, we include initial conditions, namely that the initial number of susceptible persons should be 98,900,000. We want to have 100,000 infectious persons at the start and 1 million recovered persons at the start. Your task is now to model that 5 million people are vaccinated at the beginning. Here are four options of which you have to pick the right one. A--change the 5 to a 4 in both equations, B--subtract 100 persons per day from the S equation, C-- change these initial conditions to 93,900,000 persons and 6 million persons, D--add 5 million persons per day to the R equation.
Now, comes the answer--option D would mean that we have a huge amount of immigration of recovered persons--no. Option B would mean that the number of susceptible persons is reduced by 100 persons per day that looks like immigration but not like vaccination. With option A, we would be changing this coefficient. If you remember where this coefficient came from, it had to do with the probability of being infected, mainly more precised with the probability that a susceptible person is infected by an infectious person and it had to do with the number of contacts. It doesn't have to do with vaccination. It's option C, we move 5 million persons from the susceptible compartment to the recovered compartment. Of course, we're cheating a little there. These 5 million persons have not recovered from the illness. They are immune as are all other recovered persons.
Let's have another look at the SIR Model to learn about an important phenomenon and that's Herd Immunity. On the right-hand side of the I equation, the factor I appears in every term. If we factor it out, this is what happens. This helps us to investigate what happens when a huge percentage of all persons is vaccinated. We model vaccination by taking persons from susceptible to recovered, which means that the number of susceptible person decreases if the size of the population says constant. When the number of susceptible persons shrinks, the first term in the left parenthesis decreases and eventually if you vaccinate enough people, the complete expression is going to become negative. The rate of change is negative, this means that the number of the infectious persons is decreasing over the cause of time, more susceptible people get infected. The number of susceptible person decreases which means that this first term here gets even smaller which means that the expression in the right parenthesis gets even more negative. The number of susceptible persons can only decrease when enough people have been vaccinated, we can never have an epidemic. This is called Herd Immunity. Not everybody is vaccinated, but if the number of people who are vaccinated is large enough, we can not have an epidemic nonetheless.
Assume that we have a population size of 100 million, which percentage of that population needs to be immune so that we are at the threshold of Herd Immunity and enter an integer number of percent here.
Now, the answer--that's a 60% and here is why. At the threshold of Herd Immunity, the expression in the red parenthesis has to be 0. If we vaccinate even more persons, the number of susceptible persons shrinks further, and this expression becomes negative, and now we can solve for the number of susceptible persons. And we get 1/5 days divided by this coefficient--the days answer and we get the number of 40 million persons--that's 40% of the population of 100 million. If the number of infectious persons is negligible at that point of time, the remaining 60% have to be in the R compartment of the model. They have to be immune. If more than 60% of the population are in the R compartment, this expression is negative and we can never have an epidemic outbreak. So 60% is the threshold of Herd Immunity.
So this is the model that we started with. There's a compartment of susceptible persons, a compartment of infectious persons, and a compartment of recovered persons. Depending on the number of infectious persons, susceptible persons become infectious, and infectious persons automatically recover after awhile. Now, we're introducing another effect latency. Persons who get infected should not immediately become infectious and we model this by introducing another compartment--exposed persons. Depending on the number of infectious persons and susceptible persons, susceptible persons become exposed--that is infected and after awhile, become infectious. This should obviously be called the SEIR model.
Now, that we have the mechanical idea behind the SEIR model, turn that into a set of equations. For each term that may appear in these equations, enter the corresponding letter into one of these boxes. Leave that box empty if it's not needed. 1/1 day, 1/5 days, 5*10⁻⁹ per day and person, the number of susceptible persons, the number of exposed persons, the number of infectious persons, the number of recovered persons, the product of the numbers of infectious and susceptible persons, and the product of the number of exposed and infectious persons.
In our model, the number of susceptible persons always decreases. It cannot have any component here with a + sign and leave this blank. And similarly, the number of recovered persons always increases. It cannot have something with a - sign here. A certain fraction of the infectious persons is taken per day to the recovered compartment. What worked best is 1/5 days, b times well the number of infectious persons f. We need the same here with a - sign. Everybody who leaves the infectious compartment joins the recovered compartment. The similar process is in place with e and i. Exposed persons become infectious afterwards. The rate should be a, 1/1 day seems to be plausible latency of 1 day times the number of exposed persons. Certain fraction of the exposed persons per day gives the exposed compartment and joins the infectious compartment. That same number that leaves here has to joint here, a,e and here we have the infectious process. The number of exposed, infected that is, persons increases per day by 5*10⁻⁹/day times the product of the number of infectious persons times the number of susceptible persons. So, it's c, h and what we gained here, we have to use here, c, h. There's one thing to keep in mind though. This latency period of one day is an average. It's not a strict number similar to this recovery period of 5 days. These are average numbers. The small percentage of exposed infected disease persons will be taken pretty quickly to the infectious state and the small percentage will take quite a while to get to the infectious state. This times of 1 day and 5 days are not rigid. They are averages.
Now, it's time to solve the SEIR model by computer. We've provided the constants that we've been using all of the time. We are providing erase to write the results. We're providing the initial values. And your job is to actually solve this system of differential equations.
One way of solving this is to find out which terms appear in several places-- for instance, this one--the number of persons taken from the susceptible to the exposed compartment and then the number of persons taken from the exposed to the infectious compartment, and the number of persons taken from the infectious to the recovered compartment. If you write it like that, it's easier to read and you can immediately see that the total number is conserved. The number of persons we lose here we gain here, the number of persons we lose here we gain here, and the number of persons we lose here is gain here. So up to very tiny round off errors of the total number and result you'll see that the infection takes quite some time to evolve. The maximum of the exposed compartment comes first and then comes the maximum of the infectious compartment, and as we go, you'll see that the number of susceptible persons shrinks drastically and the number of recovered persons grows by the same amount.
This quiz builds on the one before. Now, the Euler method is in place and we want to experiment with the step size. If you look at the solution curves, you could get the impression that a step size of half a day is way too small. Why should we be stepping in half days. What about steps of say 5 days. Shouldn't that be working given these curves. Let's step size in and leave that step size in for which the solution just doesn't explode and accuracy of 1 decimal place is enough.
This is what happens when we set the step size to 2 days, which look pretty reasonable with the original curves. This is what we get--the numbers explode to 1.4 * 10⁷⁹--that doesn't look too good. In case you're wondering about the conservation of the total number, some of these variables get huge in a positive direction and some of these get huge at the negative direction. Actually, the total number stays pretty much constant at the price of some of these becoming negative. With a step size of 1.9 days, which is the solution, we don't see an explosion, the number stays limited but the behavior is, of course, highly implausible. Now, we look into what's going wrong here.
In the previous quiz, you saw that the solution evolves on a pretty long time scale, but the step size that we are required to use and the solver is pretty short. That's a phenomenon that can occur with certain types of differential equations. It's called stiffness--it's of course annoying if we are forced to use a very small step size even though the solution develops on much longer time scales. There is no generally accepted definition of what stiffness precisely means. But if you think about this situation, you may get an idea about where the term comes from. Imagine that we have a pendulum, a mass attached below and this is a stiff spring and the pendulum swings left to right and can oscillate up and down. Then we have two different time scales, the time scale of the swinging motion and the time scale of the oscillation of the spring and this time scale of the oscillation is going to be much smaller than the time scale of the swinging motion. The observer, however, mostly sees the swinging the motion and can't see the oscillation up and down. This oscillation, however, will force us to use a very tiny time step in a forward Euler solver.
So in which of these three situations would you expect stiff differential equations-- the motion of a spacecraft around a planet, cloth simulation where you, for instance, simulate how one shakes a tablecloth, or a 1-dimensional harmonic oscillator, which strictly moves up and down?
The vital ingredient for stiffness is that the model includes to vary different time scales. That's not the case for this spacecraft and that's not the case for the 1D harmonic oscillator. But with the cloth simulation, we do have to vary different time scales. Imagine this cloth to be composed of tons of tiny, tiny springs. This is actually how cloth simulation is being done. If you push the cloth in this direction, you're working against the stiff spring, meaning you're working on a very small time scale. If you push the cloth up, you're rotating that spring, meaning you're working on a slow time scale. This is a very annoying problem for the simulation of pieces of clothing, in video games, and for special effects in movies. You do not want to use time steps of microseconds for something that seems so simple.
To study what's going wrong when trying to solve stiff differential equations numerically, we'll typically looks at this type problem of very simple differential equation which by the way is not stiff but serves to help us understand what's happening here. The derivative of x with respect to time should minus a constant times x itself. This constant is supposed to be a positive number that's called Dahlquist' test equation and it's easy to figure out what the solution should be. The rate of change is proportional to the current value with the negative proportionality constant. This means that per time step, we are losing a specific percentage of the current value. If we start here, we lose that specific percentage and lose that particular percentage and again and again and again. This is exponential decay, so we know that as t goes to infinity, any solution will approach 0. This decay to 0 takes place regardless of the value of k. As long as k is positive and this proportionality constant, any solution is going to decay to 0 regardless of the value of k. If k is larger, this decay is faster. If k is smaller, this decay is slower. And of course, we want that our numerical solution has the same behavior when we simulate this equation in the computer, we want to see the decay to 0. Technical term here is stability. How small does the step size have to be for the solver to stay stable?
Now we're going to analyze the forward Euler method with the help of Dahlquist test equation. In which range should the step size h lie to maintain stability. Should that step size be smaller than k²/2, k being our constant. Should the step size be smaller than the inverse of that constant. Should it be smaller than twice that constant, or should it be smaller than 2 over that constant.
Let's look at what the forward Euler method does to this equation. X after the first step will be the initial value of x plus the time step times the rate of change at the initial time plus the rate of change at times 0. What's the rate of change at times 0? It's -kx(0). We can factor out x or 0 and end up with this expression. That's the result after the first times step. The initial value gets multiplied by 1 minus step size times the constant that appears in the test equation. Now, let's do the second step--it's x at the end of the previous step plus h times the rate of change at that time or actually our estimate of that rate of change at that time which would be -kx(h), again we can factor out and we can plug in the value of x(h) from the first step, so we end up with (1-hk) times now x(h), (1-hk) x(0) which is if we combine these two, (1-hk)² x(0) and you see how this is going to proceed if you go to (3h), this will become (4h)^4 and so on. The question was whether or not x is going to approach 0. This boils down to whether or not (1-hk)ⁿ approaches 0 as n times to infinity. So what happens to the nth power of (1-hk), we are interested in how large h the step size can be? If the step size h gets large, the expression inside the parenthesis becomes negative, we are subtracting from 1. If we are subtracting more than 1, the inner term here gets negative. So the question here is what happens to powers of negative numbers to some huge power. As an example, what happens if the inner expression equals -2 and we take the 1st, 2nd, 3rd, 4th and so on power. (-2)¹ is -2. (-2)² is 4. (-2)³ is -8. (-2)⁴ is 16 and so on which means that this expression, the nth power would explode and this is precisely not what we want for a stable solution. So let's try a different number. So what would happen if the inner expression was -0.9, and we took powers of (-0.9). (-0.9)¹ is 0.9. (-0.9)² is 0.81. (-0.9)³ is -0.723. (-0.9)⁴ is 0.6561. Each time, we are losing 10% in magnitude but you see that the sign keeps changing, minus, plus, minus, plus, so this indeed decays to 0 but in a strange fashion. It oscillates. You may remember this behavior from the SEIR model that we simulated before. If this expression would equal -0.99, it would still see decay with this sort of oscillation. If it would equal -0.999, this would still work but become very, very slow. So the threshold would be that (1-hk) is in the expression equals -1. The threshold, it's 2/k. So the condition is h is smaller than 2/k for this power to decay to 0 and for the forward Euler method to be stable.
In the quiz, we analyze the stability of the forward Euler method. Now, we introduce some changes to turn the forward Euler method to backward Euler method. This method is called an implicit solver. We are going to see why. The backward Euler method works like this to go from the initial value of x to the value after the first step, it added the step size times the rate of change at times 0. The backward Euler method works with the rate of change at the end of the first step, x(h). Actually, our estimate of the rate of change at the end of the first step. Keep in mind that I'm writing x(h) but I'm referring to my estimates, the numerical solution of the exact solution. This means that the value at the end of the first step appears on the left-hand side and on the right-hand side. This is why this is called an implicit method, we have to solve for x(h). It's hidden in this equation. It's not explicit so let's solve for x(h). When this term over to the left-hand side, then we have x at the end of the first step plus time step times constant times x at the end of the first step equals the initial value, we can factor out x(h) and find that x(h) equals the initial value divided by 1+hk. This is something like 1/1+100 if you choose a really large h which would be 1/101. If you form powers of these, they are going to decay pretty quickly and there is no oscillation whatsoever, so we have no problem at all with this expression. It's always going to work no matter how large the steps as is. This is very handy when one tries to solve stiff differential equations.
Your job is now is to implement the backward Euler method for the SEIR model. We provide the forward Euler method as a starting point. So this end, we have reshuffled the order--EISR to make things easier. We can start with the R equation using what you have computed for E and I and S before and then you can continue with S and with I. Doing the backward Euler method for E, it's a little harder. We provided the solution in case you don't want to spend that time. If you start from the bottom, you can check your result whether or not it make sense by running the software for that single line. When you're done, set the time step to 5 days to see that this method can handle such a long time step.
Let's start with the output of that software. For reference, this is the result of the forward Euler method with a step size of half a day. You already saw that the Euler method had step sizes of 2 or more. But there's a price to pay for the stability of the backward Euler method. The code gets really complex as opposed to the code of the forward Euler method. You really do need pencil and paper to derive these equations. That's what we're going to do now. The equation for R, the number of recovered persons, is the easiest one. We'll start with that and to save space let me write under score 1 for the current quantities and subscript 2 for the quantities of the next step. The forward Euler method would turn this into the following. The value of R one step later is the current value of R plus step size times the current rate of change, which is 1/5 days times the current value of I. I₁ in my new notation. This would be the forward Euler method but now we're concerned with the backward Euler method, which means that we have to use the rate at the end of that step meaning I₂ has to go in here. And now we're done. We know the current value of R. We have computed the next value of I in the lines before, and that's all we need to compute the next value of R. Next up is S. If we were using the forward Euler methods, this would work as follows. The next value of S is its current value plus the time step times its current rate of change meaning -5x10⁻⁹ per day and person times I at the current step and S at the current step if we were doing forward Euler. But now that we're doing backward Euler, this has to be I one step after that and S one step after that. We bring this down to the left hand side because we're interested in S₂ and we factor out S₂. So the factor of 1 times S₂ and now it's plus this coefficient times I₂ times S₂. 5x10⁻⁹ divided by day person I₂ and all that remains on the right hand side is S₁ and now that is simple to solve. S₂ equals S₁ divided by the very same expression that we've got here. I'm lazy. I don't spell it out. This completes the computation of S₂ for which we need I₂. But we have computed that in the line before. And now the equation for the rate of change of I, the number of infectious persons. Those persons that is infected persons become infectious with a certain time constant and infectious persons recover with another time constant. If we were implementing forward Euler, it would be doing the following. I₂ the value of the next step equals the current value plus the time step times the current rate of change, which is 1/1 day E₁ minus 1/5 days I₁ but now we want to implement the backward Euler method, which means that we have to use the values for the next time step and this rate of change. I bring this term with the I₂ over to the left hand side and factor out the I₂. This leaves me with one I₂ plus h times 1/5 days. On the right hand side, I₁ remains and h times 1/1 day and E₂ remains. And now we can divide by this expression and find I₂ equals the right hand side divided by this expression I₁ plus h times 1/1 day times E₂ divided by this expression. This completes the equation for I. The main message here is that implicit methods come at a price. Solving these equations can get very ugly. If you are convinced at this moment, you may want to skip the rest of the segment. Otherwise, stay with me to see how we can compute E₂, which by the way is needed to be able to compute I₂. A trick to find the value of E for the next time step is to look at these two differential equations and parallel. For the first one, the backward Euler method says the value for the next time step is the value for the current time step plus the time step times our estimate of the rate of change for the next time step which is 5x10⁻⁹ per day and person times I₂ S₂ backward Euler minus 1/1 day E₂ backward Euler again and not E₁. And for S the new value is the old one plus the time step times minus our coefficient 5x10⁻⁹ per day and person times I₂ S₂ backward Euler again and that's it. If you now form the sum of these two equations, you'll get E₂ plus S₂ equals E₁ plus S₁. This ugly expression cancels with that ugly expression, and we are left with h times minus 1/1 day E₂. So this is what we have reached by now. E₁ the current value is known. S₁ the current value is known. We want to determine E₂, the next value. At this point, we do not yet know S₂. But luckily enough, we derived an equation for S₂ to be for namely S₂ equals S₁ divided by 1 plus the step size times 5x10⁻⁹ divided by days and persons times I₂. And we have an equation for I₂ that says I₂ equals I₁ plus step size times 1/1 day E₂ divided by 1 plus step size times 1/5 days. If you use this equation to express I₂ and then insert this result as S₂, you end up with an equation that contains E₂ and E₁ and S₁ and it can be solved for E₂. E₂ appears in three places and with a little edge ??? you can turn the resulting equation into a quadratic equation which is the standard solution formula for quadratic equations. The result is what you see in the code. Hopefully, this is now really convincing that working with implicit servers requires quite some effort with pen and paper.
Let's summarize facts about implicit solvers. The simplest of which is the backward Euler method. The big benefit is the stability. We can use large step sizes. And stiff differential equations don't force us to use overly tiny step sizes, which results in speed. Some types of simulation are almost impossible without implicit solvers. Think of cloth simulation and the simulation of water vase. These types of simulation may be forbiddingly slow without the speed offered by implicit solvers. The big disadvantage of implicit solvers is that they are well implicit. You get equations that we have to solve, which is hard to do or can even be impossible to do with standard functions. Mathematicians have come up with quite a range of methods to address this issue. For instance, in a predictor-corrector approach, you use an explicit solver such as the forward Euler method to predict the value for the next step and then use an implicit method to correct it, actually Heun method, that we learned about in a previous unit is of that type. All you can use is a rate of methods to numerically solve this implicit equation, not with pen and paper. And the second disadvantage of implicit solvers, most of them tend to lose energy in physical simulations. We saw that the forward Euler method, an explicit solver, increased the energy. Implicit solvers tend to work the other way round and to lose energy.
Now, let's have a look at the zoo of methods that we've seen so far and then let's add one further method to that collection. The forward Euler method work by starting from the current value and then incrementing that by step size h times the rate of change at the current point. I'm writing f for the rate of change here. Heun method also known as the improved or the modified Euler method advances by the time step times the average of the rate of change at the beginning plus the rate of change at the position predicted by the forward Euler method. With the forward Euler method, the error of a numerical solution grows linearly with step size. If you double the step size, the global error also doubles approximately at least. The forward Euler method is a method of order 1. This function, the error grows like h¹. Heun method, however, is the solver of order 2. The error grows like h squared. If you take half the step size, the error is shrinking to one quarter, which is much more efficient. For the backward Euler method, you advance by the rate of change at the new position, which makes this equation difficult--it's an implicit equation. To complete our zoo of methods, we add a method that looks like Heun method but is implicit like the backward Euler method. This new rule is called the trapezoidal rule. Like Heun method, it's a method of order 2. And if you compare these equations, they look pretty similar. The trapezoidal rule advances by the average of the rate of change at the beginning and the rate of change at the end, which means that this is an implicit equation.
Now apply the trapezoidal rule to Dahlquist test equation. How large can we chose the step size and still have numerical solution that approaches 0, x(t) goes to infinity. Do we have to chose a step size that's less than 1/k, k being the constant of Dahlquist test equation, or a step size that's smaller than 2/k or a step size that's smaller than 1/k², or could we pick any positive value for the step size h.
Let me summarize some observations that we made in this unit about the complexity of models. For instance, in the SIR and SEIR models, we worked with fractional population numbers. We didn't care whether or not the numbers of persons are integers. In the SEIR model, we used an average amount of latency. We did not build in a fix delay and we assumed perfect mixing. There is no geographical distribution, there is no age distribution, there are no isolated towns. One of the advantages of the simple models we looked at in this unit is that we were able to derive general statements such as the existence of Herd Immunity. Of course, we have to be careful. A general statement about a cartoonish model may not be that helpful, and the simple equations that we've got from these models are helpful in may ways. For instance, when it comes to working with implicit methods such as the big foot Euler method.
In this unit, we look at models with flows between different compartments. The obvious use was in epidemiology where we had susceptible persons who became ill and eventually recovered, but what about economy. Money may flow from one account to the other or goods may flow from one country to the other. In chemistry, you may have several types of molecules that combined to form other types of molecules. In biology, the members of a species are born young, become adolescent, and then mature. In the next unit, we are going to look at examples from biology, in particular, to study equilibria that is situations in which one falls is exactly balance by another falls.
In this problem, you're going to use the SIR model to look at a situation in which people's immunity doesn't last. In other words, after someone who has been infected moves to the recovered population. Here, she can become susceptible again. We're going to call this waning immunity and we'll set the constant over here waning time to equal the amount of time that it takes for a person's immunity to fade away. Most of the information we've given you in the supplied code you've seen before, but there are couple of things that are particular to this unit. You'll see another transmission coefficient, which is a measure of how frequently each day, each infectious person might infect someone else. You could also think of it as the likelihood that an interaction between an infected person and a susceptible person went up with the susceptible person becoming infected. Just make sure that you know that its unit is 1 per person per day and the first thing that we want you to do is to define waning time right here. So eventually we have a system in stead state where the number of recovered people is equal to 2 times the number of infected people. Think about what this means in terms of the time that people should spend infected versus the time that they should spend recovering. Just for reference, I've include the equations for the standard SIR model right here. Remember that this equation don't include waning immunity. First think about how you need to change them in order to reflect the situation at hand. Once you figure that out, you'll need to translate your mathematical expressions into code. We've already created list for S, I, and R for you. And we set the initial values for susceptible, infected, and recovered populations right here. Next comes a for loop, part of which you'll have to fill in. We've given you S to I, right here, which is a number of people moving from susceptible to infected during one-time step. We've also given you I to R, which is a number of people going from infected to recovered in a given time step. What you need to do is to define R to S. The number of people who in one-time step become susceptible again after they've recovered. Just to be clear, all three of these are measured in number of people, not in number of people per time. Make sure that you check that your units are correct. Your last task is to modify the recursive functions right here for S, I, and R.
Now, let's look at the solution for a problem concerning waning immunity. First things first, the constant waning time should be defined as 2 times the infectious time. We can show this mathematically using the equations we have for the derivatives of S, I, and R with their respective time. We know that after a long period of time you want to attain a steady state situation and there's a number of people in each portion to the population-- susceptible, infected, and recover stays constant. Since they want to find out how long people should spend in the recovered stage, we start with the time derivative of R and set that equal to zero. Since you know that Rdot now has an extra term added to it or actually subtracted from it showing the number of people that are leaving the recovered population and going back to the infected population. We can set these to terms right here equal to zero as well. Just to note, I've used CINF to stand for the infectious time and CWAN to stand for the waning time. Now with just a little bit of Algebra, we come up with the answer that R=2I. Since we know that we want the number of recovered people to be twice the number of infected people, we can plug in this extra information to the equation above and end up with the answer that the waning time is equal to twice the infectious time. This put us directly into the next part of the problem. We defined R to S in the same way that I to R is defined except that we replace infectious time as the waning time and I step with R step as you can see right here. Then moving on to our recursive relations, for each element in a step plus one position, we need to take into account the value of the previous element in the number of people added and subtracted from the population during each time step. We know that the one thing that has changed in this model from the standard SIR model is that people are now moving from the recovered population back to the susceptible population. This is why we've needed to add in this extra term R to S, which we subtract from R and add to S. Now, let's run the program and see what we get. Here we see we end up with this fancy graph, which has three different series for the three different parts of the population we're looking at. Remembering how we set the initial values for S, I, and R. Remember that the blue line here stands for the susceptible population. The green line stands for the infected population, and the red line stands for the recovered population. And the maximum time that we're looking at 60 days, you can see that the red line is graphing twice as many people over here as the green line is marking, which is exactly the answer that we wanted to end up with. Congratulations on successfully completing the first problem of Unit 3.
Welcome to the second problem of Unit 3. If you thought that the whole situation has started to look grim in the last problem when we dealt with waning immunity, this time circumstances had become even more dire. There's no immunity at all. And so with the typical SIR model, people simply switch back and forth between susceptible and infected. We're assuming that more people are lost in the scenario. So that means that the total number of people which I'm going to call N stays constant. Since every person is either in the S group or the I group, we know that in any given moment S+I=N. Here once again are the standard SIR model equations. You'll need to change this which reflect the new situation similar to what you have to do in the last problem. Start by reducing this system at differential equations to a single equation that contains only I and N but not S. Then solve this equation numerically using the trapezoidal method, which in case you forgot I will show you right now. Remember that the trapezoidal rule finds the value of a quantity after an additional time step by taking their original quantity and adding to it h times the average of its rate of change at the two points in consideration. Now use notation here that probably make sense but just to be clear x sub-n is just equal to x measured at the n time step and x do- n is just equal to the time derivative of x at the n time step. Now once you go through using the trapezoidal method, I want you to solve for S and then convert your equations for I and S into Python. Right now the supplied code contains expressions for S and I that come from the forward Euler method. Before you fill in the program with the trapezoidal rule run it as it is right now. Try it first with h set to 10 as it is here and see what happens then change h to equal 0.5, run the program, and compare the difference between the plots. When you change the code as we've asked you too, see how each of these step sizes works with the trapezoidal method instead. After you're done with everything, set h to 10 and press submit.
Now for the solution to the second problem of Unit 3, let's start it by adding in a constraint to the system with this line right here about n. This shows, I've explained earlier, that the total population n remains constant in all times. The people just switching between the S and I groups. Now let's move on to the core of the problem. Switching the contents to the for loop to use the trapezoidal rule instead of the forward Euler method. We begin with the standard expressions for Ṡ and İ, but this time showing that people can move from the infected population back to the susceptible population directly. Just for simplicity, I'm going to call our constants a and b right now and I'll change them later to show what's actually in the code. Now we can get a single equation that encapsulates the information in all three of these equations if we just substitute N-I for S in the equation for İ. This is what this line right here shows. If we do a little bit of rearranging and simplifying, we get a quadratic equation for I-dot. Now this make things a lot simpler altogether because now we have one equation and one independent variable instead of three equations in two independent variables. Now it's time to move on and use the trapezoidal rule. You'll notice here that I'm using a t in the subscripts instead of an n this time. That's just so that we don't get confused between big N's and the little n's. So we know that for the t+1 time step, I is going to equal this expression just in the trapezoidal rule. However, we just created an equation for I-dot that depends only on I. So we can plug that expression in to get this expanded equation right here. Since many of the terms in this expression have factors of h/2 and a outside of them, I'm just going to pull out a factor of ha/2. With a little bit of rearranging, we come up with this long expression. Now we can think of this equation as made up of two parts-- one part that involves I sub-t+1 and another part that just involves I sub-t. Since we want to think of I t+1 as the only independent variable, we can lump everything else under other variable names. Looking at the equation, we have one term with (I subt+1)² . A couple of terms with I t+1¹ and several terms that don't have I t+1 at all and we're going to call everything that's in front of I t+1¹p. And we're going to call everything that doesn't involve at all q and this is going to give us the simplified quadratic expression that shows very clearly how this is a quadratic equation using I sub-t+1. Using the quadratic formula and rearranging a little bit, we get this expression right here using p and q to express I sub-t+1. Now this actually has two possible solutions-- one with a plus sign here and then one with the minus sign here. So how do we decide which one to pick? Well we know that this h to times that goes to zero. We want I sub-t+1 and I sub-t to get closer and closer together since the curve is going to become smoother and smoother. Since p has a term of 1/h inside of it, this means that p is going to infinity as this happens. If we pick the minus sign over here, they're going to have minus infinity, which means that I-sub-t + 1 is definitely going to go to negative infinity and this is definitely not what we want. So just to be safe, we will pick the plus sign here. Now we return to the code where we've translated our expressions for p and q into Python. To do this, we have to remember that the coefficients a and b that we were using when we're writing things by hand actually each have different names in the code. A is equal to a transmission coefficient and b is equal to one over the infectious time. Using this information, you've just done a simple translation of p and q into Python. This makes writing out the expressions for I step+1 and S step+1 very easy. For I step+1, we just write in the solution that we got from using the quadratic formula using of course the plus sign right here in the variables p and q since they are just written right here anyway. Because we know that I plus S is always equal to n, our expression for S is simply equal to n-I at the time step. Now, let's run this first with step size 0.5 and see what happens. This looks very nice. We have two really smooth curves with the infectious population is green and the susceptible population is blue just like in the last problem. We see that since each curve in the plot approaches a constant value at long times we're approaching a steady state situation. Now let's run the program instead with step size 10. So this is what we get from our model if we're using h=10 instead of h=0.5. Even though now that the curves is particularly smooth, we can see that they obey the same overall pattern as we saw in our graph using the smaller step size. The most important thing is that we still are approaching a stead state situation at the end of the graph. Now think back to what happen when you run the supplied code like we ask you to. When you run the supplied code with h=0.5, you should have gotten a graph that looks like this. This definitely looks a lot like what we saw when we run the solution code using h=0.5 as well. However, let's peek at what happen when you use h=10 using the forward Euler method. You can see that there are a lot of differences between this and the graph that we just saw and also between this and the graph with the solution code when we had h=10. The most important thing to notice though is what's happening over here in this region. You can see that instead of each curve approaching one set value as time goes on they keep crossing paths at each other. In other words, since I and S here don't obey nearly the same trend that they do with smaller step sizes, this shows that the forward Euler method doesn't allow for stability with large step sizes nearly as well as implicit methods like the trapezoidal method do. Hope that was an interesting investigation with the difference between implicit methods and other kinds of methods and get ready for the last problem of Unit 3.
Last problem of Unit 3, we find a population faced with a malaria epidemic which is definitely very scary for everyone involved. Luckily, we're here to help them out. In order to slow the spread of the disease, we're going to make people start using mosquito nets in their houses, so they can avoid being bitten as often. We're going to introduce the nets after 100 days and they immediately decrease the probability of a person getting bitten by 90% as we can see from the constant that we've defined right here called bite reduction by net. Bite per day in mosquito tells us that every person gets bitten by each mosquito, 0.1*per day. Now when a mosquito bites a human, there's a chance of two different things happening that are a particular interest to us. Now, if a mosquito has malaria, it can give a disease to human or if a human has malaria, then he or she can give it to the mosquito. The probability of either of these things happening is shown right here. Once infected, it takes a human approximately 70 days to recover. At this point, it's important to notice that mosquitoes only live for 10 days. Clearly, humans also have a finite lifetime but compared to mosquitoes, we might as well be immortal even if we do end up with malaria. Sets that we're going to pretend for the purposes of this problem, that we can model a scenario like this, we need two separate bubble diagrams-- one for the mosquito population and one of the human population. The mosquito life cycle is kind of a sad one--as soon as they are born, they fall into the S part of the population, the susceptible part. From here, they can become infected and then eventually they die. There's no option for them to become susceptible again. A couple of other things to keep in mind about mosquitoes. For one thing, they can't get malaria until after they're born. Second, only female mosquitoes can transfer malaria. So just pretend that all the mosquitoes that we're talking about in this problem are female. Now this problem requires you to take into account all the factors affecting the infected population and to include them using the forward Euler method. Think critically, keep track of your units, and good luck. If you have any trouble at all, visit the forums.
The first thing that we need to deal in this problem is the issue of the mosquito nets. Well we know that at first there are no mosquito nets in place at all, so we're going to create this variable called net factor and initialize it to 1. We know that after 100 days, net factor is reduced by bite reduction by net, this constant that we created up here. We add in this if statement right here saying that after 100 days, net factor is reduced by 0.9 which is bite reduction by net. Remember that this is basically the probability that a person is going to be bitten. It makes sense that if you add in the mosquito net, the probability of you being bitten will go down. Now for implementing the forward Euler method with the infected population, both human and mosquito. For the number of infected humans at the next time step, we of course, need to start with the number of infected humans at the previous step. To this, we're going to add on the number of people that are being added to the infected population. We are going to subtract the number of people that are moving out of the infected into the susceptible population. Let's first think about the term that's being added to the number of infected humans. We take the probability of a person being bitten which is in that factor times the number of bites per day and per mosquito that a person gets also multiply that by the number of infected mosquitoes and then we multiply this by the fraction of the total human population that is still susceptible also multiply this by the probability of a mosquito transmitting malaria to a human. On the other side of the coin, we have factors that move people away from the infected population and back into a susceptible population. This is a very simple term. We simply take the number of people that were infected at the previous time step and divide this by the time that it takes a person to recover. Remember that all of that was multiplied by age and added to the number of people infected at the previous time step. Now, let's take a look at how the infected mosquito population was changing. Since we're using the forward Euler method, we started with the number of infected mosquitoes at the previous step and then add to this age times a bunch of other things. Now to find out how many mosquitoes are moving into the infected group, we start with the probability of a mosquito biting a person, which of course might give the mosquito malaria as well, multiply this again by the bites per day of mosquito, multiply this by the number of susceptible mosquitoes and the number of infected humans. We then divide this by the total number of humans and multiply the transmission probability from human to mosquito. If you remember the drawing that I showed earlier on in the video introducing this problem, I told you that the only way the mosquitoes can no longer be infected is to die. They don't move back to being susceptible. This means that the only thing that is going to decrease the number of infected mosquitoes is death of mosquitoes. So, we subtract from this entire expression if 1 over the mosquito lifetime times another infected mosquitoes at the previous step. I know that that was a lot of constants and variables to keep track of, but great job if you're able to fill all this up. Now, let's take a peek at what our final graph for Unit 3 looks like. Now in this graph, the blue curve shows the fraction of humans that were infected and the green curve shows the fraction of the number of mosquitoes that are infected. At first glance, it might look like there are many many more humans infected the mosquitoes. But remember that because this is fraction, this doesn't actually show us plain numbers. In fact, we started with 10¹⁰ mosquitoes, but only 10⁸ humans. So even if a much smaller percent of the mosquito population was infected from the human population, the mosquitoes are able to infect such a great portion of humans because there are so many more mosquitoes than there are humans to it. You can see that the number of both humans and mosquitoes infected, at first, increases sharply and then levels off and stays at a very high number. However, as soon as 100 days hits, and we introduce the mosquito nets., the number of infected in each population starts to decrease certainly. We can see that a mosquito net idea was definitely very effective for helping slow the spread of malaria in the human population as well as the mosquito population. Great job on all three problems in Unit 3. I hope you found this discussion of the SIR model and its variations interesting.