1
00:00:03,670 --> 00:00:05,250
PROFESSOR: Many
mathematical models

2
00:00:05,250 --> 00:00:07,750
involve high order derivatives.

3
00:00:07,750 --> 00:00:11,760
But the MATLAB ODE
solvers only work

4
00:00:11,760 --> 00:00:16,149
with systems of first
order ordinary differential

5
00:00:16,149 --> 00:00:17,720
equations.

6
00:00:17,720 --> 00:00:21,230
So we have to rewrite the
models to just involve

7
00:00:21,230 --> 00:00:23,380
first order derivatives.

8
00:00:23,380 --> 00:00:28,690
Let's see how to do that with a
very simple model, the harmonic

9
00:00:28,690 --> 00:00:30,290
oscillator.

10
00:00:30,290 --> 00:00:33,220
x double prime plus x equals 0.

11
00:00:33,220 --> 00:00:36,270
This involves a second
order derivative.

12
00:00:36,270 --> 00:00:39,040
So to write it as a
first order system,

13
00:00:39,040 --> 00:00:42,580
we introduced the vector y.

14
00:00:42,580 --> 00:00:45,740
This is a vector with
two components, x,

15
00:00:45,740 --> 00:00:47,660
and the derivative of x.

16
00:00:47,660 --> 00:00:52,830
We're just changing
notation to let y

17
00:00:52,830 --> 00:00:56,360
have two components,
x and x prime.

18
00:00:56,360 --> 00:01:00,340
Then the derivative
of y is the vector

19
00:01:00,340 --> 00:01:03,400
with x and x double prime.

20
00:01:03,400 --> 00:01:06,990
So the differential
equation now becomes

21
00:01:06,990 --> 00:01:11,750
y2 prime plus y1 equals zero.

22
00:01:11,750 --> 00:01:14,910
Do you see how we've just
rewritten this differential

23
00:01:14,910 --> 00:01:16,180
equation.

24
00:01:16,180 --> 00:01:26,750
so y2 prime is playing
of x double prime?

25
00:01:26,750 --> 00:01:29,515
Once you've done that,
everything else is easy.

26
00:01:32,290 --> 00:01:41,790
The vector system is now
y1, y2 prime is y2 minus y1.

27
00:01:41,790 --> 00:01:45,670
The first components
says y1 prime is y2.

28
00:01:45,670 --> 00:01:48,850
That's just saying
that the derivative

29
00:01:48,850 --> 00:01:51,580
of the first component
is the second.

30
00:01:51,580 --> 00:01:54,070
Here's the differential
equation itself.

31
00:01:54,070 --> 00:02:00,300
Y2 prime is minus y1 is the
actual harmonic oscillator

32
00:02:00,300 --> 00:02:01,370
differential equation.

33
00:02:04,370 --> 00:02:08,330
When we write this as an
autonomous function for MATLAB,

34
00:02:08,330 --> 00:02:11,130
here's the autonomous function.

35
00:02:11,130 --> 00:02:15,040
f is an autonomous
function of t and y,

36
00:02:15,040 --> 00:02:17,800
that doesn't depend upon t.

37
00:02:17,800 --> 00:02:22,340
First it's a vector
now, a column vector.

38
00:02:22,340 --> 00:02:25,600
The first component of f is y2.

39
00:02:25,600 --> 00:02:30,390
And the second component is -y1.

40
00:02:30,390 --> 00:02:34,450
The first component here is
just a matter of notation.

41
00:02:34,450 --> 00:02:38,717
All the content is in the second
component, which expresses

42
00:02:38,717 --> 00:02:39,800
the differential equation.

43
00:02:43,810 --> 00:02:46,120
Now for some
initial conditions--

44
00:02:46,120 --> 00:02:50,600
suppose the initial conditions
are that x of 0 is 0,

45
00:02:50,600 --> 00:02:53,990
and x prime of 0 is 1.

