MATLAB for Math 308 (Ordinary Differential Equations) MTWRF 10:35 am BLOC 128 Instructor: Jean Marie Linhart MATLAB Basics Help with direction fields Help with phase planes Differential Equation Solvers Plotting graphs

## MATLAB Basics

MATLAB Support from MathWorks Sometimes it is easier to Google Search rather than going directly to the MathWorks site.

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.

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

• Type dfield8 at the MATLAB prompt to bring up the dfield8 dialog box.
• In the boxes for "The differential equation" put in y' = y -x^2 = C
• Fill in the independent variable (you know what this is!)
• We leave "parameters and expressions" blank.
• Choose an appropriate range for x and y. Experiment! Start with -2<= x <= 2 and -4 <= y <= 4
• Hit "proceed" when you are ready to see the direction field.

To plot an isocline on the direction field:

• Go to the graph of the direction field (dfield8 Display), and choose Options-->Plot Level Curves (this opens up a dialog box)
• In the "dfield8 Level sets" dialog box, enter the isocline. Ours is y - x^2.
• If we want to plot curves for the isocline with C = 0, choose the radio button "Enter a vector of level values" and enter 0. This will give the isocline y = x2 with slope = C = 0 along the isocline.
• hit "proceed" to plot the level curve.
• This is the isocline along with the slope (direction field) is zero. Can you see this? Maybe not so easily.

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:

• Go to the graph of the direction field (dfield8 Display)
• Choose Options-->Keyboard Input
• In the Keyboard input dialog, input some data along your isocline, e.g. (0,0), (1,1), (2,4), (-1,1) ... and hit Compute.
• You can also click with your mouse on the graph to make solutions.
• Now you should easily be able to see that where these curves cross your isocline y=x2, they have zero slope.

To remove isoclines from your plot:

• Go to the graph of the direction field (dfield8 Display)
• Choose Edit-->Erase all level curves

To remove solutions from your plot:

• Go to the graph of the direction field (dfield8 Display)
• Choose Edit-->Erase all solutions

## 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:

• Type pplane8 at the MATLAB prompt to bring up the dfield8 dialog box.
• In the boxes for "The differential equations" put in x' = y -x^2 and y' = x - y^2.
• We leave "parameters and expressions" blank.
• Choose an appropriate range for x and y. Since our critical points are at (0,0) and (1,1), we want our range to include these points. Experiment! For this problem a good start is something like -4 <= x <= 4 and -4 <= y <= 4 which will encompass both and a little bit more.
• Hit "proceed" when you are ready to see the phase plane.

To plot an equilibrium value, especially one you've already computed:

• Go to the graph of the phase plane (pplane8 Display).
• Choose Solutions-->Find an equilibrium point
• Click on the graph near where the equilibrium point should be.
• The equilibrium point will appear on your graph in red.
• A dialog box with information about the equilibrium point will also appear.

To plot a trajectory (also called an orbit) to the DE system on the phase plane:

• Go to the graph of the phase plane (pplane8 Display)
• You can click on any location on your graph to get a solution for that point.
• If you want to graph a solution through a particular point:
• Choose Solutions-->Keyboard Input
• In the Keyboard input dialog, input the point you want. Try (-1,0), (0.5, 1) (4,-1) for example ... and hit Compute.

To remove solutions from your plot:

• Go to the graph of the phase plane (pplane8 Display)
• Choose Edit-->Erase all solutions

## 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.
• Go to File-->New --> Function M-file
• Replace output_args with yprime, untitled with yp, and input_args with x, y. You are defining yprime = yp(x,y).
• Replace %UNTITLED Summary of the function goes here with % yp: y' = y Here % indicates a comment.
• If you want you can replace the next % line with another comment.
• Below the comment lines add
yprime = y;
This is your function. The ; indicates that you do not want this command to show its output.
• Save the file as yp.m. Remember: yp is the name of your function -- the file name must agree with the name of the function you gave in the M-file!.
• Let me emphasize this. Always check your function!
• yp(2,4)
• yp(1,5)
• yp(10,2)
• yp(2,10)
Did you get what you expected? If not, you need to go back and correct your work.
• I can now use my function yp with any of the ODE solvers available to me 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.

• Go to File --> New --> Function M-file
• Replace output_args with xp, untitled with xprime, and input_args with t, x. You are defining xp = xprime(t,x).
• Replace
%UNTITLED Summary of the function goes here
with
% xprime: x' = f(t,x) where x is vector.
Here % indicates a comment.
• If you want you can replace the next % line with another comment describing what we are doing. Or make even more comment lines! It is good to have a thorough explanation of what's going on in the comments so a future programmer can figure out what is going on.
• Below the comment lines add
xp = zeros(2,1);
xp(1) = x(2);
xp(2) = x(1) + sin(t);
The first line initializes the output to be a 2x1 vector of zeros. The second line defines the first component of the vector. The third line defines the second component of the vector.
• Save the file as xprime.m (xprime is the name of your function -- this must agree with the name of the function you gave in the M-file. Now in your MATLAB workspace, check your function.
• Let me emphasize this. Always check your function! The arguments to xprime are a number, t, and a vector x, with two entries, we input x as [a, b] where a = x(1) and b = x(2).
• xprime(0,[1,1]) should return a column vector with 1, 1 entries.
• xprime(pi, [0,0]) should return a column vector with 0, 0 entries.
• xprime(pi/2, [1,0]) should return a column vector with 0, 2 entries
Did you get what you expected? If not, you need to go back and correct your work.
• Now I can use my function xprime with any of the ODE solvers available to me in MATLAB!

## 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);
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".

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