WEBVTT

00:00:04.080 --> 00:00:06.280
PROFESSOR: Software
that implements

00:00:06.280 --> 00:00:09.920
modern numerical methods
has two features that

00:00:09.920 --> 00:00:15.940
aren't present in codes like
ODE4 and classical Runge-Kutta.

00:00:15.940 --> 00:00:21.860
The methods in the software
can estimate error and provide

00:00:21.860 --> 00:00:24.930
automatic step size control.

00:00:24.930 --> 00:00:28.230
You don't specify
the step size h.

00:00:28.230 --> 00:00:32.530
You specify an
accuracy you want.

00:00:32.530 --> 00:00:36.600
And the methods
estimate the errors

00:00:36.600 --> 00:00:41.460
as they go along and adjust
the step size accordingly.

00:00:41.460 --> 00:00:43.400
And they provide
a fully accurate

00:00:43.400 --> 00:00:45.210
continuous interpolant.

00:00:45.210 --> 00:00:48.150
They don't just
provide the solution

00:00:48.150 --> 00:00:51.160
at the discrete set of points.

00:00:51.160 --> 00:00:56.950
They provide a function that
defines the solution everywhere

00:00:56.950 --> 00:00:58.740
in the interval.

00:00:58.740 --> 00:01:04.890
And so you can plot it,
find zeroes of the function,

00:01:04.890 --> 00:01:09.300
provide a facility called
event handling, and so on.

00:01:09.300 --> 00:01:12.680
Larry Shampine is an authority
on the numerical solution

00:01:12.680 --> 00:01:14.840
of ordinary
differential equations.

00:01:14.840 --> 00:01:19.250
He is the principal
author of this textbook

00:01:19.250 --> 00:01:22.490
about solving ODEs with MATLAB.

00:01:22.490 --> 00:01:28.180
He's a, now, emeritus professor
at the Southern Methodist

00:01:28.180 --> 00:01:29.940
University in Dallas.

00:01:29.940 --> 00:01:33.550
And he's been a long time
consultant to the MathWorks

00:01:33.550 --> 00:01:38.550
about the development
of our ODE Suite.

00:01:38.550 --> 00:01:42.330
Shampine and his student,
Przemyslaw Bogacki,

00:01:42.330 --> 00:01:45.890
published this method in 1989.

00:01:45.890 --> 00:01:51.110
And it's the basis
for ODE23, the first

00:01:51.110 --> 00:01:57.100
of the methods we will use
out of the MATLAB ODE Suite.

00:01:57.100 --> 00:01:59.970
The basic method is order three.

00:01:59.970 --> 00:02:02.780
And the error estimate is
based on the difference

00:02:02.780 --> 00:02:06.710
between the order three method
and then the underlying order

00:02:06.710 --> 00:02:08.060
two method.

00:02:08.060 --> 00:02:12.170
There are four slopes involved.

00:02:12.170 --> 00:02:14.620
The first one is the
value of the function

00:02:14.620 --> 00:02:16.890
at the start of the interval.

00:02:16.890 --> 00:02:19.940
But that's based on
something called FSAL,

00:02:19.940 --> 00:02:24.930
first same as last, where
that slope is most likely

00:02:24.930 --> 00:02:27.750
left over from
the previous step.

00:02:27.750 --> 00:02:30.510
If the previous
step was successful,

00:02:30.510 --> 00:02:35.590
this function value is the
same as the last function

00:02:35.590 --> 00:02:38.530
value from the previous step.

00:02:38.530 --> 00:02:42.640
That slope is used to step into
the middle of the interval,

00:02:42.640 --> 00:02:45.330
function is evaluated there.

00:02:45.330 --> 00:02:48.480
That slope is used to
step 3/4 of the way

00:02:48.480 --> 00:02:53.730
across the interval and a
third slope obtained there.

00:02:53.730 --> 00:02:58.420
Then these three values
are used to take the step.

00:02:58.420 --> 00:03:01.960
yn plus 1 is a
linear combination

00:03:01.960 --> 00:03:04.980
of these three function values.

00:03:04.980 --> 00:03:09.850
Then the function is evaluated
to get a fourth slope

00:03:09.850 --> 00:03:11.820
at the end of the interval.

