The first two commands you should learn if you are not already familiar with MATLAB are help and diary

A BRIEF GUIDE TO USING MATLAB IN MA2132

It is strongly recommended that you start using MATLAB as soon as you can. Not only because MATLAB assignments are part of your grade, but also as a tool to check your worksheet and practice homework solutions.

The first two commands you should learn if you are not already familiar with MATLAB are help and diary. The help command is used to give information about a MATLAB command. The diary command allows you to send your MATLAB output to a text file in the work folder, where you can write question numbers, edit out mistakes, add comments, put your solutions into a more readable form etc. This notes file has been made using the diary command and then edited in MS Word: MATLAB commands are in normal text, MATLAB output is indented and additional comments are in boldface. For example, in MATLAB type

help diary

and the MATLAB output is

DIARY Save text of MATLAB session.

DIARY filename causes a copy of all subsequent command window input

and most of the resulting command window output to be appended to the

named file. If no file is specified, the file 'diary' is used.

DIARY OFF suspends it.

DIARY ON turns it back on.

DIARY, by itself, toggles the diary state.

Use the functional form of DIARY, such as DIARY('file'),

when the file name is stored in a string.

There are several additional files you should download into your work folder which are needed in MA2132 and these are detailed elsewhere on the MA2132 website. For example, one of the first topics in the course is direction fields for first order differential equations. Once you have downloaded the m-file dfield6 (or dfield5 if you have MATLAB 5) you can graph direction fields by typing

dfield6

which opens a window in which you can input the differential equation. This in turn will open a graphics window with a graph of the direction field, and in which you can click on initial points to see specific particular solutions. You can also input initial conditions from the file menu. Similarly, you can download the file pplane6 and then use that command to graph vector fields for planar autonomous systems (this topic comes around halfway through the course).

**Note that you can not save graphics files within the file created by diary, you should print graphs separately.**

The main command for solving differential equations is dsolve. To find the general solution of the first order equation y’=t-y,

y=dsolve('Dy=t-y')

y =

t-1+exp(-t)*C1

To solve the initial value problem y’=t-y,y(0)=1

y=dsolve('Dy=t-y','y(0)=1')

y =

t-1+2*exp(-t)

The solution is stored as y, and you can evaluate y at t=3 using the subs (short for substitute) command

t=3;y3=subs(y)

y3 =

2.0996

Here, we have used y3 so as not to overwrite the function y.

To check or verify answers, you may need to integrate or differentiate which you can do with the functions int and diff. Note that if you do not assign a name to an operation, the output is stored as ans (short for answer) until another un-named operation is output which overwrites ans

int(y)

ans =

1/2*t^2-t-2*exp(-t)

diff(y)

ans =

1-2*exp(-t)

Sometimes an initial value problem has more than one solution. In this case the output is a vector

y=dsolve('Dy=-t/y','y(1)=0')

y =

[ (-t^2+1)^(1/2)]

[ -(-t^2+1)^(1/2)]

The two solutions are stored as y(1) and y(2). The MATLAB command plot is powerful but not necessarily easy or practical to use for ‘quick’ problems. Instead you could use the ezplot command to graph the second solution

ezplot(y(2))

One other thing to note about the dsolve command is that the default variable name is t. If you want to use a different variable name, you must specify that letter as an argument in dsolve. In this example z is a function of q.

z=dsolve('Dz-q^2*z=0','q')

z =

exp(1/3*q^3)*C1

Here is a more involved example. A rabbit population P(t) satisfies the logistic equation, where t is measured in months, and has an initial population of 50. After one month the population is 150, and after two months it is 250.

First we solve the initial value problem to find P in terms of the parameters K and r0

P=dsolve('DP=r0*(1-P/K)*P','P(0)=50')

P =

K/(1+1/50*exp(-r0*t)*(K-50))

Now use the values for P(1) and P(2) to find K and r0 using the solve command and store the answers in a vector

[K,r0]=solve('150=K/(1+1/50*exp(-r0)*(K-50))','K/(1+1/50*exp(-r0*2)*(K-50))=250')

