Tip:
Highlight text to annotate it
X
Welcome to Calculus, I'm Professor Ghrist.
We're about to begin lecture 48 on numerical O.D.E.'s.
Let's retreat from the digital back to the analog world and see how the two relate.
We have used our intuition for smooth calculus
to tell us things about the discrete world.
In this lesson, we're going to see how the discrete calculus helps
us solve smooth problems, in particular,
the problem of solving ordinary differential equations.
In this chapter, we've been using our
intuition from differential calculus to understand discrete calculus.
But the inference can run backwards as well, and in this
lesson we're going to use our discreet understanding to solve problems
in differential calculus.
In particular, we're going to consider differential equations.
Recall some of the main classes.
In the simplest case one just integrates both sides of the
differential equation. Others are not so easy, but if
the equation is say, separable, then again, we can integrate and solve
in more complex examples still, say,
a linear differential equation, we can apply an integrating factor.
This is great, but it will not help us solve the general
case, and indeed, at some point, calculus fails to help and
we have to turn to a numerical method for solution.
We're going to focus on the initial value problem for an ODE of the form
dx/dt equals f of x,t.
We can think of that differential equation as defining a slope field in the xt
plane, that tells us how x changes with t, or at what rate.
Now given some
initial condition, x not, at an initial time T
naught. What we want to do is approximate the
final value, x star, at some final time,
t star. Now, we don't know the solution x of t,
and we're not going to be able to use calculus to get that analytic solution.
So what we're going to do is try to
approximate that final value, x star, as best we can.
To approximate, we're going to use sequences.
Let's begin, first of, all, along the time axis.
We're going to write a sequence, t sub n, that
goes from t0, t1, all the way to t, capital n,
where t not is the initial time, and t capital
N is our final time, t star. This is going to be a uniform
grid, meaning that we have a step size, h, equal
to the final time, minus the initial time, divided by capital
N. We are likewise going to approximate
the x values by a sequence x sub n where
x not is the initial condition on x and hopefully
our final x value x N is close enough to x
star that we have a good approximation.
The simplest way to do this is called Euler's method.
Perhaps you've seen it before, but we're going
to look at it from a discrete calculus perspective.
We're going to take our differential equation and discretize it using
sequences, and instead of derivatives,
differences, a literal quotient of differences.
What does this mean? While if we consider at
the nth time step that on the right we have f evaluated at
x sub n comma t sub n, on the
left we have the quotient of differences. Delta x is
x n plus 1 minus x n. Delta t is t n plus
1 minus t. And of course,
since we're using a uniform grid, that delta t is simply equal
to hr step size. And now, with a little bit of
rearrangement, we obtain Euler's method,
this is an update rule. It tells you that if
you know X sub N and T sub N, you can obtain X sub N plus
1 as XN plus H times F evaluated at XN coma TN.
That's the method, but what does it mean geometrically?
Well, if we are at time step T sub N and we have our approximation X sub N, we'd
like to get to an approximation TN plus 1. We don't know
the true solution that passes through that point, but we do
know the value of F there. That is telling us the slope, and what
Euler's Method is really doing is projecting forward 2 time T
sub N plus 1. At heart, Eulers Method is linearization.
So why does it work? Well,
it works because we repeat this at each time step, linearizing and approximating.
As we move from timestep to timestep, let's look at a simple example.
A differential equation that we already know how to solve, dx dt equals x.
We'll choose an initial time of 0, an initial value x knot
equal to 1. A final time T star equal to 1
is X star. Now we already know that the solution of
this equation is E to the T, so X star really should be E,
but let's see what happens when we chose a step size of H,
that is final time minus initial time over capital N.
In other words, 1 over capital N. Euler's Method
tells us that xn plus 1 equals xn plus h times f.
Add xn, tn. What's f in this case?
Oh, it's the right hand side of the ODE. In this case, just the x value.
So we can simplify Euler's update
rule to xn plus 1 equals
xn plus hxn. Factoring out the xn, we get quantity
1 plus h times xn. Now we've seen
this before, this is really a recursion relation, E applied
to the sequence X equals quantity 1 plus H
times X. Now, you may recall from our lectures
on discreet calculus that we know how to solve such
simple recurrence relations. Sequence X is really
X naught times quantity one plus H to
the N. It is that coefficient of 1 plus H that
determines the solution. This means that our final value X sub
capital N is X naught times quantity 1 plus h to the capital n.
If we substitute in 1 over capital n for h, then we obtain something
that should look familiar. Quantity 1 plus
1 over capital n to the capital n power. You may
recall that for n large, this is a very good approximation
to e, our desired final value. What
happens when we do an example, the answer to which we don't
already know. Consider dx, dt equals cosine
t plus sine of x, with an initial time.
Zero, an initial value of zero, and a final time of pi.
I don't know what the final value should be.
This is not a separable equation, and so we're going to fire up the computer
program in let's say five time steps and see what it says.
We will break up teh time interval into five equal
sub integrals and then use the method to
approximate what X sub five should be. And we get a
value of about 2.4. Now, that's with 5
time steps. If we use a greater number of time
steps, Euler's Method will return, hopefully,
a more accurate approximation to that final value.
And as we move towards, let's say, 80 times steps, we seem to be converging.
The final answer looks like it ought to be about 2.25, 2.24, something like that.
It's hard to say without running more examples.
Now why is it that Euler's
method works and seems to converge? How quickly does
it converge? Well, we're going to think from a Taylor's
series point of view, consider what happens when we
expand the solution X of T about T naught. X of
T naught of H is by, a Taylor series,
X naught plus H. Times the derivative of x with respect
to T, evaluated at T not, plus something in big O of H squared.
Now what's that derivative evaluated at T naught?
That is simply F at X naught, T naught.
And here, we see Euler's method hiding inside of
the Taylor Series. Euler's method is really just Taylor
expansion and dropping all of the terms of order to and greater.
Now, we say, because of this, that Euler's method is a first order method.
And that the error you make at each step is in big O of
h squared. Where h is our time step.
Remember.
Now, what is the net error when you do capital N such steps?
Well, it's capital n times something in big O of h squared.
You might think, ha, ha, capital n is a constant,
so I can forget about it, since we're doing big O.
No, capital N depends
on age. It's really some constant, the change in
time, divided by H. And so being careful,
we see that the net error is in big O of H
We'll look at two other methods for doing numerical ODEs.
The first is called the midpoint method. It begins easily enough
just like the Euler method, one projects forward.
Call this change in x that we get from Euler's method kappa.
Now, this is the midpoint method. So what we're going to do is take the time
value that is in between t sub n and t n plus one,
and try to approximate what the differential equation is doing there.
We will use evaluation at that midpoint to obtain a
new slope and pull that back to X of n, and then project
that forward at that midpoint slope to get x and plus 1.
Now, when we write all this out algebraically, it involves
this constant kappa and it looks a little bit complicated.
You don't have to memorize this, but what you do need to know is that this
is a second-order method. That means that
we're really doing a Taylor expansion, and keeping all the terms up
to, and including degree to the step error in the
midpoint method is in big O of h cubed, meaning that the net
error is in big O of h squared. This is a better
approximation at the expense of being algebraically more complex.
Now with that in mind, you might wonder, wouldn't it be possible
to do higher order approximations using midpoints, pulling it
back, and then evaluating again, pushing it forward, and pulling it back and
averaging everything all up together. Indeed, this is possible, and this
is a wonderful method called the Runge-Kutta method.
There's a lot of equations that go into this.
You don't have to memorize it. What you need to know is that first, this
exists And second, it's very helpful even if it seems
rather mysterious. This is a 4th order method.
Very good approximation.
Your step error is in big O of h to the 5th,
and your net error is in big O of h to the 4th.
Let's compare all these methods together in the
one simple example that we know how to solve.
Dx dt equals x with the initial time 0,
initial value 1 of final time 1, and our final
value, x star, which is known to be. We've
already seen what happens in the general case in Euler,
as H is going to 0. Let's be dumb, let's use
1 step. We're going to take our step size H to be
equal to 1 In Euler's method,
what we get, well, we get, simply that
X1 equals 1 plus H times X not, since Xnot
is equal to 1, we get an approximation of 1 plus 1.
That's not a terrible approximation to e, given that we've
only done one step, but that's all we get from Euler.
Now, with the midpoint method, our kappa value.
We already know what that is.
That is equal to 1.
What happens when we take the formula for the midpoint method?
Well, we get a value of X1 equal to 1 plus
h times quantity x nought plus one-half
kappa. That means we get 1 plus 1 plus one-half.
That's a better approximation to A. Finally when we go all the
way and write down all of those formulae for the Runge-Kutta Fourth Order Method.
Then with a little bit of algebra, you'll see we get 1 plus 1
plus a half plus 1 6th plus 1 24th. And now
you see the pattern. These numerical methods are really
just giving higher and higher terms from the Taylor series.
Moral of our story is that, by looking
at things from a Taylor series point of view, you can understand errors.
We've now seen just a small taste of how
discrete or digital methods help us solve smooth calculus problems.
But this is, by no means, the only possible example.
In our next lesson, we'll consider how to
use the discreet calculus to solve problems involving
definite integrals.