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.