K =

300

r0 =

log(5)

Substitute the values for K and r0 into the formula for P, the subs command remembers the values of letters found previously in the MATLAB session

P=subs(P)

P =

300/(1+5*exp(-log(5)*t))

You can use the simple and simplify commands to find a simpler form for P

P=simple(P)

P =

300/(1+5/(5^t))

The pretty command will (or should) give output as we would write it on paper

pretty(P)

300

--------

1 + ----

To find the population after three months

t=3;P3=subs(P)

P3 =

288.4615

and to plot the solution P(t) on the interval [0,10].

ezplot(P,[0,10])

The dsolve command can also be used to find the solutions of higher order differential equations. D2y is the second derivative of y, and in general Dny is the nth derivative of y where n is a positive whole number. To solve the IVP

y”+16y=3sin(4t), y(0)=0, y’(0)=0

y=dsolve('D2y+16*y=3*sin(4*t)','y(0)=0','Dy(0)=0')

y =

-3/32*cos(4*t)^2*sin(4*t)+(3/32*cos(4*t)*sin(4*t)-3/8*t)*cos(4*t)+3/32*sin(4*t)

y=simple(y)

y =

-3/8*cos(4*t)*t+3/32*sin(4*t)

In the second half of the course we are concerned with solving first order linear systems and to do this we must analyze matrices. You can find plenty of information concerning matrices in MATLAB at the MA2012 website, but here is a reminder of how to enter a matrix.

A=[-1 2 1; 0 -1 0; -1 -3 -3]

A =

-1 2 1

0 -1 0

-1 -3 -3

We are mainly interested in finding the eigenvalues and eigenvectors of A, and you can do this as follows

[v,d]=eig(A)

v =

0.7071 -0.7071 0.4082

0 0 0.4082

-0.7071 0.7071 -0.8165

d =

-2.0000 0 0

0 -2.0000 0

0 0 -1.0000

where v is a matrix whose columns are the eigenvectors and d is a matrix whose diagonal entries are the eigenvalues for A.

However, for our purposes, this probably isn’t the best command to use – firstly because it doesn’t address the issue of eigenvectors for repeated eigenvalues, and secondly because it normalizes the eigenvectors to have unit length whereas vectors with integer and/or fractional entries are more suitable (or easier) for our computations.

Instead, it is recommended that you use the eig command just to find the eigenvalues (which we will call lambda)

eig(A)

ans =

-2.0000

-1.0000

and then use the null command with the ‘r’ argument (this outputs vectors with integer and fraction entries) to find a basis for the nullspace of (A – lambda I) which is the same thing as finding eigenvectors for the eigenvalue lambda

N1=null(A+2*eye(3),'r')

N1 =

-1

This shows us that the eigenvalue -2 has algebraic multiplicity 2, and geometric multiplicity 1 because the basis of the nullspace (which is the space of eigenvectors) has dimension 1. To two generalized eigenvectors for -2 use the command

N2=null((A+2*eye(3))^2,'r')

N2 =

1 0

0 0

0 1

Then find the eigenvector for the eigenvalue -1

N3=null(A+eye(3),'r')

N3 =

-0.5000

1.0000

MATLAB by default assumes that letters represent numbers. If you want a letter to represent a variable, use the syms command.

syms t

Then, for example, to find the exponential of a square matrix use the expm command

E=expm(A*t)

E =

[ t*exp(-2*t)+exp(-2*t), -exp(-2*t)+exp(-t)+t*exp(-2*t), t*exp(-2*t)]

[ 0, exp(-t), 0]

[ -t*exp(-2*t), -t*exp(-2*t)+2*exp(-2*t)-2*exp(-t), -t*exp(-2*t)+exp(-2*t)]

The solution of the homogeneous, first order linear system IVP x’=Ax,x(0)=v is x=(e^tA)v. This gives us a quick way of finding/checking solutions of systems in MATLAB

v=[1;2;1]

v =

SysSol=E*v