00:03:11.820 --> 00:03:17.750
And then, these four slopes
are used to estimate the error.

00:03:17.750 --> 00:03:24.550
The error estimate here is the
difference between yn plus 1

00:03:24.550 --> 00:03:29.400
and another estimate
of the solution

00:03:29.400 --> 00:03:32.680
that's obtained from
a second order method

00:03:32.680 --> 00:03:36.000
that we don't actually evaluate.

00:03:36.000 --> 00:03:39.730
We just need the difference
between that method and yn

00:03:39.730 --> 00:03:42.560
plus 1 to estimate the error.

00:03:42.560 --> 00:03:49.070
This estimated error is compared
with a user-supplied tolerance.

00:03:49.070 --> 00:03:53.170
If the estimated error
is less than a tolerance,

00:03:53.170 --> 00:03:55.370
then the step is successful.

00:03:55.370 --> 00:04:00.000
And this fourth
slope, s4, becomes

00:04:00.000 --> 00:04:03.010
the s1 of the next step.

00:04:03.010 --> 00:04:06.120
If the answer is bigger
than the tolerance,

00:04:06.120 --> 00:04:08.600
then the error
could be the basis

00:04:08.600 --> 00:04:10.800
for adjusting the step size.

00:04:10.800 --> 00:04:13.330
In either case,
the error estimate

00:04:13.330 --> 00:04:18.410
is the basis for adjusting the
step size for the next step.

00:04:18.410 --> 00:04:23.830
This is the Bogacki-Shampine
Order 3(2) Method

00:04:23.830 --> 00:04:27.970
that's the basis for ODE 23.

00:04:34.060 --> 00:04:41.030
Let's look at some very simple
uses of ODE23 just to get

00:04:41.030 --> 00:04:42.240
started.

00:04:42.240 --> 00:04:44.690
I'm going to take the
differential equation

00:04:44.690 --> 00:04:46.950
y prime is equal to y.

00:04:46.950 --> 00:04:51.500
So I'm going to
compute e to the t.

00:04:51.500 --> 00:04:58.720
And just call ODE23 on
the interval from 0 to 1,

00:04:58.720 --> 00:05:01.320
with initial value 1.

00:05:01.320 --> 00:05:03.180
No output arguments.

00:05:03.180 --> 00:05:08.210
If I call it ODE23, it
just plots the solution.

00:05:08.210 --> 00:05:09.000
Here it is.

00:05:09.000 --> 00:05:10.850
It just produces a plot.

00:05:10.850 --> 00:05:15.620
It picks a step size,
goes from 0 to 1,

00:05:15.620 --> 00:05:24.440
and here it gets the final
value of e-- 2.7 something.

00:05:24.440 --> 00:05:28.330
If I do supply output arguments.

00:05:28.330 --> 00:05:34.880
I say t comma y equals
ODE23, it comes back

00:05:34.880 --> 00:05:37.780
with values of t and y.

00:05:37.780 --> 00:05:41.810
ODE23 picks the
values of t it wants.

00:05:41.810 --> 00:05:43.810
This is a trivial problem.

00:05:43.810 --> 00:05:47.900
It ends up picking
a step size of 0.1.

00:05:47.900 --> 00:05:50.700
After it gets started, it
chooses an initial step size

00:05:50.700 --> 00:05:57.580
of .08 for whatever
error tolerances.

00:05:57.580 --> 00:06:05.470
And the final value of y is
2.718, which is the value of e.

00:06:05.470 --> 00:06:10.190
So these are the two
simple uses of ODE23.

00:06:10.190 --> 00:06:15.130
If you don't supply any output
arguments, it draws a graph.

00:06:15.130 --> 00:06:18.460
If you do supply output
arguments, t and y,

00:06:18.460 --> 00:06:22.610
it comes back with
the values of t and y

00:06:22.610 --> 00:06:27.160
choosing the values of
t to meet the error.

00:06:27.160 --> 00:06:31.090
The default error tolerances
is 10 to the minus 3.

00:06:31.090 --> 00:06:35.880
So this value is going to
be accurate to three digits.

00:06:35.880 --> 00:06:38.000
And sure enough
that's what we got.

00:06:42.379 --> 00:06:43.920
Now let's try
something a little more