46
00:02:53,990 --> 00:02:57,940
In terms of the vector
y, that's y1 of 0,

47
00:02:57,940 --> 00:03:01,340
the first component of y is 0.

48
00:03:01,340 --> 00:03:04,230
And the second component is 1.

49
00:03:04,230 --> 00:03:10,140
Or in vector terms, the
initial vector is 0, 1.

50
00:03:10,140 --> 00:03:15,105
That implies they solution
is sine t and cosine t.

51
00:03:18,130 --> 00:03:22,350
When we write the initial
condition in the MATLAB,

52
00:03:22,350 --> 00:03:27,480
it's the column vector 0, 1.

53
00:03:27,480 --> 00:03:29,870
Let's bring up the
MATLAB command window.

54
00:03:29,870 --> 00:03:32,050
Here's the
differential equation.

55
00:03:32,050 --> 00:03:34,700
y1 prime is y2.

56
00:03:34,700 --> 00:03:37,010
And y2 prime is -y1.

57
00:03:37,010 --> 00:03:40,120
Here's the harmonic oscillator.

58
00:03:40,120 --> 00:03:42,350
We're going to
integrate from 0 to 2pi,

59
00:03:42,350 --> 00:03:45,230
because they're trig functions.

60
00:03:45,230 --> 00:03:51,140
And I'm going to ask for output
in steps of 2 pi over 36,

61
00:03:51,140 --> 00:03:54,440
which corresponds
to every 10 degrees

62
00:03:54,440 --> 00:03:58,370
like the runways at an airport.

63
00:03:58,370 --> 00:04:02,410
I'm going to need an
initial condition.

64
00:04:02,410 --> 00:04:06,680
y0 not is 0.

65
00:04:06,680 --> 00:04:15,170
I need a column vector, 0,
1, for the two components.

66
00:04:15,170 --> 00:04:18,300
I'm going to use
ODE45, and if I call it

67
00:04:18,300 --> 00:04:25,876
with no output arguments, ODE45
of the differential equation f,

68
00:04:25,876 --> 00:04:32,970
t span the time interval,
and y0 the initial condition.

69
00:04:32,970 --> 00:04:37,100
If I call ODE45 with
no output arguments,

70
00:04:37,100 --> 00:04:40,470
it just plots the
solution automatically.

71
00:04:40,470 --> 00:04:45,970
And here we get a graph
of cosine t starting at 1,

72
00:04:45,970 --> 00:04:50,390
and sine t starting at 0.

73
00:04:50,390 --> 00:04:53,910
Now if I go back to
the command window,

74
00:04:53,910 --> 00:05:01,400
and ask to capture
the output in t and y,

75
00:05:01,400 --> 00:05:08,390
I then get vectors of output.

76
00:05:08,390 --> 00:05:15,040
37 steps, vector t,
and two components

77
00:05:15,040 --> 00:05:20,790
y, the two columns
containing sine and cosine.

78
00:05:20,790 --> 00:05:24,130
Now I can plot them
in the phase plane.

79
00:05:24,130 --> 00:05:28,610
Plot ya against y2.

80
00:05:28,610 --> 00:05:37,390
If I regularize the axes, I
get a nice plot of a circle

81
00:05:37,390 --> 00:05:43,970
with small circles every
10 degrees, as I said,

82
00:05:43,970 --> 00:05:46,035
like the runways at an airport.

83
00:05:48,570 --> 00:05:52,660
The Van der Pol oscillator
was introduced in 1927

84
00:05:52,660 --> 00:06:00,080
by Dutch electrical engineer, to
model oscillations in a circuit

85
00:06:00,080 --> 00:06:02,730
involving vacuum tubes.

86
00:06:02,730 --> 00:06:07,620
It has a nonlinear damping term.

87
00:06:07,620 --> 00:06:11,200
It's since been used
to model phenomena

88
00:06:11,200 --> 00:06:13,470
in all kinds of
fields, including