SysSol =

[ 4*t*exp(-2*t)-exp(-2*t)+2*exp(-t)]

[ 2*exp(-t)]

[ -4*t*exp(-2*t)+5*exp(-2*t)-4*exp(-t)]

We can also find the solution of the system using the dsolve command. Remember that t is the default variable

[x,y,z]=dsolve('Dx=-x+2*y+z','Dy=-y','Dz=-x-3*y-3*z','x(0)=1','y(0)=2','z(0)=1')

x =

4*t*exp(-2*t)-exp(-2*t)+2*exp(-t)

y =

2*exp(-t)

z =

-4*t*exp(-2*t)+5*exp(-2*t)-4*exp(-t)

Sometimes matrices have complex eigenvalues and eigenvectors

B=[-1 3; -3 -1]

B =

-1 3

-3 -1

eig(B)

ans =

-1.0000 + 3.0000i

-1.0000 - 3.0000i

N=null(B-(-1+3*i)*eye(2),'r')

N =

0 - 1.0000i

1.0000

To find the numerical solutions of first order equations (and systems) with the Euler method and Runge-Kutta methods you need to download the m-files eul, rk2 and rk4 to your work folder as detailed elsewhere on the MA2132 website. Suppose we want to study the IVP y’=-y+t, y(0)=1. The first thing to do is create an m-file for the function f(t,y)=-y+t. First go to the file menu, then new, then m-file (or alternatively ctrl-n), and type

function yprime=f(t,y);

yprime=-y+t;

and save. This will by default save the function under the name f (of course if you use a different letter or name, it will save it under that name) into your MATLAB work folder. To use the Euler method, type

[teuler,yeuler]=eul('f',[0,1],1,.1)

This gives the t and y values of the Euler method for the function f. The second argument [0,1] is the interval

for which we want to know the approximate values of y, with 0 being the initial value for t. The third argument 1 is the initial value for y, so in other words y(0)=1. The fourth argument .1 is the stepsize h. The output is as follows (so for example, y(.8) is approximately .6609)

teuler =

0.1000

0.2000

0.3000

0.4000

0.5000

0.6000

0.7000

0.8000

0.9000

1.0000

yeuler =

1.0000

0.9000

0.8200

0.7580

0.7122

0.6810

0.6629

0.6566

0.6609

0.6748

0.6974

To use the second order Runge-Kutta method (also called the Improved Euler method) use the command rk2

[trk2,yrk2]=rk2('f',[0,1],1,.1)

trk2 =

0.1000

0.2000

0.3000

0.4000

0.5000

0.6000

0.7000

0.8000

0.9000

1.0000

yrk2 =

1.0000

0.9100

0.8381

0.7824

0.7416

0.7142

0.6988

0.6944

0.7000

0.7145

0.7371

To find the exact solution (and hence to see how accurate the numerical methods are) use the dsolve command

y=dsolve('Dy=-y+t','y(0)=1')

y =

t-1+2*exp(-t)

Another way to study the accuracy of these numerical methods is graphically. The plot command uses vectors to plot co-ordinates and the following command sets up two vectors: the first contains the values of t from 0 to 1 in steps of .1 (discover how the : command works), and the second contains the exact y values of the solution at those values of t.

t=0:.1:1;yexact=subs(y)

yexact =

Columns 1 through 5

1.0000 0.9097 0.8375 0.7816 0.7406

Columns 6 through 10

0.7131 0.6976 0.6932 0.6987 0.7131

Column 11

0.7358

Then you can use the plot command to graph the exact solution, the Euler solutions (marked with the ‘o’ symbol) and the RK2 solutions (marked with the ‘+’ symbol) on the same graph

plot(t,yexact,t,yeuler,'o', t, yrk2, '+')

The commands above are by no means exhaustive, nor are these methods the only “right” or “best” ways to answer questions. MATLAB is an extensive software package, and there are usually several ways to solve any problem. As always it is recommended that you practice with as many problems as possible and find your own preferred ways to do things. And don’t forget the help command!!