WEBVTT
00:00:03.980 --> 00:00:05.740
CLEVE MOLER: Hello.
00:00:05.740 --> 00:00:08.970
I'm Cleve Moler,
one of the founders
00:00:08.970 --> 00:00:12.800
and chief mathematician
at the MathWorks.
00:00:12.800 --> 00:00:16.100
This series of videos is about
solving ordinary differential
00:00:16.100 --> 00:00:19.530
equations in MATLAB.
00:00:19.530 --> 00:00:23.440
We can begin by recalling
the definition of derivative.
00:00:23.440 --> 00:00:26.010
The derivative of a
function at a point
00:00:26.010 --> 00:00:31.620
is the slope of the
tangent line to the graph
00:00:31.620 --> 00:00:33.865
of the function at that point.
00:00:36.380 --> 00:00:39.150
Our numerical
approximations will
00:00:39.150 --> 00:00:44.910
rely upon the slope of
the secant to the graph.
00:00:44.910 --> 00:00:51.980
That's a line through two points
separated by a distance h.
00:00:51.980 --> 00:00:57.700
We'll have a lot to say about
the step size h as we go along.
00:00:57.700 --> 00:01:03.990
What's important to realize
is that as h goes to 0,
00:01:03.990 --> 00:01:06.950
the slope of the
secant approaches
00:01:06.950 --> 00:01:08.080
the slope of tangent.
00:01:11.810 --> 00:01:16.910
The wiggly equals sign means
approximately equal to.
00:01:16.910 --> 00:01:22.030
T0 is the point where we are
finding the approximation.
00:01:22.030 --> 00:01:26.170
The value of the
derivative at t0
00:01:26.170 --> 00:01:31.710
is approximately equal to
the slope of the secant.
00:01:31.710 --> 00:01:34.010
The slope of the
secant is the change
00:01:34.010 --> 00:01:37.840
in the y value over the
change in the t value.
00:01:37.840 --> 00:01:40.100
The change in y value
is the difference
00:01:40.100 --> 00:01:44.470
between the two values of y.
00:01:44.470 --> 00:01:48.690
The change in the t
value is the step size h.
00:01:48.690 --> 00:01:51.820
If we rewrite this,
we get the value of y
00:01:51.820 --> 00:01:56.200
at the point t0 plus
h is approximately
00:01:56.200 --> 00:02:01.020
equal to the value of
y at t0 plus h times
00:02:01.020 --> 00:02:05.030
the value of y prime at t0.
00:02:05.030 --> 00:02:09.270
This is the basis for our
first numerical method,
00:02:09.270 --> 00:02:09.959
Euler's method.
00:02:13.960 --> 00:02:19.670
Leonhard Euler was a 18th
century Swiss mathematician.
00:02:19.670 --> 00:02:25.500
Probably the most influential
mathematician of his era.
00:02:25.500 --> 00:02:29.130
He made important
contributions to a wide range
00:02:29.130 --> 00:02:34.710
of fields of mathematics,
physics, and astronomy.
00:02:34.710 --> 00:02:37.740
He invented the notion of
a function, for example.
00:02:41.190 --> 00:02:44.700
The differential equation
is given by this function
00:02:44.700 --> 00:02:48.300
f of two variables t and y.
00:02:48.300 --> 00:02:52.960
And the task in general
is to find the function y
00:02:52.960 --> 00:02:56.070
whose derivative is equal to f.
00:02:56.070 --> 00:03:00.570
Now, there's lots of functions y
whose derivative is equal to f.
00:03:00.570 --> 00:03:06.725
And so there's an initial
condition, a point t
00:03:06.725 --> 00:03:11.660
naught or t0, and a value y0.
00:03:11.660 --> 00:03:17.160
And the initial
condition is that y at t0
00:03:17.160 --> 00:03:19.580
should be equal to y0.
00:03:19.580 --> 00:03:22.220
Here's some examples.
00:03:22.220 --> 00:03:31.380
The compound interest problem is
just the interest rate times y.
00:03:31.380 --> 00:03:35.800
Here, the function of t and y
doesn't actually depend upon t.
00:03:35.800 --> 00:03:38.130
And it's linear in y.
00:03:38.130 --> 00:03:42.380
The initial condition
is at time 0,
00:03:42.380 --> 00:03:46.940
there's a specified
value of y, like $100.
00:03:46.940 --> 00:03:51.910
That's the compound
interest problem.
00:03:51.910 --> 00:03:54.290
Here's the logistic equation.
00:03:54.290 --> 00:03:55.940
Nonlinear equation.
00:03:55.940 --> 00:04:00.550
Here, f of t and y, again,
doesn't depend upon t.
00:04:00.550 --> 00:04:04.870
And it's a constant times y
minus another constant times
00:04:04.870 --> 00:04:06.520
y squared.
00:04:06.520 --> 00:04:08.770
That's the logistic equation.
00:04:08.770 --> 00:04:14.250
And again, the value
is specified at 0.
00:04:14.250 --> 00:04:19.750
Let's say y at 0 is equal to 1.
00:04:19.750 --> 00:04:22.040
Here's another
nonlinear equation.
00:04:22.040 --> 00:04:26.290
F of t and y is t
squared plus y squared.
00:04:26.290 --> 00:04:29.010
It's not possible to
find an analytic solution
00:04:29.010 --> 00:04:30.000
to this equation.
00:04:30.000 --> 00:04:33.170
We'll use these numerical
methods to find a solution
00:04:33.170 --> 00:04:34.760
to this equation.
00:04:34.760 --> 00:04:38.390
Initial condition y
at 0 is equal to 0.
00:04:38.390 --> 00:04:41.040
That's an example of
a function of t and y.
00:04:44.930 --> 00:04:48.490
Euler's method actually isn't
a practical numerical method
00:04:48.490 --> 00:04:50.330
in general.
00:04:50.330 --> 00:04:52.430
We're just using it
to get us started
00:04:52.430 --> 00:04:56.810
thinking about the ideas
underlying numerical methods.
00:04:56.810 --> 00:05:00.950
Euler's method involves
a sequence of points
00:05:00.950 --> 00:05:05.550
t sub n separated by
a fixed step size h.
00:05:05.550 --> 00:05:11.250
And then y sub n is an
approximation to the value
00:05:11.250 --> 00:05:13.740
of the solution at t sub n.
00:05:18.070 --> 00:05:22.950
The approximation comes from
the slope of the secant.
00:05:22.950 --> 00:05:30.710
The ratio of the difference
of the values of y
00:05:30.710 --> 00:05:35.060
and to the step size h.
00:05:35.060 --> 00:05:40.370
The differential equation
says that this ratio
00:05:40.370 --> 00:05:47.200
should be the value of
the function at t sub n.
00:05:47.200 --> 00:05:52.720
And if we rearrange
this equation,
00:05:52.720 --> 00:05:56.200
we get Euler's method.
00:05:56.200 --> 00:06:04.040
That yn plus 1 is yn plus
h times the function f
00:06:04.040 --> 00:06:07.770
evaluated at t
sub n and y sub n.
00:06:07.770 --> 00:06:08.850
This is Euler's method.
00:06:13.110 --> 00:06:18.250
We're now ready for our
first MATLAB program, ODE1.
00:06:18.250 --> 00:06:22.490
It's called ODE1 because
it's our first program
00:06:22.490 --> 00:06:25.340
and because it
evaluates the function f
00:06:25.340 --> 00:06:30.090
that defines the differential
equation once per step.
00:06:30.090 --> 00:06:33.760
There are five input arguments.
00:06:33.760 --> 00:06:38.640
The first is f, a function
that defines the differential
00:06:38.640 --> 00:06:39.850
equation.
00:06:39.850 --> 00:06:43.390
This is something called
an anonymous function.
00:06:43.390 --> 00:06:45.740
I'll talk more about
that in a moment.
00:06:45.740 --> 00:06:50.440
The other four are
scalar numerical values.
00:06:50.440 --> 00:06:55.340
The first three define the
interval of integration.
00:06:55.340 --> 00:07:01.140
We're going from t0 in
steps of h to t final.
00:07:01.140 --> 00:07:06.860
The fifth input argument
is y0, the initial value.
00:07:06.860 --> 00:07:09.600
The output is a vector.
00:07:09.600 --> 00:07:15.220
Vector y out is the
values of the solution
00:07:15.220 --> 00:07:20.210
at the points in the interval.
00:07:20.210 --> 00:07:25.130
We start by putting y0,
the initial value, into y
00:07:25.130 --> 00:07:29.290
and then putting y
into the output vector.
00:07:29.290 --> 00:07:33.120
The body of the
function is a four loop,
00:07:33.120 --> 00:07:41.610
t goes from T0 not in steps of H
up to one step short of t final
00:07:41.610 --> 00:07:46.710
and then the final passage
through the body of the code
00:07:46.710 --> 00:07:49.730
takes t up to t final.
00:07:49.730 --> 00:07:54.180
We evaluate the
function f at t and y.
00:07:54.180 --> 00:07:59.210
That gives us a slope
s, s is for slope.
00:07:59.210 --> 00:08:01.710
Here's the Euler step.
00:08:01.710 --> 00:08:07.300
Take the current value of
y, add h, times the slope.
00:08:07.300 --> 00:08:09.980
That gives us this
new value of y.
00:08:09.980 --> 00:08:13.150
And then y is appended to y out.
00:08:13.150 --> 00:08:19.190
This MATLAB construction
with the square brackets
00:08:19.190 --> 00:08:22.990
takes a vector y, adds
another value to it,
00:08:22.990 --> 00:08:27.970
making it one element longer
and puts the resulting y out
00:08:27.970 --> 00:08:29.890
back in y out.
00:08:29.890 --> 00:08:31.530
This is the entire code.
00:08:31.530 --> 00:08:32.600
This is it.
00:08:32.600 --> 00:08:36.619
This is ODE1 that
implements Euler's method.
00:08:39.720 --> 00:08:43.909
The first argument to any
of the MATLAB ODE solvers
00:08:43.909 --> 00:08:47.370
is the name of a function that
specifies the differential
00:08:47.370 --> 00:08:48.790
equation.
00:08:48.790 --> 00:08:52.070
This is known as
a function handle.
00:08:52.070 --> 00:08:55.710
The easiest way to
get a function handle
00:08:55.710 --> 00:09:00.400
is to make use of an
anonymous function created
00:09:00.400 --> 00:09:02.915
with the ampersand or at sign.
00:09:04.370 --> 00:09:07.060
All of the
differential equations
00:09:07.060 --> 00:09:13.940
involve anonymous functions
of two variables, t and y.
00:09:13.940 --> 00:09:20.070
And so we have f equals
at parenthesis t comma y
00:09:20.070 --> 00:09:22.060
closed parenthesis.
00:09:22.060 --> 00:09:27.590
This is followed by any
expression involving
00:09:27.590 --> 00:09:29.980
either t or y.
00:09:29.980 --> 00:09:33.440
And many of them
don't depend upon t.
00:09:33.440 --> 00:09:41.230
So here is an anonymous function
defining our interest problem.
00:09:41.230 --> 00:09:49.340
And we can just evaluate this
like any ordinary function.
00:09:49.340 --> 00:09:55.870
When y is equal to
1, f of 1 is 0.06.
00:09:55.870 --> 00:10:03.200
Here's an example of a function
that depends upon both t and y.
00:10:03.200 --> 00:10:10.670
The functions can involve
constants that have values.
00:10:10.670 --> 00:10:16.140
So here, we can
define two constants.
00:10:16.140 --> 00:10:18.940
And then we can use
those two constants
00:10:18.940 --> 00:10:25.270
to define the logistic
equation f of a times
00:10:25.270 --> 00:10:28.340
y minus b times y squared.
00:10:28.340 --> 00:10:33.190
Again, this is an autonomous
equation that doesn't actually
00:10:33.190 --> 00:10:34.990
depend upon t.
00:10:38.810 --> 00:10:42.370
Let's see how Euler's
method and ODE1
00:10:42.370 --> 00:10:44.560
work on this simple example.
00:10:44.560 --> 00:10:48.990
Y prime equals 2y with the
initial condition y of 0
00:10:48.990 --> 00:10:54.085
equals 10 on the interval
t between 0 and 3.
00:10:58.360 --> 00:11:06.310
We define the anonymous function
f of t and y is equal to 2y.
00:11:06.310 --> 00:11:10.190
The initial condition
is t0 equals 0.
00:11:10.190 --> 00:11:13.340
We're going to take
a step size of 1.
00:11:13.340 --> 00:11:18.570
Go to t final equals
3 starting at y0
00:11:18.570 --> 00:11:24.310
equals 10 and here's
our call to ODE1.
00:11:28.170 --> 00:11:32.740
We have an animation
that shows these steps.
00:11:32.740 --> 00:11:37.750
Start at t0 equals
0 and y0 equals 10.
00:11:37.750 --> 00:11:41.170
Here's our first point.
00:11:41.170 --> 00:11:44.800
We evaluate the function there.
00:11:44.800 --> 00:11:46.590
We get a slope of 20.
00:11:46.590 --> 00:11:49.110
That's 2 times 10.
00:11:49.110 --> 00:11:56.280
We take an Euler step of
length 1 across the first step.
00:11:56.280 --> 00:12:01.130
That brings us to
the second point, 30.
00:12:01.130 --> 00:12:03.660
Evaluate the function there.
00:12:03.660 --> 00:12:05.370
Two times 30 is 60.
00:12:05.370 --> 00:12:07.120
That's our slope.
00:12:07.120 --> 00:12:11.370
Take the second
step to get to y2.
00:12:11.370 --> 00:12:13.260
Y2 is 90.
00:12:13.260 --> 00:12:15.820
Evaluate the function there.
00:12:15.820 --> 00:12:19.050
Get 2 times 90 is 180.
00:12:19.050 --> 00:12:22.030
That gives us a slope.
00:12:22.030 --> 00:12:24.590
Take a step across the
interval with that slope
00:12:24.590 --> 00:12:26.900
would get us to a third point.
00:12:26.900 --> 00:12:29.710
The third point
is 270 and that's
00:12:29.710 --> 00:12:31.570
the end of the integration.
00:12:31.570 --> 00:12:38.420
So that's three Euler steps
to get from t0 to t final.
00:12:41.990 --> 00:12:43.800
Euler's method is
actually the same
00:12:43.800 --> 00:12:46.450
as computing compound interest.
00:12:46.450 --> 00:12:49.620
So let's do a compound
interest problem.
00:12:49.620 --> 00:12:52.770
Define the interest rate.
00:12:52.770 --> 00:12:57.960
Define our anonymous function
using that interest rate.
00:12:57.960 --> 00:13:01.020
Start at time 0.
00:13:01.020 --> 00:13:04.650
Go in steps of a month.
00:13:04.650 --> 00:13:08.010
Go for 10 years.
00:13:08.010 --> 00:13:13.990
Start with $100.
00:13:13.990 --> 00:13:19.000
And here's our
result of using ODE1
00:13:19.000 --> 00:13:21.580
to compute compound interest.
00:13:21.580 --> 00:13:24.970
That's 121 numbers.
00:13:24.970 --> 00:13:33.040
MATLAB actually has a format for
looking at dollars and cents.
00:13:33.040 --> 00:13:41.840
And so here they are as dollars
and cents starting with $100
00:13:41.840 --> 00:13:44.700
and compounding every month.
00:13:44.700 --> 00:13:52.860
We get up to just over $180.
00:13:52.860 --> 00:13:54.730
I'm going to plot that.
00:13:54.730 --> 00:13:59.950
So I want a time vector months.
00:13:59.950 --> 00:14:03.590
And I actually want to compare
it with simple interest.
00:14:03.590 --> 00:14:07.560
So here's how you
compute simple interest.
00:14:07.560 --> 00:14:09.100
$0.50 a month.
00:14:09.100 --> 00:14:11.570
And now let's plot those two.
00:14:19.910 --> 00:14:27.310
So the straight line is simple
interest getting up to $160.
00:14:27.310 --> 00:14:32.400
And the blue line is
the compound interest.
00:14:32.400 --> 00:14:35.650
There is a slight
upward curvature,
00:14:35.650 --> 00:14:39.270
getting us up to $180.
00:14:39.270 --> 00:14:42.090
There's a dot every
month here as we
00:14:42.090 --> 00:14:46.030
show the results of Euler's
method, which as I said
00:14:46.030 --> 00:14:49.845
is the same thing as
computing compound interest.
00:14:54.300 --> 00:14:57.260
Finally, here's an exercise.
00:14:57.260 --> 00:15:01.590
Find the differential equation
that produces linear growth
00:15:01.590 --> 00:15:07.360
and rerun this example
using ODE1 twice,
00:15:07.360 --> 00:15:10.290
once to compute the
compound interest
00:15:10.290 --> 00:15:13.510
and once to compute
the simple interest.