00:06:43.920 --> 00:06:49.050
challenging to see the automatic
error-controlled step size

00:06:49.050 --> 00:06:51.690
choice in action.

00:06:51.690 --> 00:06:54.870
Set a equal to a quarter.

00:06:54.870 --> 00:07:00.940
And then set y0 equal to 15.9.

00:07:00.940 --> 00:07:05.440
If I would set it to 16,
which is 1 over a squared,

00:07:05.440 --> 00:07:08.190
I'd run into a singularity.

00:07:08.190 --> 00:07:11.240
Now the differential
equation is y

00:07:11.240 --> 00:07:18.580
prime is equal to 2 (a
minus t) times y squared.

00:07:18.580 --> 00:07:25.140
I'm going to integrate this
with the ODE23 on the interval

00:07:25.140 --> 00:07:33.010
from 0 to 1 starting at y0, and
saving the results in t and y,

00:07:33.010 --> 00:07:36.530
and then plotting them.

00:07:36.530 --> 00:07:42.810
So here's my plot command,
and there is the solution.

00:07:42.810 --> 00:07:47.770
So there is a near
singularity at a.

00:07:47.770 --> 00:07:49.490
It nearly blows up.

00:07:49.490 --> 00:07:51.730
And then it settles back down.

00:07:51.730 --> 00:07:56.270
So the points are
bunched together

00:07:56.270 --> 00:07:59.800
as you go up to the
singularity and come back down,

00:07:59.800 --> 00:08:04.340
but then get farther apart
as the solution settles down.

00:08:04.340 --> 00:08:11.050
And the ODE solver is
able to take bigger steps.

00:08:11.050 --> 00:08:14.060
To see what steps
were actually taken,

00:08:14.060 --> 00:08:19.985
let's compute the difference
of t, and then plot that.

00:08:25.620 --> 00:08:30.100
So here are the step
sizes that were taken.

00:08:30.100 --> 00:08:34.580
And we see that
a small step size

00:08:34.580 --> 00:08:40.980
was taken near the almost
singularity at that 0.25.

00:08:40.980 --> 00:08:44.380
And then as we get towards
the end of the interval,

00:08:44.380 --> 00:08:46.810
a larger step size is taken.

00:08:46.810 --> 00:08:49.820
And then, finally,
the step size just

00:08:49.820 --> 00:08:55.920
to reach the end of the interval
is taken as the last step.

00:08:55.920 --> 00:09:00.920
So that's the automatic
step size choice of ODE23.

00:09:06.230 --> 00:09:11.860
BS23 has a nice
natural interpolant

00:09:11.860 --> 00:09:14.590
that goes along with
it that's actually

00:09:14.590 --> 00:09:16.700
been known for over 100 years.

00:09:16.700 --> 00:09:21.510
It's called Hermite
Cubic Interpolation.

00:09:21.510 --> 00:09:25.850
We know that two points
determine a straight line.

00:09:25.850 --> 00:09:30.520
Well, two points and two
slopes determine a cubic.

00:09:30.520 --> 00:09:35.960
On each interval we have the
values of y and yn plus 1.

00:09:35.960 --> 00:09:38.710
We also have two
slopes, namely this.

00:09:38.710 --> 00:09:44.050
We have the derivatives at the
end points, yn prime and yn

00:09:44.050 --> 00:09:49.450
plus 1 prime, that's the
values of the differential

00:09:49.450 --> 00:09:51.490
equation at those points.

00:09:51.490 --> 00:09:55.000
So those four values
determine a cubic

00:09:55.000 --> 00:09:57.120
that goes through
those two points

00:09:57.120 --> 00:09:59.830
and has those two slopes.

00:09:59.830 --> 00:10:06.100
This cubic allows the software
to evaluate the solution

00:10:06.100 --> 00:10:10.390
at any point in the interval
without additional cost

00:10:10.390 --> 00:10:15.380
as defined by addition
evaluations of the function f.

00:10:15.380 --> 00:10:19.230
This can be used to draw
graphs of the solution,

00:10:19.230 --> 00:10:22.020
nice smooth graphs
of the solution,

00:10:22.020 --> 00:10:27.830
find zeroes of the solution,
do event handling, and so on.

00:10:27.830 --> 00:10:31.360
Another feature
provided by ODE23.