WEBVTT

00:00:03.304 --> 00:00:06.910
INSTRUCTOR: I want to
illustrate the important notion

00:00:06.910 --> 00:00:14.600
of stiffness by running ode45,
the primary MATLAB ODE solver,

00:00:14.600 --> 00:00:18.400
on our flame example.

00:00:18.400 --> 00:00:21.580
The differential
equation is y prime is y

00:00:21.580 --> 00:00:25.070
squared minus y
cubed, and I'm going

00:00:25.070 --> 00:00:29.770
to choose a fairly--
an extremely small

00:00:29.770 --> 00:00:33.620
initial condition, 10
to the minus sixth.

00:00:33.620 --> 00:00:37.500
The final value of t
is 2 over y naught,

00:00:37.500 --> 00:00:45.170
and I'm going to impose a
modest accuracy requirement, 10

00:00:45.170 --> 00:00:46.970
to the minus fifth.

00:00:46.970 --> 00:00:55.440
Now let's run ode45
with its default output.

00:00:55.440 --> 00:00:59.420
Now, see it's taking-- it's
moving very slowly here.

00:00:59.420 --> 00:01:01.950
It's taking lots of steps.

00:01:01.950 --> 00:01:05.019
So I'm take- pressing
the stop button here.

00:01:05.019 --> 00:01:06.830
It's working very hard.

00:01:06.830 --> 00:01:17.550
Let's zoom in and see why it's
taking so many steps, very

00:01:17.550 --> 00:01:20.480
densely packed steps here.

00:01:27.690 --> 00:01:30.960
This is stiffness.

00:01:30.960 --> 00:01:36.870
It's satisfying the accuracy
requirements we imposed.

00:01:36.870 --> 00:01:42.310
All these steps are within
10 to the minus sixth of one,

00:01:42.310 --> 00:01:45.080
but it's taken very
small steps to do it.

00:01:45.080 --> 00:01:51.690
These steps are so small
that the graphics can't even

00:01:51.690 --> 00:01:54.210
discern the step size.

00:01:54.210 --> 00:01:56.400
This is stiffness.

00:01:56.400 --> 00:01:58.460
It's an efficiency issue.

00:01:58.460 --> 00:02:00.950
It's doing what we asked for.

00:02:00.950 --> 00:02:04.480
It's meeting the
accuracy requirements,

00:02:04.480 --> 00:02:08.770
but it's having to take
very small steps to do it.

00:02:11.420 --> 00:02:13.460
Let's try another
ODE solver-- ode23.

00:02:17.180 --> 00:02:22.380
Just change this to 23
and see what it does.

00:02:22.380 --> 00:02:28.300
It's also taking very small
steps for the same reason.

00:02:28.300 --> 00:02:36.420
If we zoom in on here, we'll
see the same kind of behavior.

00:02:39.260 --> 00:02:41.400
But it's taking very
small steps in order

00:02:41.400 --> 00:02:45.900
to achieve the desired accuracy.

00:02:45.900 --> 00:02:50.460
Now let me introduce
a new solver, ode23s.

00:02:50.460 --> 00:02:53.240
The s for stiffness.

00:02:53.240 --> 00:02:56.700
This was designed to
solve stiff problems.

00:02:56.700 --> 00:03:00.570
And boom, it goes
up, turns the corner,

00:03:00.570 --> 00:03:08.750
and it takes just a few steps
to get to the final result.

00:03:08.750 --> 00:03:10.606
There it turns the
corner very quickly.

00:03:14.890 --> 00:03:18.140
We'll see how ode23s
works in a minute,

00:03:18.140 --> 00:03:22.050
but first let's try
to define stiffness.

00:03:22.050 --> 00:03:24.360
It's a qualitative
notion that doesn't

00:03:24.360 --> 00:03:28.040
have a precise
mathematical definition.

00:03:28.040 --> 00:03:30.540
It depends upon the
problem, but also

00:03:30.540 --> 00:03:36.020
on the solver and the
accuracy requirements.

00:03:36.020 --> 00:03:39.410
But it's an important notion.

