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.