1
00:00:04,080 --> 00:00:06,280
PROFESSOR: Software
that implements

2
00:00:06,280 --> 00:00:09,920
modern numerical methods
has two features that

3
00:00:09,920 --> 00:00:15,940
aren't present in codes like
ODE4 and classical Runge-Kutta.

4
00:00:15,940 --> 00:00:21,860
The methods in the software
can estimate error and provide

5
00:00:21,860 --> 00:00:24,930
automatic step size control.

6
00:00:24,930 --> 00:00:28,230
You don't specify
the step size h.

7
00:00:28,230 --> 00:00:32,530
You specify an
accuracy you want.

8
00:00:32,530 --> 00:00:36,600
And the methods
estimate the errors

9
00:00:36,600 --> 00:00:41,460
as they go along and adjust
the step size accordingly.

10
00:00:41,460 --> 00:00:43,400
And they provide
a fully accurate

11
00:00:43,400 --> 00:00:45,210
continuous interpolant.

12
00:00:45,210 --> 00:00:48,150
They don't just
provide the solution

13
00:00:48,150 --> 00:00:51,160
at the discrete set of points.

14
00:00:51,160 --> 00:00:56,950
They provide a function that
defines the solution everywhere

15
00:00:56,950 --> 00:00:58,740
in the interval.

16
00:00:58,740 --> 00:01:04,890
And so you can plot it,
find zeroes of the function,

17
00:01:04,890 --> 00:01:09,300
provide a facility called
event handling, and so on.

18
00:01:09,300 --> 00:01:12,680
Larry Shampine is an authority
on the numerical solution

19
00:01:12,680 --> 00:01:14,840
of ordinary
differential equations.

20
00:01:14,840 --> 00:01:19,250
He is the principal
author of this textbook

21
00:01:19,250 --> 00:01:22,490
about solving ODEs with MATLAB.

22
00:01:22,490 --> 00:01:28,180
He's a, now, emeritus professor
at the Southern Methodist

23
00:01:28,180 --> 00:01:29,940
University in Dallas.

24
00:01:29,940 --> 00:01:33,550
And he's been a long time
consultant to the MathWorks

25
00:01:33,550 --> 00:01:38,550
about the development
of our ODE Suite.

26
00:01:38,550 --> 00:01:42,330
Shampine and his student,
Przemyslaw Bogacki,

27
00:01:42,330 --> 00:01:45,890
published this method in 1989.

28
00:01:45,890 --> 00:01:51,110
And it's the basis
for ODE23, the first

29
00:01:51,110 --> 00:01:57,100
of the methods we will use
out of the MATLAB ODE Suite.

30
00:01:57,100 --> 00:01:59,970
The basic method is order three.

31
00:01:59,970 --> 00:02:02,780
And the error estimate is
based on the difference

32
00:02:02,780 --> 00:02:06,710
between the order three method
and then the underlying order

33
00:02:06,710 --> 00:02:08,060
two method.

34
00:02:08,060 --> 00:02:12,170
There are four slopes involved.

35
00:02:12,170 --> 00:02:14,620
The first one is the
value of the function

36
00:02:14,620 --> 00:02:16,890
at the start of the interval.

37
00:02:16,890 --> 00:02:19,940
But that's based on
something called FSAL,

38
00:02:19,940 --> 00:02:24,930
first same as last, where
that slope is most likely

39
00:02:24,930 --> 00:02:27,750
left over from
the previous step.

40
00:02:27,750 --> 00:02:30,510
If the previous
step was successful,

41
00:02:30,510 --> 00:02:35,590
this function value is the
same as the last function

42
00:02:35,590 --> 00:02:38,530
value from the previous step.

43
00:02:38,530 --> 00:02:42,640
That slope is used to step into
the middle of the interval,

44
00:02:42,640 --> 00:02:45,330
function is evaluated there.

45
00:02:45,330 --> 00:02:48,480
That slope is used to
step 3/4 of the way

46
00:02:48,480 --> 00:02:53,730
across the interval and a
third slope obtained there.

47
00:02:53,730 --> 00:02:58,420
Then these three values
are used to take the step.

48
00:02:58,420 --> 00:03:01,960
yn plus 1 is a
linear combination

49
00:03:01,960 --> 00:03:04,980
of these three function values.

50
00:03:04,980 --> 00:03:09,850
Then the function is evaluated
to get a fourth slope

51
00:03:09,850 --> 00:03:11,820
at the end of the interval.