00:03:39.410 --> 00:03:43.860
We say that a problem is stiff
if the solution being sought

00:03:43.860 --> 00:03:48.600
very slowly, but there are
nearby solutions that very

00:03:48.600 --> 00:03:49.910
rapidly.

00:03:49.910 --> 00:03:54.250
So the numerical method
must take small steps

00:03:54.250 --> 00:03:56.655
to obtain satisfactory results.

00:04:02.960 --> 00:04:05.710
Stiff methods for ordinary
differential equations

00:04:05.710 --> 00:04:08.020
must be implicit.

00:04:08.020 --> 00:04:12.260
They must involve formulas
that involve looking backward

00:04:12.260 --> 00:04:14.960
from the forward timestep.

00:04:14.960 --> 00:04:18.779
The prototype of these methods
is the backward Euler method,

00:04:18.779 --> 00:04:21.930
or the implicit Euler method.

00:04:21.930 --> 00:04:27.390
This formula, it involves--
defines y n plus 1,

00:04:27.390 --> 00:04:30.760
but doesn't tell us
how to compute it.

00:04:30.760 --> 00:04:34.946
We have to solve this
equation for y n plus 1.

00:04:34.946 --> 00:04:37.980
And I'm not going to go into
detail about how we actually

00:04:37.980 --> 00:04:38.480
do it.

00:04:38.480 --> 00:04:41.435
It involves something
like a Newton method

00:04:41.435 --> 00:04:44.790
that would-- requires
knowing the derivative,

00:04:44.790 --> 00:04:47.370
or an approximation to
the derivative of f.

00:04:50.130 --> 00:04:51.840
But this gives you
an idea of what you

00:04:51.840 --> 00:04:54.155
can expect in stiff methods.

00:04:57.270 --> 00:04:59.830
I like to make an
analogy with taking

00:04:59.830 --> 00:05:01.700
a hike in one of
the slot canyons

00:05:01.700 --> 00:05:04.040
we have here in the Southwest.

00:05:04.040 --> 00:05:08.330
Explicit methods
like ode23 and 45

00:05:08.330 --> 00:05:11.100
take steps on the
walls of the canyon

00:05:11.100 --> 00:05:16.300
and go back and forth across
the sides of the canyon,

00:05:16.300 --> 00:05:19.330
make very slow progress
down the canyon.

00:05:19.330 --> 00:05:25.080
Whereas implicit
methods, like ode15s,

00:05:25.080 --> 00:05:30.060
look ahead down the
canyon and look ahead

00:05:30.060 --> 00:05:34.255
to where you want to go and make
rapid progress of the canyon.

00:05:39.560 --> 00:05:48.090
The stiff solver, ode23s,
uses an implicit second-order

00:05:48.090 --> 00:05:55.220
formula and an associated
third-order error estimator.

00:05:55.220 --> 00:05:58.460
It evaluates the partial
derivatives of f with respect

00:05:58.460 --> 00:06:03.470
to both t and f at each
step, so that's expensive.

00:06:03.470 --> 00:06:07.610
It's efficient at
crude error tolerances,

00:06:07.610 --> 00:06:11.540
like graphic accuracy.

00:06:11.540 --> 00:06:15.920
And it has relatively
low overhead.

00:06:19.990 --> 00:06:26.220
By way of comparison,
the stiff solver ode15s,

00:06:26.220 --> 00:06:30.070
can be configured to use
either the variable order

00:06:30.070 --> 00:06:32.870
numerical
differentiation formula,

00:06:32.870 --> 00:06:41.820
NDF, or the related to backward
differentiation formula BDF.

00:06:41.820 --> 00:06:45.280
Neither case it
saves several values

00:06:45.280 --> 00:06:49.020
of the function
over previous steps.

00:06:49.020 --> 00:06:54.030
The order varies automatically
between one and five,

00:06:54.030 --> 00:06:59.530
it evaluates the partial
derivatives less frequently,

00:06:59.530 --> 00:07:06.140
and did see efficient at
higher tolerances then 23s.