WEBVTT

00:00:04.470 --> 00:00:06.570
PROFESSOR: The most
frequently used ODE solver

00:00:06.570 --> 00:00:09.865
in MATLAB and Simulink is ODE45.

00:00:09.865 --> 00:00:14.700
It is based on method published
by British mathematicians JR

00:00:14.700 --> 00:00:19.230
Dormand and PJ Prince in 1980.

00:00:19.230 --> 00:00:22.760
The basic method is order five.

00:00:22.760 --> 00:00:29.350
The error correction uses a
companion order four method.

00:00:29.350 --> 00:00:33.610
The slope of tn is, first
same as last left over

00:00:33.610 --> 00:00:37.020
from the previous
successful step.

00:00:37.020 --> 00:00:44.070
Then there are five more slopes
from function values at 1/5 h,

00:00:44.070 --> 00:00:52.240
3/10h, 4/5h, 8/9h and
then at tn plus 1.

00:00:52.240 --> 00:00:55.710
These six slopes, linear
combinations of them,

00:00:55.710 --> 00:00:59.540
are used to produce yn plus 1.

00:00:59.540 --> 00:01:05.900
The function is evaluated
at tn plus 1 and yn plus 1

00:01:05.900 --> 00:01:07.950
to get a seventh slope.

00:01:07.950 --> 00:01:10.940
And then linear
combinations of these

00:01:10.940 --> 00:01:14.760
are used to produce
the error estimate.

00:01:14.760 --> 00:01:17.630
Again, if the error
estimate is less

00:01:17.630 --> 00:01:22.610
than the specified
accuracy requirements

00:01:22.610 --> 00:01:24.800
the step is successful.

00:01:24.800 --> 00:01:31.750
And then that error estimate is
used to get the next step size.

00:01:31.750 --> 00:01:36.460
If the error is too big,
the step is unsuccessful

00:01:36.460 --> 00:01:42.390
and that error estimate
is used to get the step

00:01:42.390 --> 00:01:44.800
size to do the step over again.

00:01:44.800 --> 00:01:47.890
If we want to see the actual
coefficients that are used,

00:01:47.890 --> 00:01:51.550
you can go into
the code for ODE45.

00:01:51.550 --> 00:01:53.850
There's a table with
the coefficients.

00:01:53.850 --> 00:01:59.190
Or you go to the Wikipedia page
for the Dormand-Prince Method

00:01:59.190 --> 00:02:01.160
and there is the
same coefficients.

00:02:04.590 --> 00:02:08.150
As an aside, here is an
interesting fact about higher

00:02:08.150 --> 00:02:10.780
order Runge-Kutta methods.

00:02:10.780 --> 00:02:15.010
Classical Runge-Kutta required
four function evaluations

00:02:15.010 --> 00:02:18.800
per step to get order four.

00:02:18.800 --> 00:02:24.490
Dormand-Prince requires six
function evaluations per step

00:02:24.490 --> 00:02:27.060
to get order five.

00:02:27.060 --> 00:02:32.760
You can't get order five with
just five function evaluations.

00:02:32.760 --> 00:02:36.560
And then, if we were to try
and achieve higher order,

00:02:36.560 --> 00:02:41.690
it would take even more
function evaluations per step.

00:02:41.690 --> 00:02:47.136
Let's use ODE45 to
compute e to the t.

00:02:47.136 --> 00:02:50.910
y prime is equal to y.

00:02:50.910 --> 00:02:54.660
We can ask for
output by supplying

00:02:54.660 --> 00:02:58.700
an argument called tspan.

00:02:58.700 --> 00:03:02.710
Zero and steps of 0.1 to 1.

00:03:02.710 --> 00:03:05.750
If we supply that as
the input argument

00:03:05.750 --> 00:03:07.700
to solve this
differential equation

00:03:07.700 --> 00:03:11.700
and get the output
at those points,

00:03:11.700 --> 00:03:14.220
we get that back as the output.

00:03:14.220 --> 00:03:17.860
And now here's the
approximations to the solution

00:03:17.860 --> 00:03:21.370
to that differential
equation at those points.

00:03:21.370 --> 00:03:30.500
If we plot it, here's the
solution at those points.

00:03:30.500 --> 00:03:33.700
And to see how
accurate it is, we

00:03:33.700 --> 00:03:42.140
see that we're actually getting
this result to nine digits.

00:03:42.140 --> 00:03:45.435
ODE45 is very accurate.

00:03:49.740 --> 00:03:52.520
Let's look at step size
choice on our problem

00:03:52.520 --> 00:03:57.300
with near singularity,
is a quarter.

00:03:57.300 --> 00:04:01.270
y0 is close to 16.

00:04:01.270 --> 00:04:08.390
The differential equation is
y prime is 2(a-t) y squared.

00:04:08.390 --> 00:04:15.230
We let ODE45 choose its own
step size by indicating we just

00:04:15.230 --> 00:04:19.970
want to integrate from 0 to 1.

00:04:19.970 --> 00:04:26.000
We capture the output
in t and y and plot it.

00:04:29.870 --> 00:04:34.360
Now, here, there's a
lot of points here,

00:04:34.360 --> 00:04:41.310
but this is misleading
because ODE45, by default,

00:04:41.310 --> 00:04:44.150
is using the refine option.

00:04:44.150 --> 00:04:47.530
It's only actually
evaluating the function

00:04:47.530 --> 00:04:50.440
at every fourth
one of these points

00:04:50.440 --> 00:04:56.880
and then using the interpolant
to fill in in between.

00:04:56.880 --> 00:05:03.380
So we actually need a
little different plot here.

00:05:06.380 --> 00:05:09.240
This plot shows a little
better what's going on.

00:05:09.240 --> 00:05:13.820
The big dots are the
points that ODE45

00:05:13.820 --> 00:05:18.180
chose to evaluate the
differential equation.

00:05:18.180 --> 00:05:23.050
And the little dots are filled
in with the interpolant.

00:05:23.050 --> 00:05:26.570
So the big dots are
every fourth point.

00:05:26.570 --> 00:05:33.740
And the refine option says that
the big dots are too far apart

00:05:33.740 --> 00:05:36.100
and we need to fill it
in with the interpolant.

00:05:36.100 --> 00:05:41.060
And so this is the continuous
interpolant in action.

00:05:41.060 --> 00:05:46.600
The big dots are more
closely concentrated

00:05:46.600 --> 00:05:49.010
as we have to go
around the curve.

00:05:49.010 --> 00:05:55.430
And then, as we get farther
away from the singularity

00:05:55.430 --> 00:05:57.970
the step size increases.

00:05:57.970 --> 00:06:04.280
So this shows the
high accuracy of ODE45

00:06:04.280 --> 00:06:07.745
and the automatic step
size choice in action.

00:06:11.750 --> 00:06:13.870
Here's an exercise.

00:06:13.870 --> 00:06:21.700
Compare ODE23 and ODE45 by using
each of them to compute pi.

00:06:21.700 --> 00:06:28.530
The integral 4 over 1 plus
t squared from 0 to 1 is pi.

00:06:28.530 --> 00:06:32.700
You can express that as
a differential equation,

00:06:32.700 --> 00:06:36.040
use each of the routines to
integrate that differential

00:06:36.040 --> 00:06:40.900
equation and see how close
they get to computing pi.