
Title:
0336 Seir Backwards Solution

Description:

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.