89
00:06:13,470 --> 00:06:16,740
geology and neurology.

90
00:06:16,740 --> 00:06:20,980
It exhibits chaotic behavior.

91
00:06:20,980 --> 00:06:24,790
We're interested in it
for numerical analysis

92
00:06:24,790 --> 00:06:28,860
because, as the
parameter mu increases,

93
00:06:28,860 --> 00:06:34,500
the problem becomes
increasingly stiff.

94
00:06:34,500 --> 00:06:37,100
To write it as a first
order system for use

95
00:06:37,100 --> 00:06:42,840
with the MATLAB ODE solvers,
we introduce the vector y,

96
00:06:42,840 --> 00:06:46,890
containing x and x prime.

97
00:06:46,890 --> 00:06:52,460
So y prime is x prime
and x double prime.

98
00:06:52,460 --> 00:06:54,360
And then the
differential equation

99
00:06:54,360 --> 00:07:02,850
is written so that the first
component of y prime is y2.

100
00:07:02,850 --> 00:07:05,960
And then the
differential equation

101
00:07:05,960 --> 00:07:10,540
is written in the
second component of y.

102
00:07:10,540 --> 00:07:15,730
Here's the nonlinear
damping term minus y1.

103
00:07:15,730 --> 00:07:22,350
When mu is 0, this becomes
the harmonic oscillator.

104
00:07:22,350 --> 00:07:26,655
And here it is as the
anonymous function.

105
00:07:32,080 --> 00:07:36,160
Let's run some experiments with
the Van der Pol oscillator.

106
00:07:36,160 --> 00:07:39,850
First of all, I have to
define the value of mu.

107
00:07:39,850 --> 00:07:41,860
And I'll pick a
modest value of mu.

108
00:07:41,860 --> 00:07:43,770
Mu equals 100.

109
00:07:43,770 --> 00:07:49,860
And now with mu set, I can
define the anonymous function.

110
00:07:49,860 --> 00:07:52,100
It involves a value of mu.

111
00:07:52,100 --> 00:07:56,690
And here is the Van
der Pol equation

112
00:07:56,690 --> 00:07:59,820
in the second component of f.

113
00:07:59,820 --> 00:08:04,350
I'm going to gather statistics
about how hard the ODE

114
00:08:04,350 --> 00:08:06,270
solver works.

115
00:08:06,270 --> 00:08:09,210
And for that, I'm
going to use ODE set,

116
00:08:09,210 --> 00:08:11,940
and tell it I want
to turn on stats.

117
00:08:16,760 --> 00:08:18,265
I need an initial condition.

118
00:08:23,970 --> 00:08:28,030
Now I'm going to run
ODE45 on this problem.

119
00:08:28,030 --> 00:08:32,740
And I'm specifying just
a starting value of t,

120
00:08:32,740 --> 00:08:35,990
and a final value of t.

121
00:08:35,990 --> 00:08:41,490
ODE45 is going to pick
its own time steps.

122
00:08:41,490 --> 00:08:47,280
And I know with t
final equals 460,

123
00:08:47,280 --> 00:08:51,570
it's going to integrate
over it about two

124
00:08:51,570 --> 00:08:52,735
periods of the solution.

125
00:08:56,450 --> 00:09:01,890
Now we can watch
it go for a while.

126
00:09:01,890 --> 00:09:04,380
It's taking lots of steps.

127
00:09:04,380 --> 00:09:07,030
And it's beginning
to slow down as it

128
00:09:07,030 --> 00:09:08,270
takes more and more steps.

129
00:09:16,440 --> 00:09:21,060
Now this begins to get painfully
slow as it runs into stiffness.

130
00:09:25,660 --> 00:09:27,970
I'll turn off the
camera for a while here,

131
00:09:27,970 --> 00:09:30,010
so you don't have to
watch all these steps.

132
00:09:35,130 --> 00:09:38,140
We're trying to
get up here to 460.

