|
MATLAB for Math 308 (Ordinary
Differential Equations) MTWRF 10:35 am BLOC 128 | ||||||
MATLAB Support from MathWorks Sometimes it is easier to Google Search rather than going directly to the MathWorks site.
MATLAB setup for Math 308 students
MATLAB naturally understands vectors and matrices. Type
x = [1 2 3]
at the MATLAB prompt, and MATLAB will understand that x is a
row vector with 3 elements as given, and it will echo the row vector
back to you.
If you type
x = [1;2;3]
Matlab now understands that x is a column vector with those 3
elements, and MATLAB will echo the column vector back to you.
If you don't want output echoed back to you, end your line with a
semi-colon:
;
x = [1;2;3];
If you type
A = [1, 0, -1; 0, 1, 1; 1, -1, 0];
MATLAB now understands A is a 3x3 matrix where the first row is
1, 0, -1, the second row is 0, 1, 1, and the final row is 1, -1, 0.
Since there is a semi-colon at the end of this command, MATLAB will
not echo the matrix back to you. To see it you could just type
M
x = -3:.1:3; fills x with a vector -3, -2.9, etc
x = [1,2,3; -3, -2, -1; 2, -1, 3] makes a 3x3 matrix
The colon operator and submatrix indexing:
x(:,1) gives the first column of x; read this as all rows and
column 1.
x(1:2, 1) gives the first two elements onf the first column of
x; read this as row 1-2 of column 1.
x(1,:) gives the first row of x; read this as row 1, all columns.
If you want MATLAB to give you elementwise operations try using a
period . before the mathematical operator.
.* is elementwise multiplication
./ is elementwise division
Warning: I often say "differential field" when I mean "direction field" -- they are, technically, not the same thing.
To plot direction fields we are using dfield8 in class. In the Polking book it is called dfield6. I will use dfield8 here and trust you can translate.
Type dfield8 at the MATLAB prompt to bring up the dialog box.
Let's say we want to look at the direction field for the differential equation (DE)
An isocline to the DE is a set of points on which all solutions of the DE have the same slope. For this equation it means to find the points in (x,y) space where
To plot an isocline on the direction field:
To make it easier to see that the slope is zero along the isocline, we are going to plot some solutions to the differential equation.
To plot a solution to the DE on the direction field:
To remove isoclines from your plot:
To remove solutions from your plot:
To plot direction fields for 2-d systems of differential equations, we really need to be working with an autonomous system. To do this we will use pplane8 in class. In the Polking book it is called pplane6.
Type pplane8 at the MATLAB prompt to bring up the dialog box.
Let's say we want to look at the direction field for the system of differential equations (DE)
To get the phase plane in MATLAB:
To plot an equilibrium value, especially one you've already computed:
To plot a trajectory (also called an orbit) to the DE system on the phase plane:
To remove solutions from your plot:
To use many of these methods, it is wise to know how to define a function in an M-file. An M-file is a MATLAB file containing commands or a function or a program.
I am going to define a function M-file for
You can also create functions that return a matrix or a vector with MATLAB.
Let's say I want to program a function to return
x1' = x2;
x2' = x1 + sin(t)
The reason for this choice is that these equations comes from the
ODE system
Now I write a function M-file for this system. We will call this function xprime(t,x). Here x is a 2x1 column vector.
Now I am going to use the Euler method. I will use my scalar
function yp() that I defined above; I could also
use xprime my vector function. The Euler method is
implemented as the eul routine with syntax
[t,y] = eul(@yp, tspan, y0, stepsize)
The stepsize argument is optional. Let us say I want to evaluate this on [0,1] Recall e1 = e = 2.71828....
If I omit the stepsize argument, it will use 100 equally spaced points to calculate, e.g. in this case the stepsize will be 0.01
We will start with a stepsize of 0.1, and plot our results with MATLAB's plot command.
At the MATLAB prompt enter the following commands. Omit the ";" at the end
of the commands if you want to see output.
MATLAB thinks in vectors and in matrices, so the output here t will be
a vector containing the steps of size 0.1 from 0 to 1. The y vector
will contain the approximations for y(x) where y(0) = 1. Of course
you'll want to visualize your data with a plot.
[t,y] = eul(@yp,[0,1],1,.1);
yexact = exp(t);
plot(t,y,t,yexact,'.');
The last entry here, the '.', tells MATLAB to plot yexact vs. t using dots
at the t values instead of connecting the points with lines.
Let's say I want to compare doing this with a stepsize of 0.1 and
a stepsize of 0.05. Here's what I would do to approximate with the Euler
method and plot my two approximations.
[t1,y1] = eul(@yp, [0,1],1,.1);
[t05,y05] = eul(@yp,[0,1],1,.05);
plot(t1,y1,t05,y05);
If I want to plot the exact solution calculated at the t1 points, I can use
the command
yexact = exp(t1);
And then I plot
plot(t1,y1, t05,y05,t1,yexact);
or
plot(t1,y1,t05,y05,t1,yexact,'.');
to plot the exact solution with dots instead of lines.
You should observe the smaller the stepsize the better the final result.
The Runge-Kutta routines rk2() and rk4() have the same syntax as eul(). So
[t,yrk4] = rk4(@yp,[0,1],1,.1);
Uses the Fourth Order Runge-Kutta method on [0,1] with a step size
of 0.1 to calculate.
Using the vector function xprime() with eul() or rk4 works just as before. Now the Euler method will return t, our time vector and x, a matrix with two columns for x1 and x2.
Let's say I want to evaluate this from [0,4 pi] with a stepsize of
4 pi/25. Recall the initial condition was y(0) = 0, y'(0) = -0.5, which
is x1 = 0 and x2 = 1, or x0 = [0,1].
We enter
[t, x] = eul(@xprime, [0, 4*pi], [0,-0.5], 4*pi/25);
Our solution is now in 2 columns in x. We can do component plots of
each component versus time
plot(t,x(:,1),t, x(:,2), '.');
Here x(:,1) says give me all rows and the first column. x(:,2) says
all rows and the 2nd column. Unfortunately
neither of these graphs look like they should...
I figure my stepsize is too big, so I'm going to try this again
with a much smaller stepsize.
[t, x] = eul(@xprime, [0, 4*pi], [0,-0.5], 4*pi/250);
Let's try plotting the components again
plot(t,x(:,1),t, x(:,2), '.');
But that doesn't look right.
I am going to try this again with the Runge-Kutta 4 routine, which
I know is a lot more accurate:
[t, x] = rk4(@xprime, [0, 4*pi], [0,-0.5], 4*pi/250);
plot(t,x(:,1),t, x(:,2));
Finally! That looks mostly like what I was expecting.
What is going on here? This is numerical stability.
We can run into stability problems with numerical routines. Short of doing a detailed stability analysis, we should always use our common sense for what a solution should look like, and evaluate what happens as we change the stepsize. Be very careful before drawing conclusions when you don't know what your equation should look like or how good your numerical routine is.
One last ODE solver I want to introduce is MATLAB's built in ode45
routine. This is a Runge-Kutta routine. Its syntax is
[t,y] = ode45(@yp, timevector, y0)
The timevector is a vector of time values. To generate it,
we can use the linspace(a,b, N) routine that generates N evenly
spaced points between a and b.
times = linspace(0,4*pi, 101);
[t, x] = ode45(@xprime, times, [0, -0.5]);
plot(t, x(:,1), t, x(:,2))
If you want to know about MATLAB's ODE solvers try Googling on
"MATLAB ODE solver".
If I have my data in a matrix x and I want to plot the first
two columns vs t
plot(t, x(:,1),'o', t, x(:,2) ,'.');
Or I can look at how they trace through 3-d space, with t for
time as the z-axis:
plot3(x(:,1),x(:,2), t);
If I want to do a function plot, I can evaluate the function and
put the data into matrices. For example, to get overlaid plots of y
= x2, y = x2 -1, y = x2 -2, etc. on
-3 <= x <= 3
x = -3:.1,:3;
y1 = x.*x;
y2 = x.*x - 1;
y3 = x.*x - 2;
plot(x,y1,x,y2,x,y3);
You can get overlaid 3d plots as well. Now you use plot3. Let's get
a 3-d spiral. Let x1 = sin(t), y1 = cos(t) and z1 = t, and let x2 = cos(t),
y2 = sin(t) and z2 = z1 = t.
t = 0:pi/25:4*pi;
x1 = sin(t);
y1 = cos(t);
x2 = cos(t);
y2 = sin(t);
plot3(x1, y1, t, x2, y2, t);
Note that the above is equivalent to:
plot3(x1, y1, t, y1, x1, t);
I can throw some options in to get a fancier graph. I use hold to
keep plotting things in the same figure window.
plot3(x1, y1, t, 'b-'); where 'b-' says to use blue lines
hold;
plot3(y1, x1, t, 'r-'); where 'r-' says to use red lines
hold off; releases the graph.
(C) 2009/2011 Jean Marie Linhart, all rights reserved.