MATLAB for Math 308 (Ordinary Differential Equations) MTWRF 10:35 am BLOC 128
Instructor: Jean Marie Linhart              


Help with direction fields

Help with phase planes

Writing a function in MATLAB

Writing a vector function in MATLAB

Differential Equation Solvers

Plotting graphs


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

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.

Matrix Arithmetic

* is matrix/vector multiplication
+ and - give matrix/vector addition/subtraction
/ doesn't always do what you might think it should; caution! It works fine with scalars.

If you want MATLAB to give you elementwise operations try using a period . before the mathematical operator.
.* is elementwise multiplication
./ is elementwise division

Quick function plotting

ezplot('x*x') will give a plot of y=x2
ezplot('sin(x)') will give a plot of y = sin(x)

Direction Fields

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)

y' = y -x2
and plot its isoclines.

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

y' = y - x2= C
The isoclines are the equations
y = x2 + C
These are parabolas, all parallel to y = x2

To get the direction field in MATLAB

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:

Phase Plane/Direction Field Plots for 2-d Systems

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)

x' = y -x2
y' = x -y2
By hand we can calculate that this has 2 equilibrium points, one at (0,0) and another at (1,1); equilibrium points (also called critical points) are points where both x' and y' are equal to 0.

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:

Creating a Function in MATLAB

1-d functions

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

y' = y
We will call this function yp(x,y). Note that
y = ex
solves this ODE. We can get ex with the exp() function in MATLAB.

Vector functions

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

y'' - y = sin(t) with y(0) = 0; y'(0) = -1/2
I know the solution to this equation is y(t) = -sin(t)/2. I wrote this as a system of first order ODEs
x1 = y
x2 = x'1
x'2 = x1 + sin(t)

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.

Euler, Advanced Euler, Runge-Kutta Methods

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);
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);

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".

Plotting with MATLAB

MATLAB thinks in vectors and matrices, the natural objects to plot are vectors. If x, y,and t and you want to see plots of x and y vs t: plot(t, x, t, y);
If I prefer circles and dots to interconnected lines, I can use plot(t, x,'o', t, y,'.');

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;

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
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.