52
00:03:11,820 --> 00:03:17,750
And then, these four slopes
are used to estimate the error.

53
00:03:17,750 --> 00:03:24,550
The error estimate here is the
difference between yn plus 1

54
00:03:24,550 --> 00:03:29,400
and another estimate
of the solution

55
00:03:29,400 --> 00:03:32,680
that's obtained from
a second order method

56
00:03:32,680 --> 00:03:36,000
that we don't actually evaluate.

57
00:03:36,000 --> 00:03:39,730
We just need the difference
between that method and yn

58
00:03:39,730 --> 00:03:42,560
plus 1 to estimate the error.

59
00:03:42,560 --> 00:03:49,070
This estimated error is compared
with a user-supplied tolerance.

60
00:03:49,070 --> 00:03:53,170
If the estimated error
is less than a tolerance,

61
00:03:53,170 --> 00:03:55,370
then the step is successful.

62
00:03:55,370 --> 00:04:00,000
And this fourth
slope, s4, becomes

63
00:04:00,000 --> 00:04:03,010
the s1 of the next step.

64
00:04:03,010 --> 00:04:06,120
If the answer is bigger
than the tolerance,

65
00:04:06,120 --> 00:04:08,600
then the error
could be the basis

66
00:04:08,600 --> 00:04:10,800
for adjusting the step size.

67
00:04:10,800 --> 00:04:13,330
In either case,
the error estimate

68
00:04:13,330 --> 00:04:18,410
is the basis for adjusting the
step size for the next step.

69
00:04:18,410 --> 00:04:23,830
This is the Bogacki-Shampine
Order 3(2) Method

70
00:04:23,830 --> 00:04:27,970
that's the basis for ODE 23.

71
00:04:34,060 --> 00:04:41,030
Let's look at some very simple
uses of ODE23 just to get

72
00:04:41,030 --> 00:04:42,240
started.

73
00:04:42,240 --> 00:04:44,690
I'm going to take the
differential equation

74
00:04:44,690 --> 00:04:46,950
y prime is equal to y.

75
00:04:46,950 --> 00:04:51,500
So I'm going to
compute e to the t.

76
00:04:51,500 --> 00:04:58,720
And just call ODE23 on
the interval from 0 to 1,

77
00:04:58,720 --> 00:05:01,320
with initial value 1.

78
00:05:01,320 --> 00:05:03,180
No output arguments.

79
00:05:03,180 --> 00:05:08,210
If I call it ODE23, it
just plots the solution.

80
00:05:08,210 --> 00:05:09,000
Here it is.

81
00:05:09,000 --> 00:05:10,850
It just produces a plot.

82
00:05:10,850 --> 00:05:15,620
It picks a step size,
goes from 0 to 1,

83
00:05:15,620 --> 00:05:24,440
and here it gets the final
value of e-- 2.7 something.

84
00:05:24,440 --> 00:05:28,330
If I do supply output arguments.

85
00:05:28,330 --> 00:05:34,880
I say t comma y equals
ODE23, it comes back

86
00:05:34,880 --> 00:05:37,780
with values of t and y.

87
00:05:37,780 --> 00:05:41,810
ODE23 picks the
values of t it wants.

88
00:05:41,810 --> 00:05:43,810
This is a trivial problem.

89
00:05:43,810 --> 00:05:47,900
It ends up picking
a step size of 0.1.

90
00:05:47,900 --> 00:05:50,700
After it gets started, it
chooses an initial step size

91
00:05:50,700 --> 00:05:57,580
of .08 for whatever
error tolerances.

92
00:05:57,580 --> 00:06:05,470
And the final value of y is
2.718, which is the value of e.

93
00:06:05,470 --> 00:06:10,190
So these are the two
simple uses of ODE23.

94
00:06:10,190 --> 00:06:15,130
If you don't supply any output
arguments, it draws a graph.

95
00:06:15,130 --> 00:06:18,460
If you do supply output
arguments, t and y,

96
00:06:18,460 --> 00:06:22,610
it comes back with
the values of t and y

97
00:06:22,610 --> 00:06:27,160
choosing the values of
t to meet the error.

98
00:06:27,160 --> 00:06:31,090
The default error tolerances
is 10 to the minus 3.

99
00:06:31,090 --> 00:06:35,880
So this value is going to
be accurate to three digits.

100
00:06:35,880 --> 00:06:38,000
And sure enough
that's what we got.

101
00:06:42,379 --> 00:06:43,920
Now let's try
something a little more