133
00:09:38,140 --> 00:09:42,880
And I'll turn it back on
as we get close to the end.

134
00:10:07,620 --> 00:10:10,980
OK, well, the camera's been
off about three minutes.

135
00:10:10,980 --> 00:10:13,920
And you can see how
far we've gotten.

136
00:10:13,920 --> 00:10:16,050
We're nowhere near the end.

137
00:10:16,050 --> 00:10:20,450
So I'm going to press
the stop button here.

138
00:10:20,450 --> 00:10:29,660
And we'll go back to
the command window.

139
00:10:29,660 --> 00:10:35,440
And oh, so we didn't
get to the end here.

140
00:10:38,320 --> 00:10:42,166
Let me back off on
the time interval

141
00:10:42,166 --> 00:10:47,309
and try this value here.

142
00:10:47,309 --> 00:10:48,100
See how that works.

143
00:10:53,320 --> 00:10:56,650
So this is going to
still take lots of steps.

144
00:10:56,650 --> 00:11:04,810
But we'll be able to-- This
will go over about one period.

145
00:11:04,810 --> 00:11:07,400
We'll actually get
to the end here.

146
00:11:24,020 --> 00:11:26,750
I'll leave the camera
on until we finish.

147
00:11:44,480 --> 00:11:47,780
OK so that took a
little under a minute.

148
00:11:47,780 --> 00:11:52,390
And it took nearly 15,000 steps.

149
00:11:52,390 --> 00:11:54,355
So let's run it
with a stiff solver.

150
00:12:04,150 --> 00:12:06,550
There.

151
00:12:06,550 --> 00:12:11,670
So it took half a second,
and only 500 steps.

152
00:12:11,670 --> 00:12:18,800
So there's a modest
example of stiffness here.

153
00:12:18,800 --> 00:12:25,420
So let's examine the
Van der Pol equation

154
00:12:25,420 --> 00:12:28,520
using the stiff solver.

155
00:12:28,520 --> 00:12:30,930
Let's capture the
output and plot it.

156
00:12:30,930 --> 00:12:33,890
Because that plot
wasn't very interesting.

157
00:12:33,890 --> 00:12:36,780
I want to plot it a
couple of different ways.

158
00:12:36,780 --> 00:12:42,600
And again, I want to go
back up to the-- capture

159
00:12:42,600 --> 00:12:43,620
a couple periods.

160
00:12:51,130 --> 00:12:53,450
Let's plot one of the
current components

161
00:12:53,450 --> 00:12:54,695
as a function of time.

162
00:12:57,380 --> 00:12:58,130
And here it is.

163
00:12:58,130 --> 00:13:00,410
Here's the Van der Pol equation.

164
00:13:00,410 --> 00:13:03,840
And you can see it has
an initial transient,

165
00:13:03,840 --> 00:13:10,770
and then it settles into
this periodic oscillation

166
00:13:10,770 --> 00:13:20,400
with these really
steep spikes here.

167
00:13:20,400 --> 00:13:26,960
And even this stiff
solver is working hard

168
00:13:26,960 --> 00:13:28,610
at these rapid changes.

169
00:13:28,610 --> 00:13:30,960
We've got a fair
number of points

170
00:13:30,960 --> 00:13:35,720
in here, as it is it
chooses the step size.

171
00:13:35,720 --> 00:13:42,550
And now, let's go back
to the command line

172
00:13:42,550 --> 00:13:47,502
and do a phase plane plot.

173
00:13:51,230 --> 00:13:55,650
So here's the phase plane
plot of this oscillator

174
00:13:55,650 --> 00:13:56,750
with damping.

175
00:13:56,750 --> 00:14:02,840
And it's nowhere near a circle,
which it would be if mu was 0.

176
00:14:02,840 --> 00:14:09,390
And this is the characteristic
behavior of the Van der Pol

177
00:14:09,390 --> 00:14:10,940
oscillator.