Welcome back.
We're gonna switch gear a little bit now
and move away from sort of mathematical treatment
of this continuous model
and think about how to implement them computationally.
So as I said in the previous video,
one of the big strengths of differential equations
is that they give you, without any additional alteration,
a complete view of the flow diagram
because that's literally
what those differential equations are
is the flow or the speed of the system
at any point in variable space.
So if you tell me the state, x of a system,
the time derivative of x of t gives me the direction
of the movement of the system, as well as its speed,
how fast is it going there?
I have three different notations right now.
So you're probably most used
to the explicit time derivative.
As I said, I like using Newton's notation.
It's just a little simpler,
this dot representing a time derivative,
and really the way we write those differential equation
is having x dot or s dot or i dot, whatever it is,
equal to some function of the state of the system x,
maybe also explicitly of time, t.
So these differential equations are always written
as x dot equals f of x of t.
That's our differential,
our system of differential equations.
So in this video, what I'd like to describe is
what you'll be doing often,
which is not necessarily mathematical analysis
over these differential equations,
but just tracking them in time.
So you put the system in some initial state,
you have some initial conditions, x at t equals zero,
and you wanna follow this flow of dynamics.
So let's say that our system here
has only one degree of freedom.
This x of t is just one variable,
and we're interested in knowing from x0.
An initial condition of T equals zero,
the state of the system for all t, for all time.
Now there might be some true evolution for this system,
which might look like this,
but that would be the true x of t.
The simplest way to follow it numerically
is discretized time, the same way we're doing
in discretized model.
So at t equal zero,
you can evaluate your differential equation, x dot of t.
So you evaluate this function, f of x of t here,
and what it gives you is the instantaneous derivative
at that point.
So then one thing you can do once you have that is say
x at time, t, plus h where h is some time step,
simply gonna be equal to where I am
plus the projection over h of the slope,
f of x of t, my differential equation.
And then that might put me here.
And then I reevaluate my equation,
project myself, reevaluate project, project, project.
And as you can see,
I don't expect to get exactly the right distribution.
I mean, I'm doing a little bit of whatever here,
but the point is I expect to get something
that's close to the true solution, but that deviates.
So this method, the simplest method
is called Euler's method.
And we're gonna say that its local error
is of order, it's a big O here, order h squared.
The reason for that is we know that there's a true x of t,
and then there's the x of t from Euler's method.
And we can write x of t plus h
by developing around x of t,
the Taylor series of the true solution.
So we know that it's gonna be x of t plus h
f of x of t plus...
While we're doing that step,
we know that the system continues evolving,
continues following the dynamics.
So there is some curvature,
curvature f prime of x of t plus h cubed.
If you remember your Taylor series,
three factorials, f prime prime of x of t, and so on.
And so the reason why we say that Euler's method
as a local error term in h squared
is because it captures those first term
of the Taylor series,
but ignores the higher order term,
so ignored by Euler.
And so at every time step h that we do,
the biggest error term is gonna be this one here.
So if the curvature of f prime is equal to zero,
if we have a linear growth of a system, let's say,
so f prime is gonna be zero.
Euler's method is gonna do awesome
because this error term is equal to zero,
but for any like non-zero error,
really the key point is that
the biggest term that we have is proportional to h squared,
and that means two things.
That means just as we move from discrete to continuous.
If h goes to zero, if we have a very good computer
and we're very patient and we say I'm gonna do,
you know, infinitely many step
of tiny, tiny, tiny, tiny durations,
then Euler's method is gonna be fine,
but otherwise, you know, it's important to remember
that you're going to have those local error terms.
And often, the difference here, let me put it in blue,
so the difference here is the local error.
And often as in my cartoon, this error sort of grows
with time steps, right?
It might shrink sometimes.
You might be lucky,
but overall you expect it to grow.
And so, we're often gonna say that the global error
of our numerical integrator is the sum.
Let me rewrite that more clearly here.
Global error is gonna be the sum
over all time steps of the local error.
And we had the fact...
So we had this idea that the local error
was an h square for Euler.
Well, this is gonna be equal
to the duration of our numerical integration.
So we're integrating from t equals zero
to, let's say T, capital T.
So how many steps are we gonna have to do?
How many local error are we solving here,
are we summing up?
So it's capital T divided by h.
That's the number of steps that we have to do,
times local error, which is of order h squared.
So this is gonna be an h.
And what I'm trying to say here is that
there's this key thing is,
which is that Euler's method is capturing up to order h
the Taylor series of the true solution,
and it's so that its local error
and what's being ignored is h square,
but the global error,
because you're doing a lot of different steps,
the global error is still gonna be of order h
even though you're capturing the thing.
You're capturing the term of order h.
Really, I wanna highlight this distinction
between local error and global error
because it is kind of surprising to have something
that captures the derivative.
And yet as you integrate a trajectory in time,
you do expect its error to scale
with your step size linearly.
So basically if you wanna do a numerical integration
that's 10 times closer to the true solution,
if I wanted my red that's here to be 10 times closer,
I would need to do 10 times more steps.
But we can have integrator that do better.
So, I apologize here.
I expected a blank sheet.
So we wanna integrate again from t equals zero
to t equals T,
whoops, x of t.
We have something different, true x of t here.
And now we want to do better.
Well, what's the key part that Euler's method is missing?
It's the fact that as you do one step,
your derivative is already changing
during that step, right?
There might be a second derivative.
And so it's a mistake to to project it too far
into the future.
So there is another method, which I've learned
under the name of improved Euler.
Improved Euler's method.
If you Google around, in English,
I believe it's often called Heun's method.
Whoops.
Too sensitive here.
Heun's method.
And this, I'll tell you right away
is a second order method in global error.
So this method is gonna have a global error
that scales of the order of h squared,
meaning that, you know, if you wanted Euler's method
to be 100 times more accurate,
you needed 100 times more time step.
Here, you would only need 10 times more time steps
because the global error scales better with h.
And so that means that the local error
for the same reason as before has to be an h cubed
because the number of steps we're making grows
as one over h as we decrease our time step.
So the local error error and the global error,
there's one difference in order, okay.
Well, what is this method?
Let me just fix my typo here.
Heun's method.
The method itself looks super simple.
What we're gonna say, and I'm going to write it
as xh for Heun at time t plus h,
it's gonna be equal to x of t
plus h f of x of t
plus f of xe time t plus h.
So I'm switching up my notation here a little more.
So let me explain what I mean,
so that the prediction from Heun's method
is gonna be a projection from x of t over a time step h
of the average between the measured derivative.
So this term here, this f of x of t
is the initial derivative here.
And now what I have, this f of xe, let me go back,
is gonna be the prediction from Euler's method.
So basically what you do here is you say
I look at my initial derivative,
I project myself over time h.
Okay, project myself over time h,
so I get an Euler prediction of here.
Well, now at this point,
my derivative here is much faster,
much steeper than it was at my initial point.
So I'm gonna take the average
of both, my first point and my second point.
So my average might look something like this,
and now this is what I'm gonna project,
and I'm gonna end up here.
So basically every step of Heun's method
does one step with Euler's method,
looks at the derivative and average.
The derivative that Euler's would have
for its second step with the initial step
and then does a true first step under Heun.
Let me give another example here.
So with Euler, I evaluate my derivative here.
I project into the future up to 2h,
which would be, you know, somewhere here.
Euler is gonna end up way off base here.
It's a steep derivative, might end up here.
Now would evaluate the derivative here, which is negative.
So Heun would average both,
which would be pretty much flat in this case
and end up here.
Doing much better than the blue dot, right?
So Euler is doing this and this,
and Heun is doing this and this.
And the key part is that the local error
and the global error are both one order of magnitude better.
So if you want to be 100 times more accurate,
you only have to do 10 times more steps.
Now every single evaluation takes two evaluation
of your derivatives.
So every step does take twice as many operation,
twice as many evaluation of this f function
for your derivatives.
But really you're still winning, you know, a factor of five
if you wanna be 100 times more accurate.
So you can have your time steps be 10 times smaller.
You do twice as many operations.
So that's 20 operations for your 100 fold precision,
which is still much better
than what you would do with Euler.
Now as you'll see in the homework,
Heun's method, which is this definition here,
is actually not defined as a simple arithmetic average
between f of x of t and f of xe of t plus h.
The basic idea looks like you're just taking an average,
but that's not what it is.
The goal here, and this is a hint for the homework,
is that you can ask, well, what is f of xe
at time t plus h?
You can develop that with the Taylor series.
And then you can show that this one half is actually chosen
such that this xh approximation
at time t plus h captures one more term,
so captures the h square term
in the original Taylor series for the true solution.
This is something that like sounds weird.
Like your only evaluating first order derivative.
How could you like go and get information
about the second order derivative?
But if you do work it out by using Taylor series
of both x of t, which we did right here,
so this is our Taylor series for the true solution,
and then you go do your Taylor series for this term,
for the derivative itself under Euler's method,
then you can see that you're able
to combine multiple evaluation of your derivatives
in a smart way to capture more and more
or capture higher order terms.
In a future video, so next week,
I'll go into deeper detail in the in the Runge-Kutta
or Dormand-Prince method,
which captures up to order four terms,
and there it's just gonna be a lot of arithmetics,
but we're gonna be able to see that similarly
by doing four evaluations or derivative,
we can get like a higher order term,
and there it's gonna look nothing like an average.
So the fact that you have a simple average here,
this divided by two for both evaluation
of the derivative is sort of coincidence,
but it's a nice explanation.
And as you can see, like in any cartoon,
like the one in the top left here,
it's clear that it's gonna do better.
So, apologies, by the way, if my picture here
was, in the way for some of the words up there.
But, you know, this is a topic
that we're gonna revisit time and time again
because the key part is that even if you ask,
you know, a scientific library,
please integrate my differential equations through time,
it's super important to know what order,
how much error are you getting,
what's the order of the method that it's using
by default or that you're asking it to solve
because what comes out of discrete computer
and discrete integration is not the true x of t,
but an approximation of it.
And as I showed in the cartoon,
those approximations can be way off.
So it it's gonna be important to distinguish error
due to our model specification
and error due to how we solve the model numerically.
And one way to test how well we're doing
is gonna be to compare our predictions
from the mathematical models
with simulations of the results.
So in the next video, we're gonna see how to use
some of the tools we quickly saw
to simulate discrete models,
but this time in continuous time.
So I'll see in the next one.