102
00:06:43,920 --> 00:06:49,050
challenging to see the automatic
error-controlled step size

103
00:06:49,050 --> 00:06:51,690
choice in action.

104
00:06:51,690 --> 00:06:54,870
Set a equal to a quarter.

105
00:06:54,870 --> 00:07:00,940
And then set y0 equal to 15.9.

106
00:07:00,940 --> 00:07:05,440
If I would set it to 16,
which is 1 over a squared,

107
00:07:05,440 --> 00:07:08,190
I'd run into a singularity.

108
00:07:08,190 --> 00:07:11,240
Now the differential
equation is y

109
00:07:11,240 --> 00:07:18,580
prime is equal to 2 (a
minus t) times y squared.

110
00:07:18,580 --> 00:07:25,140
I'm going to integrate this
with the ODE23 on the interval

111
00:07:25,140 --> 00:07:33,010
from 0 to 1 starting at y0, and
saving the results in t and y,

112
00:07:33,010 --> 00:07:36,530
and then plotting them.

113
00:07:36,530 --> 00:07:42,810
So here's my plot command,
and there is the solution.

114
00:07:42,810 --> 00:07:47,770
So there is a near
singularity at a.

115
00:07:47,770 --> 00:07:49,490
It nearly blows up.

116
00:07:49,490 --> 00:07:51,730
And then it settles back down.

117
00:07:51,730 --> 00:07:56,270
So the points are
bunched together

118
00:07:56,270 --> 00:07:59,800
as you go up to the
singularity and come back down,

119
00:07:59,800 --> 00:08:04,340
but then get farther apart
as the solution settles down.

120
00:08:04,340 --> 00:08:11,050
And the ODE solver is
able to take bigger steps.

121
00:08:11,050 --> 00:08:14,060
To see what steps
were actually taken,

122
00:08:14,060 --> 00:08:19,985
let's compute the difference
of t, and then plot that.

123
00:08:25,620 --> 00:08:30,100
So here are the step
sizes that were taken.

124
00:08:30,100 --> 00:08:34,580
And we see that
a small step size

125
00:08:34,580 --> 00:08:40,980
was taken near the almost
singularity at that 0.25.

126
00:08:40,980 --> 00:08:44,380
And then as we get towards
the end of the interval,

127
00:08:44,380 --> 00:08:46,810
a larger step size is taken.

128
00:08:46,810 --> 00:08:49,820
And then, finally,
the step size just

129
00:08:49,820 --> 00:08:55,920
to reach the end of the interval
is taken as the last step.

130
00:08:55,920 --> 00:09:00,920
So that's the automatic
step size choice of ODE23.

131
00:09:06,230 --> 00:09:11,860
BS23 has a nice
natural interpolant

132
00:09:11,860 --> 00:09:14,590
that goes along with
it that's actually

133
00:09:14,590 --> 00:09:16,700
been known for over 100 years.

134
00:09:16,700 --> 00:09:21,510
It's called Hermite
Cubic Interpolation.

135
00:09:21,510 --> 00:09:25,850
We know that two points
determine a straight line.

136
00:09:25,850 --> 00:09:30,520
Well, two points and two
slopes determine a cubic.

137
00:09:30,520 --> 00:09:35,960
On each interval we have the
values of y and yn plus 1.

138
00:09:35,960 --> 00:09:38,710
We also have two
slopes, namely this.

139
00:09:38,710 --> 00:09:44,050
We have the derivatives at the
end points, yn prime and yn

140
00:09:44,050 --> 00:09:49,450
plus 1 prime, that's the
values of the differential

141
00:09:49,450 --> 00:09:51,490
equation at those points.

142
00:09:51,490 --> 00:09:55,000
So those four values
determine a cubic

143
00:09:55,000 --> 00:09:57,120
that goes through
those two points

144
00:09:57,120 --> 00:09:59,830
and has those two slopes.

145
00:09:59,830 --> 00:10:06,100
This cubic allows the software
to evaluate the solution

146
00:10:06,100 --> 00:10:10,390
at any point in the interval
without additional cost

147
00:10:10,390 --> 00:10:15,380
as defined by addition
evaluations of the function f.

148
00:10:15,380 --> 00:10:19,230
This can be used to draw
graphs of the solution,

149
00:10:19,230 --> 00:10:22,020
nice smooth graphs
of the solution,

150
00:10:22,020 --> 00:10:27,830
find zeroes of the solution,
do event handling, and so on.

151
00:10:27,830 --> 00:10:31,360
Another feature
provided by ODE23.