] > Matlab primer

## A VERY BASIC MATLAB PRIMER

Date: Last updated November 2009.

Some of the examples in this tutorial are the modified versions of the items from Rudra Pratap’s, Getting started with Matlab, Oxford University Press, 2002.

### 1. A Few Elementary Calculations

>> 2+4
ans =
6

Enter $2+4$ and hit return/enter key.
>> x = 2+5
ans =
7

Assigning the value of an expression to a variable.
>> z= 3^4+log(pi)*cos(x);
>> z

81.8630

A semicolon suppresses screen output. You can recall the value of $z$ by typing $z$
>> format short e
>> z
z =
8.1863e+01

The floating point output display is controlled by the format command.
>> format long
>> z
z =
81.863014441556018
>> format short
>> z
z =
81.8630
>> quit

Quit Matlab. Alternatively you can type exit command or select Exit Matlab from the File Menu.

For additional help type in Matlab: help format

### 2. Finding Help in Matlab

There are several ways to get help in Matlab:

• help functionname provides direct on-line help on the function whose eaxct name you know. For example help help provides help on the function help.
• lookfor keyword provides a list functions with the string keyword in their description. For example,
>> lookfor prime
isprime                        - True for prime numbers.
factor                         - Prime factors.
primes                         - Generate list of prime numbers.

By clicking on isprime, factor, or primes additional help on those functions will be listed.

• helpwin provides the click and navigate help. A new help windows opens. One can also select help choices from the Help menu.
• helpdesk is the web browser-based help, including printable (PDF) documentation on the Web.
• demo provides the demonstration programs on Matlab.

### 3. Working with Arrays of Numbers

>> x = [1 2 3]
x =
1     2     3
>> size(x)
ans =
1     3

$x$ is a row vector with $3$ elements.

Size of array. Here, 1 row and 3 columns.

>> y = [4; 5; 6]
y =
4
5
6
>> size(y)
ans =
3     1

$y$ is a column vector with $3$ elements.

Size of y: 3 rows and 1 column.

>> x’
ans =
1
2
3

x’ creates a column vector from x.
>> b = [3 -1 0]
b =
3    -1     0
>> a = x+b
a =
4     1     3

Addition of two vectors of the same size. The same holds for subtraction.
>> x+y
??? Error using ==> plus
Matrix dimensions must agree.

Addition or subtraction of a row vector to a column vector is not allowed.
>> c = x.*a
c =
4     2     9
>> d = x./a
d =
0.2500    2.0000    1.0000
>> x*a
??? Error using ==> mtimes
Inner matrix dimensions must agree

One can multiply (or divide) the elements of two same-sized vectors term by term with the array operator: .* (or ./). The array operator . ̂ produces term by term exponentiation.
However, x*a produces an error.
>> X = linspace(-2,7,4)
X =
-2     1     4     7
>> X = -2:3:7
X =
-2     1     4     7

linspace(a,b,n) creates a linearly spaced vector of length $n$ from $a$ to $b$. Note that u=linspace(a,b,n) is the same as u=a:(b-a)/(n-1):b, another way of creating vectors (see also Section 7).
>> y = cos(X);
>> z = sqrt(abs(X)).*y
z =
-0.5885    0.5403   -1.3073    1.9946

Elementary math functions operate on vectors term by term.

### 4. Creating Script Files and Simple Plots

A script file is an external file that contains a sequence of Matlab statements. By typing the filename, subsequent Matlab input is obtained from the file. Script files have a filename extension of .m and are often called M-files.

Select New from File menu (alternatively, use any text editor) and type the following lines into this file:

% ELLIPSE_PARABOLA - A script file to draw ellipse (x/2)^2+(y/3)^2=1 and parabola y=x^2.
alpha = linspace(0,2*pi,100);
x = 2*cos(alpha);
y = 3*sin(alpha);
t = linspace(-2,2,100);
z = t.^2;
plot(x,y,t,z);
axis(’equal’);
xlabel(’x’)
ylabel(’y’)
title(’Ellipse (x/2)^2+(y/3)^2=1 and parabola y=x^2’)
clear alpha x y t z;

Lines starting with % sign are ignored; they are interpreted as comment lines by Matlab.

Select Save Asfrom File menu and type the name of the document as, e.g., ellipse_parabola.m. Alternatively, if you are using an external editor, use a suitable command to save the document in the file with .m extension.

Once the file is placed in the current working directory of Matlab, type the following command in the command window to execute ellipse_parabola.m script file.

>> ellipse_parabola

Execute the file. You should see the plot of the ellipse ${\left(\frac{x}{2}\right)}^{2}+{\left(\frac{y}{3}\right)}^{2}=1$ and parabola $y={x}^{2}$ (see Figure 1, below).

>> print

Print the created graph on the default printer.
>> print my_graph.eps -depsc

Save the graph in (color) Encapsulated PostScript format in the file my_graph.eps. Use -djpeg, -dtiff, or -dpng options for JPEG, TIFF, or PNG graphic formats, respectively.

The above script file also illustrates how to generate plot $x$ vs $y$ from the generated $x$ and $y$ coordinates of 100 equidistant points between $0$ and $2\pi$. For example, in the above script change the entry plot(x,y,t,z) to plot(x,y,’--’,t,z,’o’) and see the difference in graphs drawings. For more information, check Matlab help on plot command and its options by invoking help plot in the command window of Matlab. For print options enter help print in the command window of Matlab.

### 5. Creating and Executing a Function File

Function file is a script file (M-file) that adds a function definition to Matlab’s list of functions. Edit ellipse_parabola.m script file (listed in Section 4 to look like the following:

function circle_f1(r);
% CIRCLE_F1 - Function to draw circle x^2+r^2=r^2.
alpha = linspace(0,2*pi,100);
x = r*cos(alpha);
y = r*sin(alpha);
plot(x,y);
axis(’equal’);
xlabel(’x’)
ylabel(’y’)
title([’Circle of radius r =’,num2str(r)]) % put a title with the value r.

Save it in circle_f1.m file (remember to place it in the working directory of Matlab) and try it out.

Specify the input and execute the function. You should see the graph of the circle with radius 3. You can also just enter circle_f1(2.4) to draw the circle with radius $r=2.4$.

And here is slightly different version of circle_f1.m function file, called circle_f2.m This function file, in addition to drawing a circle with radius $r$, computes the area inside this circle together with the outputs for x and y coordinates.

function [A,x,y] = circle_f2(r);
% CIRCLE_F2 - Function draws circle x^2+r^2=r^2
% and output its area together with x and y coordinates.
alpha = linspace(0,2*pi,100);
x = r*cos(alpha);
y = r*sin(alpha);
A = pi*r^2;
plot(x,y);
axis(’equal’);
xlabel(’x’)
ylabel(’y’)
title([’Circle of radius r =’,num2str(r)]) %put a title with the value r.

Save it in circle_f2.m file (remember to place it in the working directory of Matlab) and execute it in one of the following ways:

>> circle_f2(5)
>> [Area] = circle_f2(5)
>> circle_f2(5);
>> [Area] = circle_f2(5);

The first two commands output the graph of a circle with radius $5$ and its area. The last two commands create the graph while suppressing the output for the area.
>> [AREA,X,Y] = circle_f2(5)

The complete output of the function is displayed: the graph of a circle with radius $5$ and its area, together with the outputs for x and y coordinates.

Finally, suppose you need to evaluate the following function many times at different values of $x$ during your Matlab session: $f\left(x\right)=exp\left(-0.3x\right)sin\left(30x\right)+\frac{{x}^{10}}{1+{x}^{10}}$. You can create an inline function and pass its name to other functions:

>> fun = inline(’exp(-0.3*x).*sin(30*x)+x.^10./(1+x.^10)’);
>> fun(2)
ans =
0.8317
>> sin(fun(1))+fun(cos(1))
ans =
-0.6363
>> x=linspace(0,10,100);
>> plot(x,fun(x))
>> x0=fzero(fun,10)
x0 =
1.0193

Use Matlab help command help fzero to check what this last command computes.

In the case you want to use function $f\left(x\right)=exp\left(-0.3x\right)sin\left(30x\right)+\frac{{x}^{10}}{1+{x}^{10}}$ in future Matlab sessions, create the following function file:

function f = my_f(x);
% MY_F evaluates function f = exp(0.3x)sin(30x)+x^10/(1+x^10)
f = exp(-0.3*x).*sin(30*x)+x.^10./(1+x.^10);

and save it in my_f.m file (remember to place it in the working directory of Matlab). You can use it in many ways:

>> my_f(2)
ans =
0.8317
>> sin(my_f(1))+my_f(cos(1))
ans =
-0.6363
>> x=linspace(0,10,100);
>> plot(x,my_f(x))
>> x0=fzero(@my_f,10)
x0 =
1.0193

Use the command help function_handle to check what is the meaning of @ symbol in the last command.

And here is Matlab function that for a given $x\in ℝ$ and given $n\in ℕ$ computes

$1+x+\frac{{x}^{2}}{2!}+\frac{{x}^{3}}{3!}+\cdots +\frac{{x}^{n}}{n!},$

function s=expseries(x,n);
% EXPSERIES: function calculates the sum
% 1+x+(x^2)/2!+(x^3)/3!+...+(x^n)/n! (up to nth power).
% Call syntax:
%         s=expseries(x,n);
% where x is real and n nonnegative integer.
k = 0:n;                         % create a vector from 0 to n
series = (x.^k)./factorial(k);   % create a vector of terms in the series
s = sum(series);                 % sum all terms of the vector ’series’.

Save the script in expseries.m file (remember to place it in the working directory of Matlab) and try it out. Alternatively, you can also download it from:

### 6. Saving/Loading Data, Recording a Session, Clearing Variables, and More on Plotting

>> save

Saves all workspace variables in the file matlab.mat.
>> save FILENAME.dat

Saves all workspace variables in the file FILENAME.dat. If FILENAME has no extension, .mat is assumed.
>> save X

Saves variable X in the file X.mat.
>> save FILENAME.dat X

Saves variable X in the file FILENAME.dat. If FILENAME has no extension, .mat is assumed.

Loads the variables saved in the default file matlab.mat.

Loads the variables saved in the file FILENAME.mat.
>> 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 diary.
>> diary on

Turns diary back on. The command diary, by itself, toggles the diary state.
>> clear X
>> clear all

Removes variable X from memory.
clear all command removes all variables, functions and MEX links. Also clear all at the command prompt removes the Java packages import list. It is a good idea to start a Matlab session with this command !

In addition to plot command, introduced in Section 4, Matlab has many plotting utilities. Here are sample uses of two of them.

fplot command

>> fplot(’exp(-x).*sin(30.*x)’,[-1,5])
>> xlabel(’x’),ylabel(’f(x)=e^{-x} sin(30x)’)
>> title(’How to use fplot’)

Plots a function between $-1$ and $5$. The last two commands are optional and only annotate the graph with the title and other details.
>> F = ’exp(-x).*sin(30.*x)’;
>> fplot(F,[0,10])
>> F = inline(’exp(-0.3*x).*sin(30*x)’);
>> fplot(F,[0,10])
>> fplot(’[my_f(0.5*x),cos(10*x),sin(x)^2]’,[-7,7])

Other uses of fplot, including multiple plots. Here, my_f is the function file defined in Section 5.
>> [X,Y] = fplot(’exp(-x).*sin(30.*x)’,[-1,1])
>> X(1)
ans =
-1
>> Y(1)
ans =
2.6857

Returns X and Y such that $Y=exp\left(-X\right)sin\left(30X\right)$ (output not provided here). No plot is drawn on the screen. One can check the length of vectors $X$ and $Y$ with length(X) and length(Y) commands. Also, X(i) or Y(k) provide $i$-th and $k$-th elements of vectors $X$ and $Y$, respectively.

ezplot command

>> F = ’exp(-x).*sin(30.*x)’;
>> ezplot(F)
>> F = inline(’exp(-0.3*x).*sin(30*x)’);
>> ezplot(F)
>> ezplot(’my_f’)
>> ezplot(’my_f(3*x^2)’)

Plots the function over the default domain $-2\pi and with the default title. Also, ezplot(fun,[min,max]) plots fun(x) over the domain: min < x < max.
>> ezplot(’x^2*y^2-sin(x^2+y^2)=1’)
>> ezplot(’t*sin(t)’,’1+cos(2*t)^2’)

Additionally, ezplot can be used for implicitly defined functions and parametrically defined curves. For parametric curves $\left(x\left(t\right),y\left(t\right)\right)$, the default domain is $0.

For other plotting utilities check Matlab help on the following commands:
ezcontour, ezcontourf, ezmesh, ezmeshc, ezplot3, ezpolar, ezsurf, ezsurfc.

### 7. Matrices and Vectors

>> A=[0 1 2; 3, 4, 5; 6 7 8]
A =
0     1     2
3     4     5
6     7     8

While entering matrices, rows are separated by semicolons and columns are separated by spaces or commas.
>> A=[0 1 2
3 4 5
6 7 8]
A =
0     1     2
3     4     5
6     7     8

>> A=[0 1 ...
2; 3 4 5; 6 7 8]
A =
0     1     2
3     4     5
6     7     8

A matrix can be also entered across multiple lines using carriage returns at the end of each row.
Three consecutive periods (ellipsis) signal continuation of the input on the next line.
>> A(3,2)
ans =
7

>> A(3,:)
ans =
6     7     8

>> A(:,2)
ans =
1
4
7

>> A(:,2:3)
ans =
1     2
4     5
7     8

Element ${A}_{ij}$ of matrix $A$ is accessed by A(i,j).
The colon by itself as a row or column index specifies all rows or columns of the matrix.
>> Z = [-2.4 1 5 2.6 0]
Z =
-2.4000  1.0000  5.0000  2.6000  0

>> D=[1;-3;sin(1);exp(1)]
D =
1.0000
-3.0000
0.8415
2.7183

$Z$ is a row vector with $5$ elements.
$D$ is a column vector with $4$ elements.
>> v=1:10:111
v =
1    11    21    31    41    51    61    71    81    91   101   111

The command v=InitValue:Increment:FinalValue
creates a vector over a given range with a specified increment,
>> x=-2:10

x =
-2    -1     0     1     2     3     4     5     6     7     8    9    10

while x=InitValue:FinalValue produces a vector over a given range with increment $1$.

The command linspace(a,b,n) creates a linearly spaced (row) vector of length $n$ from $a$ to $b$. Note that u=linspace(a,b,n) is the same as the command u=a:(b-a)/(n-1):b.

Also, logspace(a,b,n) generates a logarithmically spaced (row) vector of length $n$ from $1{0}^{a}$ to $1{0}^{b}$. For example,

>> v=logspace(0,4,5)
v =
1          10         100        1000       10000

Finally, logspace(a,b,n) is the same as 10.^(linspace(a,b,n) where .^ is term by term operator of exponentiation (see Section 7.3).

There are many ways to create special vectors, for example, vectors of zeros or ones of specific length:

>> u=zeros(1,5)
u =
0     0     0     0     0

>> u=ones(2,1)
u =
1
1

7.1. Appending or Deleting a Row or Column.

A row/column can be easily appended or deleted, provided the row/column has the same length as the length of the rows/columns of the existing matrix.

>> C=eye(3)
C =
1     0     0
0     1     0
0     0     1

The command eye(m,n) produces $m$ by $n$ matrix with ones on the main diagonal. eye(n) produces $n×n$ matrix with ones on the main diagonal. The commands zeros, ones, randwork in a similar way.
>> v=[5; 6; 7]; C=[C v]
C =
1     0     0     5
0     1     0     6
0     0     1     7

The command C=[C v] appends the column vector v to the columns of $C$.
>> C=eye(3); v=[5 6 7]; C=[C; v]
C =
1     0     0
0     1     0
0     0     1
5     6     7

The command C=[C; v] appends the row vector v to the rows of $C$.
>> A=[0 1 2 3; 4 5 6 7; 8 9 -1 -2]
A =
0     1     2     3
4     5     6     7
8     9    -1    -2

>> A(:,2)= []
A =
0     2     3
4     6     7
8    -1    -2

>> A=[0 1 2 3; 4 5 6 7; 8 9 -1 -2]
A =
0     1     2     3
4     5     6     7
8     9    -1    -2

>> A(2:3,:)=[]
A =

0     1     2     3

The command A(:,2)=[] deletes the second column of $A$.
The command A(2:3,:)=[] deletes the second and third row of $A$.

7.2. Matrix Operations.

In Matlab, as in Pascal, Fortran, or C, matrix product is just C=A*B, where $A$ is an $m×n$ matrix and $B$ is $n×k$ matrix, and the resulting matrix $C$ has the size $m×k$.

• A+B and A-B are valid if the matrices $A$ and $B$ have the same size (i.e., the same number of rows and the same number of columns).
• A*B is valid if $A$’s number of columns equals $B$’s number of rows. The command A^2, which equals A*A, makes sense if $A$ is a square matrix.
• X=A/B (right division) is valid if $A$ and $B$ have the same number of columns and it denotes the solution to the matrix equation $XA=B$. In particular, if $A$ and $B$ are square matrices, $A∕B\equiv A{B}^{-1}$, if ${B}^{-1}$ exists.
• X=A\B (left division) is valid if $A$ and $B$ have the same number of rows and it denotes the solution to the matrix equation $AX=B$. In particular, if $A$ and $B$ are square matrices, $A\setminus B\equiv {A}^{-1}B$, if ${A}^{-1}$ exists.

If $A$ is $m×p$ and $B$ is $p×n$, their product $C=AB$ is $m×n$. Matlab defines their product by C = A*B, though this product can be also defined using Matlab for loops (see Section 9) and colon notation considered in the Section 7.1. However, the internally defined vectorized form of the product A*B is more efficient; in general, such vectorizations are strongly recommended, whenever possible.

>> A=pascal(3); B=magic(3);

>> C = A*B
C =
15    15    15
26    38    26
41    70    39

>> A=pascal(3); B=magic(3);
m=3; n=3;
for i = 1:m
for j = 1:n
C(i,j) = A(i,:)*B(:,j);
end
end

>> C
C =
15    15    15
26    38    26
41    70    39

The matrix A=magic(3) belongs to the family of so called magic squares. M = magic(n) returns an $n×n$ matrix ($n\ge 3$) constructed from the integers 1 through ${n}^{2}$ such that the $n$ numbers in all rows, all columns, and both diagonals sum to the same constant. The constant sum in every row, column and diagonal is called the magic constant or magic sum $M$. The magic constants of a magic square depend only on $n$ and have the values

$15,\phantom{\rule{0.3em}{0ex}}34,\phantom{\rule{0.3em}{0ex}}65,\phantom{\rule{0.3em}{0ex}}111,\phantom{\rule{0.3em}{0ex}}175,\phantom{\rule{0.3em}{0ex}}260,\phantom{\rule{0.3em}{0ex}}\dots ,\phantom{\rule{0.3em}{0ex}}\frac{n\left({n}^{2}+1\right)}{2},\phantom{\rule{0.3em}{0ex}}\dots$

(See, for example, http://en.wikipedia.org/wiki/Magic_square .)

>> magic(3)
ans =
8     1     6
3     5     7
4     9     2

>> magic(4)
ans =
16     2     3    13
5    11    10     8
9     7     6    12
4    14    15     1

>> magic(5)
ans =
17    24     1     8    15
23     5     7    14    16
4     6    13    20    22
10    12    19    21     3
11    18    25     2     9

The commands below check that the sums in all rows, columns, and in main diagonals of A=magic(5) are indeed the same and equal to 65.

>> A=magic(5);
>> for i=1:5
rows(i)=sum(A(i,:));
columns(i)=sum(A(:,i));
end

>> rows
rows =
65    65    65    65    65

>> columns
columns =
65    65    65    65    65

>> A=magic(5); dr=0;
for i=1:5
dr=dr+A(i,5-i+1);       % right diagonal
end

>> A=magic(5); dl=0;
for i=1:5
dl=dl+A(i,i);           % left diagonal
end

>> dr
dr =
65

>> dl
dl =
65

For information on pascal(n) matrices see Section 12.1.

The scalar (dot) product of two vectors or multiplication of a matrix with a vector can be done easily in Matlab.

>> a=[1 2 3 4]; b=[5 6 7 8];

>> dot(a,b)  % Scalar (dot) product of a and b
ans =
70

>> a=[1 2 3 4]; b=[5 6 7 8];

>> a*b’  % dot product using matrix product
ans =
70

>> A=[1 2 3; 4 5 6; 7 8 9]
A =
1     2     3
4     5     6
7     8     9

>> X=[0;-1;-2]
X =
0
-1
-2

>> A*X
ans =
-8
-17
-26

7.3. Array Operations.

Array operations are done on an element-by-element basis for matrices/vectors of the same size:

• .* denotes element-by-element multiplication;
• ./ denotes element-by-element left division;
• .\ denotes element-by-element right division;
• .^ denotes element-by-element exponentiation.
>> x=[1 2]
x =
1     2

>> y=[3 4]
y =
3     4

>> x.*y
ans =
3     8

>> x.^y
ans =
1    16

>> x./y
ans =
0.3333    0.5000

>> x.\y
ans =
3     2

Comparison of matrix and array operations

>> A=[1 2; 3 4]
A =
1     2
3     4

>> B=[5 6; 7 8]
B =
5     6
7     8

>> A*B
ans =
19    22
43    50

>> A.*B
ans =
5    12
21    32

>> A=[1 2; 3 4]
A =
1     2
3     4

>> B=[5 6; 7 8]
B =
5     6
7     8

>> A/B
ans =
3.0000   -2.0000
2.0000   -1.0000

>> A\B
ans =
-3    -4
4     5

>> A./B
ans =
0.2000    0.3333
0.4286    0.5000

>> A.\B
ans =
5.0000    3.0000
2.3333    2.0000

>> A^2
ans =
7    10
15    22

>> A.^2
ans =
1     4
9    16

7.4. Matrix Functions.

• expm(A) computes the exponential of a square matrix $A$, i.e., ${e}^{A}$. This operation should be distinguished from exp(A) which computes the element-by-element exponential of $A$;
• logm(A) computes $log\left(A\right)$ of square matrix $A$. $log\left(A\right)$ is the matrix such that $A={e}^{log\left(A\right)}$;
• sqrtm(A) computes $\sqrt{A}$ of square matrix. $\sqrt{A}$ is the matrix such that $A=\left(\sqrt{A}\right)\ast \left(\sqrt{A}\right)$.

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

>> expm(A)
ans =

2.7183    1.7183    1.0862
0    1.0000    1.2642
0         0    0.3679

>> exp(A)
ans =

2.7183    2.7183    1.0000
1.0000    1.0000    7.3891
1.0000    1.0000    0.3679

Notice that the diagonal elements of the two results are equal. This would be true for any triangular matrix. But the off-diagonal elements, including those below the diagonal, are different.

For additional help type in Matlab: help matfun
or find expm, logm, or sqrtm in:
http://www.mathworks.com/access/helpdesk/help/techdoc/index.html

### 8. Character Strings

(This section can be skipped in the first reading of the tutorial.)

S = ’Any Characters’ creates a character array, or string. All character strings are entered using single right-quote characters. The string is actually a vector whose components are the numeric codes for the characters (the first 127 codes are ASCII). The length of S is the number of characters. A quotation within the string is indicated by two single right-quotes.

S = [S1 S2 ...] concatenates character arrays S1, S2, etc. into a new character array, S.

S = strcat(S1, S2, ...) concatenates S1, S2, etc., which can be character arrays or Cell Arrays of Strings.

Trailing spaces in strcat character array inputs are ignored and do not appear in the output. This is not true for strcat inputs that are cell arrays of strings. Use the S = [S1 S2 ...] concatenation syntax, to preserve trailing spaces.

>> name = [’Thomas’ ’ R. ’ ’Lee’]
name =
Thomas R. Lee

>> name = strcat(’Thomas’,’ R.’,’ Lee’)
name =
Thomas R. Lee

>> S = [’toge’ ’ther’]
S =
together

>> S = strcat(’toge’, ’ther’)
S =
together

>> S = [’toge ’ ’ther’]
S =
toge ther

S = strcat(’toge ’, ’ther’)
S =
together

A column vector of text objects is created with C=[S1; S2; ...], where each text string must have the same number of characters.

>> S1=[’alpha ’;’gamma ’]
S1 =
alpha
gamma

>> S=[S1 S2]
S =
alpha ALPHA
gamma GAMMA

>> size(S)
ans =
2    11

>> S2=[’ALPHA’; ’GAMMA’]
S2 =
ALPHA
GAMMA

>> S=strcat(S1,S2)
S =
alphaALPHA
gammaGAMMA

>> size(S)
ans =
2    10

The command char(S1,S2,...) converts strings to a matrix. It puts each string argument S1, S1, etc., in a row and creates a string matrix by padding each row with the appropriate number of blanks; in other words, the command char does not require the strings S1, S2, …to have the same number of characters. One can also use strvcat and str2mat to create strings arrays from strings arguments, however, char does not ignore empty strings but strvcat does.

>> S=char(’Hello’,’’,’World’)
S =
Hello

World

>> size(S)
ans =
3     5

>> S(1,:)
ans =
Hello

>> S(2,:)
ans =

>> S(3,:)
ans =
World

>> S=strvcat(’Hello’,’’,’World’)
S =
Hello
World

>> size(S)
ans =
2     5

>> S(1,:)
ans =
Hello

>> S(2,:)
ans =
World

>> S(3,:)
??? Index exceeds matrix dimensions.

For additional help type in Matlab: help strings or help strfun
or find strings in:
http://www.mathworks.com/access/helpdesk/help/techdoc/index.html

8.1. eval Function.

eval executes string containing Matlab expression. For example, the command

>> eval(’y=exp(3)*cos(pi/4)’)
y =
14.2026

evaluates $y=exp\left(3\right)cos\left(\pi ∕4\right)$ and is equivalent to typing $y=exp\left(3\right)cos\left(\pi ∕4\right)$ on the command prompt. eval function is often used to load or create sequentially numbered data files. Here is a hypothetical example:

>> for d=1:10
s = [’load August’ int2str(d) ’.mat’]
eval(s)
end

Loads MAT-files August1.mat to August10.mat into the Matlab workspace.

These are the strings being evaluated:

For additional help type in Matlab: help eval
or find eval in:
http://www.mathworks.com/access/helpdesk/help/techdoc/index.html

### 9. for and while Loops and if Conditionals

for x=initval:endval, statements end command repeatedly executes one or more Matlab statements in a loop. Loop counter variable x is initialized to value initval at the start of the first pass through the loop, and automatically increments by 1 each time through the loop. The program makes repeated passes through statements until either x has incremented to the value endval, or Matlab encounters a break, or return instruction,

while expression, statements, end command repeatedly executes one or more Matlab statements in a loop, continuing until expression no longer holds true or until Matlab encounters a break, or return instruction.

if conditionally executes statements. The general form of the if command is

if expression
statements
elseif expression
statements
else
statements
end

The else and elseif parts are optional.

Example 1.

The graph of function $f\left(x\right)=3{x}^{2}-{e}^{x}$

indicates that $f\left(x\right)=3{x}^{2}-{e}^{x}$ has two positive zeros.

The computation of these two positive roots can be achieved by the so called fixed-point iteration scheme. For example, the smaller positive root $p$ is the fixed point of ${g}_{1}\left(x\right)=\sqrt{\left(\frac{1}{3}\cdot {e}^{x}\right)}$. defined on $\left[0,1\right]$, i.e., such that

${g}_{1}\left(p\right)=p,\phantom{\rule{2em}{0ex}}0\le p\le 1.$

The larger positive root is the fixed point $c$ of ${g}_{2}\left(x\right)=ln\left(3{x}^{2}\right)$ defined on $\left[3,4\right]$, i.e., such that

${g}_{2}\left(c\right)=c,\phantom{\rule{2em}{0ex}}3\le c\le 4.$

These fixed points (and ultimately the positive roots of $f\left(x\right)=3{x}^{2}-{e}^{x}$) can be computed accurately by running the so called fixed-point iteration scheme in the form:

And here is the actual Matlab session for the above fixed-point iteration scheme based on ${g}_{1}$:

>> g1=inline(’((1/3)*exp(x))^(1/2)’);
>> p0=1.0; iter=1;
>> for n=1:200
p=g1(p0);
err=abs(p-p0);
if err>=0.00001
iter=iter+1;
p0=p;
else
break
end
end
>> p
p =
0.9100

Define ${g}_{1};$
Define initial ${p}_{0}$ and iteration counter;
Start for loop;
Computation of the next approximation ${p}_{i+1}={g}_{1}\left({p}_{i}\right)$;
Compute the difference between the two consecutive iterations;
If the difference is greater than ($0.00001$), advance iteration by one step;
change ${p}_{0}$ to $p$, the next approximation computed;

Terminates the execution when the error is less than $0.00001$.

The computed value of the root of $f\left(x\right)=3{x}^{2}-{e}^{x}$.

and here fixed-point iteration scheme based on ${g}_{2}$:

>> g2=inline(’log(3*x^2)’);
>> p0=4.0; iter=1;
>> for n=1:200
p=g2(p0);
err=abs(p-p0);
if err>=0.00001
iter=iter+1;
p0=p;
else
break
end
end
>> p
p =
3.7331

Define ${g}_{2};$
Define initial ${p}_{0}$ and iteration counter;
Start for loop;
Computation if the next approximation ${p}_{i+1}={g}_{1}\left({p}_{i}\right)$;
Compute the difference between the two consecutive iterations;
If the difference is greater than ($0.00001$), advance iteration by one step;
change ${p}_{0}$ to $p$, the next approximation computed;

Terminates the execution when the error is less than $0.00001$.

The computed value of the root of $f\left(x\right)=3{x}^{2}-{e}^{x}$.

The accurate value of the root of $f\left(x\right)=3{x}^{2}-{e}^{x}$ in $\left[0,1\right]$ (to within $1{0}^{-10}$) is $p=0.9100075725$.

The accurate value of the root of $f\left(x\right)=3{x}^{2}-{e}^{x}$ in $\left[3,4\right]$ (to within $1{0}^{-10}$) is $p=3.733079029$.

Note: It is interesting to know that the fixed-point iteration based on ${g}_{1}$ does not converge to the larger positive root of $f\left(x\right)=3{x}^{2}-{e}^{x}$. Similarly, the fixed-point iteration based on ${g}_{2}$ does not converge to the smaller positive root of $f\left(x\right)=3{x}^{2}-{e}^{x}$.

Example 2.

The Euler numerical method for approximating the solution to the differential equation,

${y}^{\prime }=f\left(x,y\right),\phantom{\rule{1em}{0ex}}y\left({x}_{0}\right)={y}_{0},$

has the form

In the above scheme, ${y}_{n}$, for $n=1,2,\dots \phantom{\rule{0.3em}{0ex}}$, approximates true solution at ${x}_{n}$, i.e., $y\left({x}_{n}\right)$.

The initial value problem

${y}^{\prime }+2y=3exp\left(-4x\right),\phantom{\rule{1em}{0ex}}y\left(0\right)=1,$

has the exact solution

$y\left(x\right)=2.5exp\left(-2x\right)-1.5exp\left(-4x\right).$

We will use the Euler method to approximate this solution. We will also compute the differences between the exact and approximate solutions.

>> h=0.1;
>> x=0:h:2;
>> clear y;             % clear possibly existing variable y
>> y(1)=1;              % x1=0, thus y(1) corresponds to y(x1)=1
>> f=inline(’3*exp(-4*x)-2*y’)
f =
Inline function:
f(x,y) = 3*exp(-4*x)-2*y
>> size(x)              % size of x to be used in determining the size of vector i
ans =
1    21
>> for i=1:20 y(i+1)=y(i)+h*f(x(i),y(i)); end
>> Y_exact=2.5*exp(-2*x)-1.5*exp(-4*x);
>> plot(x,y,’--’,x,Y_exact)

### 10. Interactive Input

(This section can be skipped in the first reading of the tutorial.)

The commands input, keyboard, menu, and pause provide interactive user input and can be used inside a script or function file.

>> R = input(’How many apples: ’)
How many apples: sqrt(36)+factorial(4)
R =
30

gives the user the prompt in the text string and then waits for input from the keyboard. The input can be any Matlab expression, which is evaluated, using the variables in the current workspace, and the result returned in R. If the user presses the return key without entering anything, input returns an empty matrix.
>> R = input(’How many apples: ’,’s’)
How many apples: sqrt(36)+factorial(4)
R =
sqrt(36)+factorial(4)

gives the prompt in the text string and waits for character string input. The typed input is not evaluated; the characters are simply returned as a Matlab string.

keyboard command, when placed in a script or function file, stops execution of the file and gives control to the user’s keyboard. The special status is indicated by a K appearing before the prompt. Variables may be examined or changed - all Matlab commands are valid. The keyboard mode is terminated by executing the command return and pressing the return key. Control returns to the invoking script/function file.

Command menu generates a menu of choices for user input. On most graphics terminals menu will display the menu-items as push buttons in a figure window, otherwise they will be given as a numbered list in the command window. For example,

creates a figure with buttons labeled Red, Blue and Green. The button clicked by the user is returned as choice (i.e. choice = 2 implies that the user selected Blue).

After input is accepted, the dialog box closes, returning the output in choice. You can use choice to control the color of a graph:

t = 0:.1:60;
s = sin(t);
color = [’r’,’g’,’b’]
plot(t,s,color(choice))

pause commands waits for user response. pause(n) pauses for n seconds before continuing.

pause causes a procedure to stop and wait for the user to strike any key before continuing.

pause off indicates that any subsequent pause or pause(n) commands should not actually pause. This allows normally interactive scripts to run unattended.

pause on (default setting) indicates that subsequent pause commands should pause.

pause query returns the current state of pause, either on or off.

### 11. Input/Output Commands

(This section can be skipped in the first reading of the tutorial.)

Matlab provides many C-language file I/O functions for reading and writing binary and text files. The following five I/O functions: fopen, fclose, fgetl, fscanf, and fprintf are frequently used.

• fid = fopen(filename) opens or creates the filename. The filename argument is a string enclosed in single quotation marks.
Output value fid is a scalar Matlab integer that you use as a file identifier for all subsequent low-level file input/output routines.
• fclose(fid) closes the specified file if it is open, returning 0 if successful and -1 if unsuccessful. fid is an integer file identifier associated with an open file.
fclose(’all’) closes all open files.
• tline = fgetl(fid) returns the next line of the file associated with the file identifier fid. The returned string tline does not include the newline characters with the text line. To obtain the newline characters, use fgets function.
• A = fscanf(fid, format) reads data from the text file specified by fid, converts it according to the specified format string, and returns it in matrix A. Argument fid is an integer file identifier obtained from fopen. format is a string specifying the format of the data to be read.
• count = fprintf(fid, format, A, ...) formats the data in the real part of matrix A (and in any additional matrix arguments) under control of the specified format string, and writes it to the file associated with file identifier fid. count returns a count of the number of bytes written.
Argument fid is an integer file identifier obtained from fopen. It can also be 1 for standard output (the screen) or 2 for standard error. Omitting fid causes output to appear on the screen. For example, fprintf(’A unit circle has circumference %g radians.\n’, 2*pi) displays a line on the screen:
A unit circle has circumference 6.283186 radians.

The example below reads and displays the first three lines from the file circle_f2.m (see, Section 5) and available from: http://www.csun.edu/~hcmth008/matlab_tutorial/circle_f2.m
It also tests for end-of-file and then closes circle_f2.m file.

>> fid = fopen(’circle_f2.m’);
tline = fgetl(fid); i=1;
while feof(fid)==0 & i<=5     % Alternatively, you can use ischar(tline)
% == and <= are relational operators
% EQUAL and LESS THAN OR EQUAL,
% respectively;
% & denotes the logical operator AND.
disp(tline)                   % disp(X) displays an array, without
% printing the array name. If X contains
% a text string, the string is
% displayed.
tline = fgetl(fid);
i=i+1;
end
is_end_of_file = feof(fid)    % feof(fid) command checks for end-of-file.
fclose(fid);                  % Closes circle_f2.m file.
% The result is displayed below:

function [A,x,y] = circle_f2(r);
% CIRCLE_F2 - Function draws circle x^2+r^2=r^2
% and output its area together with x and y coordinates.

is_end_of_file =
0

The value of variable is_end_of_file = 0 indicates that the end of circle_f2.m file was not reached.
Changing name of the file and replacing i<=3 by i<=n will make the above script read and display the first n lines. The entire file will be displayed if the number of lines in the file is less or equal to n.

The next script reads and displays every line from the file my_f.m (created in Section 5) and available from: http://www.csun.edu/~hcmth008/matlab_tutorial/my_f.m
It also tests for end-of-file and then closes the file my_f.m. No a priori information about the number of lines in the file is needed in the script.

>> fid = fopen(’my_f.m’);
tline = fgetl(fid);
while feof(fid)==0         % Alternatively, you can use ischar(tline) instead of feof(fid)==0.
disp(tline)
tline = fgetl(fid);
end
is_end_of_file = feof(fid) % Checks for end-of-file.
fclose(fid);               % Closes my_f.m file.
% The result is displayed below:

function f = my_f(x);
% MY_F evaluates function f = exp(0.3x)sin(30x)+x^10/(1+x^10)
f = exp(-0.3*x).*sin(30*x)+x.^10./(1+x.^10);

is_end_of_file =
1

The value of variable is_end_of_file = 1 indicates that the end of my_f.m file was reached.

The next example is more advanced. It shows how to read and display lines from given file in a non-sequential order. It uses eval and int2str commands for a customized display. The script below reads and displays the 4th, 5th, 8th, and 9th lines from circle_f2.m file. It also appends the numbers to the above lines. No a priori information about the number of lines in the file is needed in the script.

>> fid = fopen(’circle_f2.m’); % output -->
tline = fgetl(fid); i=1;       % output -->
while feof(fid)==0             % output -->
if i<=3                        % output -->
tline = fgetl(fid);
end
if 4<=i & i<=5
disp([’<’ eval(’int2str(i)’) ’>’ ’ ’ tline])
tline = fgetl(fid);
end
if i>=6 & i<=7
tline = fgetl(fid);
end
if 8<=i & i<=9
disp([’<’ eval(’int2str(i)’) ’>’ ’ ’ tline])
tline = fgetl(fid);
end
if i>=10
tline = fgetl(fid);
end
i=i+1;
end
fclose(fid);

<4> alpha = linspace(0,2*pi,100);
<5> x = r*cos(alpha);
<8> plot(x,y);
<9> axis(’equal’);

The next examples show how to use fprintf and fscanf commands. First, create a text file called exp.txt containing a short table of the exponential function. (On Windows platforms, it is recommended that you use fopen with the mode set to ’wt’ to open a text file for writing.)

x = 0:.1:1;
y = [x; exp(x)];
fid = fopen(’exp.txt’, ’wt’);
fprintf(fid, ’%6.2f %12.8f\n’, y);
fclose(fid)

%6.2f stands for fixed-point notation field of width 8 with two digits after the decimal point, while %12.8f stands for fixed-point notation field of width 12 with 8 digits after the decimal point. \n stands for new line.

Now, examine the contents of exp.txt with type filename command.

type exp.txt

0.00   1.00000000
0.10   1.10517092
0.20   1.22140276
0.30   1.34985881
0.40   1.49182470
0.50   1.64872127
0.60   1.82211880
0.70   2.01375271
0.80   2.22554093
0.90   2.45960311
1.00   2.71828183

One can now read this file back into a two-column Matlab matrix with the script:

fid = fopen(’exp.txt’, ’r’);
a = fscanf(fid, ’%f %f’, [2 inf]);
a = a’;
fclose(fid);

It has two rows now. inf indicates to read to the end of the file.

Comparing it with exp.txt and long format outputs we observe no loss of accuracy.

>> disp(a)
0    1.0000
0.1000    1.1052
0.2000    1.2214
0.3000    1.3499
0.4000    1.4918
0.5000    1.6487
0.6000    1.8221
0.7000    2.0138
0.8000    2.2255
0.9000    2.4596
1.0000    2.7183

>> type exp.txt
0.00   1.00000000
0.10   1.10517092
0.20   1.22140276
0.30   1.34985881
0.40   1.49182470
0.50   1.64872127
0.60   1.82211880
0.70   2.01375271
0.80   2.22554093
0.90   2.45960311
1.00   2.71828183

>> format long; disp(a)
0   1.000000000000000
0.100000000000000   1.105170920000000
0.200000000000000   1.221402760000000
0.300000000000000   1.349858810000000
0.400000000000000   1.491824700000000
0.500000000000000   1.648721270000000
0.600000000000000   1.822118800000000
0.700000000000000   2.013752710000000
0.800000000000000   2.225540930000000
0.900000000000000   2.459603110000000
1.000000000000000   2.718281830000000

File temp.dat is available from: http://www.csun.edu/~hcmth008/matlab_tutorial/temp.dat

Open the file using fopen and read it with fscanf. If you include ordinary characters (such as the degrees character ${}^{\circ }$ (ASCII 176) and Fahrenheit (F) symbols used here) in the conversion string, fscanf skips over those characters when reading the string. This allows for further processing (see below). The display on the right is given for a comparison.

>> fid = fopen(’temp.dat’, ’r’);
degrees = char(176);
a=fscanf(fid, [’%d’ degrees ’F’]);
fclose(fid);
disp(a’)
78    72    65   105   -10

>> sin(a)
ans =
0.5140
0.2538
0.8268
-0.9705
0.5440

>> fid = fopen(’temp.dat’, ’r’);
a=fscanf(fid, ’%c’);
%
fclose(fid);
disp(a)

>> sin(a)
??? Undefined function or method ’sin’ for input
arguments of type ’char’.

### 12. Linear Algebra

12.1. Solving Linear Square Systems.

Write the system of linear algebraic equations

$\begin{array}{llll}\hfill {x}_{1}+{x}_{2}+{x}_{3}+{x}_{4}& =1\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill {x}_{1}+2{x}_{2}+3{x}_{3}+4{x}_{4}& =-1\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill {x}_{1}+3{x}_{2}+6{x}_{3}+10{x}_{4}& =2\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill {x}_{1}+4{x}_{2}+10{x}_{3}+20{x}_{4}& =3\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$

in matrix form

$AX=b,$

with

$A=\left[\begin{array}{cccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 2\hfill & \hfill 3\hfill & \hfill 4\hfill \\ \hfill 1\hfill & \hfill 3\hfill & \hfill 6\hfill & \hfill 10\hfill \\ \hfill 1\hfill & \hfill 4\hfill & \hfill 10\hfill & \hfill 20\hfill \end{array}\right],\phantom{\rule{2em}{0ex}}X=\left[\begin{array}{c}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \\ \hfill {x}_{3}\hfill \\ \hfill {x}_{4}\hfill \end{array}\right],\phantom{\rule{2em}{0ex}}b=\left[\begin{array}{c}\hfill 1\hfill \\ \hfill -1\hfill \\ \hfill 2\hfill \\ \hfill 3\hfill \end{array}\right].$

Please note that the number of columns of $A$ equals to the number of rows of $X$, and thus the matrix multiplication of $A$ with $X$ is defined.

Enter the matrix $A$ and vector $b$, and solve for vector $X$ with A\b.

>> A = [1 1 1 1; 1 2 3 4; 1 3 6 10; 1 4 10 20]
A =
1     1     1     1
1     2     3     4
1     3     6    10
1     4    10    20

>> b = [1; -1; 2; 3]
b =
1
-1
2
3

>> X = A\b   % Solve for X
X =
15
-33
26
-7

>> c = A*X   % Check the solution
c =
1
-1
2
3

The matrix $A$ belongs to the family of so called Pascal matrices and can be entered by the command A = pascal(4).

>> pascal(1)
ans =
1

>> pascal(2)
ans =
1     1
1     2

>> pascal(3)
ans =
1     1     1
1     2     3
1     3     6

>> pascal(4)
ans =
1     1     1     1
1     2     3     4
1     3     6    10
1     4    10    20

In general, A = pascal(n) returns the Pascal matrix of order n: a symmetric positive definite matrix with integer entries taken from Pascal’s triangle that determines the coefficients which arise in binomial expansions. See, for example, http://en.wikipedia.org/wiki/Pascal’s_triangle .

12.1.1. Singular Coefficient Matrix. A square matrix $A$ is singular if it does not have linearly independent columns. $A$ is singular if and only if the determinant of $A$ is zero (see Section 12.3). If $A$ is singular, the equation $AX=b$ either has no solution or it has infinite many solutions. Actually, one can show that $AX=b$ system has a solution if and only if $b\perp N$ (i.e, the scalar product of $b$ and $N$ is zero) for all $N$ such that ${A}^{\prime }N=0$, where ${A}^{\prime }$ is the transpose (the complex conjugate transpose, if $A$ has complex entries) of $A$, i.e., the matrix where the original rows and columns are interchanged. The matrix

>> A=[1 3 7; -1 4 4; 1 10 18]
A =
1     3     7
-1  \label{sub:factors}\label{sub:factors}   4     4
1    10    18

is singular, as easy verified by typing

>> det(A)       % det(A) calculates the determinant of A
ans =
0

For b = [5;2;12], the equation $AX=b$ has a solution given by

>> b = [5;2;12]; format rat; X=pinv(A)*b, format;

X =          % "format rat" requests a rational format,
82/213   % while the last command "format" returns
-47/426   % to the default format.
301/426

>> A*X

ans =    % Solution’s verification
5.0000
2.0000
12.0000

For more information on pinv(A), see Section 12.4. In this case, command A\b returns an error condition.

One can also verify that b = [5;2;12] is perpendicular to any solution of the homogeneous system ${A}^{\prime }N=0$:

>> A’   % transpose of A

ans =
1    -1     1
3     4    10
7     4    18

>> N = null(A’,’r’)

N =               % Solves the homogeneous equation A’*N=0
-2             % with an option ’r’ requesting
-1             % a rational format
1

and

>> dot(N,b)
ans =
0

This verifies that $b\perp N$.

The general (complete) solution to $AY=b$ is given by

where $Z$ solves the homogeneous system $AZ=0$ and $c$ is any constant (scalar), as verified by

>> format rat; Z = null(A,’r’), format;

Z =                   % Solves the homogeneous equation
16/7               % A*Z=0 with an option ’r’
-11/7               % requesting a rational format.
1

>> c = 1.5; A*(c*Z+X)

ans =
5.0000
2.0000
12.0000

On the other hand, for b = [1;2;3], $b\perp ̸N$. (indeed, dot(N,b)=-0.4082) and the system $AX=b$ does not have an exact solution.

In this case, pinv(A)*b returns a least squares solution $X$, i.e.,

$X=\underset{Y}{min}\parallel AY-b\parallel ,$

where for $v=\left({v}_{1},{v}_{2},\dots ,{v}_{n}\right)$, $\parallel v\parallel ={\left(\sum _{i=1}^{n}|{v}_{i}{|}^{2}\right)}^{1∕2}$ is the so called norm of vector $v$. Matlab command for the norm of vector v is norm(v). We have

>> b = [1;2;3]

b =
1
2
3

>> format rat; X=pinv(A)*b, format;

X =                   % a least squares solution
-395/1278
1081/2556
-107/2556

>> format rat; A*X, format;

ans =          % It does not get back the original vector
2/3       % b=[1;2;3].
11/6
19/6

>> norm(A*X-b)

ans =
0.4082

12.2. Gaussian Elimination.

One of the techniques learned in linear algebra courses for solving a system of linear equations is Gaussian elimination. One forms the augmented matrix that contains both the coefficient matrix $A$ and vector $b$, and then, using Gauss-Jordan reduction procedure, the augmented matrix is transformed into the so called row reduced echelon form. Matlab has built-in function that does precisely this reduction. For A = pascal(4) and b = [1; -1; 2; 3], considered in Section 12.1, the commands C=[A b] and rref(C) calculate the augmented matrix and the row reduced echelon matrix, respectively.

>> A = pascal(4)
A =
1     1     1     1
1     2     3     4
1     3     6    10
1     4    10    20

>> b = [1; -1; 2; 3]
b =
1
-1
2
3

>> C = [A b]; rref(C)
ans =
1     0     0     0    15
0     1     0     0   -33
0     0     1     0    26
0     0     0     1    -7

Since matrix A=pascal(4) is not singular, the last column of rref(C) (One could also use rref([A b]) command here.) provides the solution X already obtained in Section 12.1 by using X = A\b command.

Here is another use of rref command. One can determine whether $AX=b$ has an exact solution by finding the row reduced echelon form of the augmented matrix [A b]. For example, enter

format rat; A = [1 3 7; -1 4 4; 1 10 18]; b = [1;2;3]; rref([A b]), format;

ans =                             % For better visibility, the command
1     0     16/7    0      % "format rat" displays the result in
% rational format, while the last
0     1     11/7    0      % command  "format" returns the
% default format.
0     0     0       1

Since the bottom row contains all zeros except for the last entry, the equation $AX=b$ does not have a solution. Indeed, the above row reduced echelon form is equivalent to the following system of equations:

$\begin{array}{lll}\hfill 1\cdot {x}_{1}+0\cdot {x}_{2}+\left(16∕7\right){x}_{3}=0& \phantom{\rule{2em}{0ex}}& \hfill \\ \hfill 0\cdot {x}_{1}+1\cdot {x}_{2}+\left(11∕7\right){x}_{3}=0& \phantom{\rule{2em}{0ex}}& \hfill \\ \hfill 0\cdot {x}_{1}+0\cdot {x}_{2}+0\cdot {x}_{3}=1& \phantom{\rule{2em}{0ex}}& \hfill \end{array}$

As indicated in Section 12.4, in this case, pinv(A)*b command returns a least squares solution.

12.3. Inverses, Determinants, and Cramer’s Rule.

For a square and nonsingular matrix $A$, the equations $AX=I$ and $XA=I$ (here $I$ is an identity matrix), have the same solution, $X$. This solution is called the inverse of $A$ and is denoted by ${A}^{-1}$. Matlab function inv(A) computes the inverse of $A$.

The inverse of A=pascal(4) is

>> A = pascal(4)
A =
1     1     1     1
1     2     3     4
1     3     6    10
1     4    10    20

>> format rat; inv(A)
ans =
4    -6     4    -1
-6    14   -11     3
4    11    10    -3
-1     3    -3     1

and thus, the solution of the system $AX=b$ with b=[1;0;2;0] is

>> b = [1;0;2;0]; format rat; X = inv(A)*b

X =
12
-28
24
-7

>> Y=A\b    % the backslash operator is another
% way to compute the same solution
% and is much more efficient for
% large size matrices
Y =
12
-28
24
-7

The determinant is a special number associated with any square matrix $A$ and is denoted by $|A|$ or $det\left(A\right)$. The meaning of a determinant is a scale factor for measure when the matrix is regarded as a linear transformation. Thus a $2×2$ matrix, with determinant 3, when applied to a set of points with finite area will transform those points into a set with 3 times the area.

The determinant of

and the Rule of Sarrus is a mnemonic for this formula:

The determinant of $n×n$ matrix $A={\left({a}_{ij}\right)}_{i,j=1,2,\dots n}$ ($n\ge 1$) can be defined by the Leibniz formula:

Here the sum is computed over all permutations $\sigma$ of the numbers $\left\{1,2,\dots ,n\right\}$. A permutation is a function that reorders this set of integers. The set of such permutations is denoted by ${S}_{n}$. The sign function of permutations in the permutation group ${S}_{n}$ returns $+1$ and $-1$ for even and odd permutations, respectively. For example, for $n=4$ and $\sigma =\left(1,4,3,2\right)$, $\text{sign}\left(\sigma \right)=-1$ (one pair switch) and $\prod _{i=1}^{4}{a}_{i\phantom{\rule{0.3em}{0ex}}\sigma \left(i\right)}={a}_{11}{a}_{24}{a}_{33}{a}_{42}$.

Furthermore, for $A=\left[\begin{array}{cc}\hfill {a}_{11}\hfill & \hfill {a}_{12}\hfill \\ \hfill {a}_{21}\hfill & \hfill {a}_{22}\hfill \end{array}\right]$, ${S}_{2}=\left\{{\sigma }_{1},{\sigma }_{2}\right\}$, with ${\sigma }_{1}=\left(1,2\right)$, ${\sigma }_{2}=\left(2,1\right)$, $\text{sign}\left({\sigma }_{1}\right)=+1$, $\text{sign}\left({\sigma }_{2}\right)=-1$, and we have

Matlab function det(A) computes the determinant of a square matrix A. For A = pascal(4), we have

>> A = pascal(4)
A =
1     1     1     1
1     2     3     4
1     3     6    10
1     4    10    20

>> det(A)
ans =
1

Actually, det(pascal(n)) is equal to 1 for any $n\ge 1$.

A square matrix $A$ has the inverse ${A}^{-1}$ if and only if $det\left(A\right)\ne 0$. In this case, $det\left({A}^{-1}\right)=1∕det\left(A\right)$.

12.3.1. Cramer’s rule. For square matrices $A$ that are invertible, the system of equations $AX=b$ with a given $n$-dimensional column vector $b$ can be solved by using so called Cramer’s rule. This rule states that the solution $X={\left({x}_{1},{x}_{2},\dots ,{x}_{n}\right)}^{T}$ is obtained by the formula:

${x}_{i}=\frac{det\left({A}_{i}\right)}{det\left(A\right)},\phantom{\rule{1em}{0ex}}i=1,2,\dots ,n,$

where ${A}_{i}$ is the matrix formed by replacing the ith column of $A$ by the column vector $b$.

Here is Matlab code for Cramer’s rule applied to A=magic(5) and b=[1;2;3;4;5]

>> A = magic(5);
>> b = [1;2;3;4;5];
>> format rat;
>> for i=1:5                          % only these three
C=A; C(:,i)=b; X(i,:)=det(C)/det(A);  % lines represent
end                                % the actual code

>> format rat; X
X =
1/78
1/78
7/39
1/78
1/78

Check that the backslash operator A\b and the command inv(A)*b provide the same solution.

Note: For large size matrices, the use of Cramer’s rule and of the inverse is not computationally efficient as compared to the backslash operator.

12.4. Pseudoinverses or Generalized Inverses.

(This section can be skipped in the first reading of the tutorial.)

The pseudoinverse ${A}^{+}$ f an $m×n$ matrix $A$ is a generalization of the inverse matrix. The terms Moore-Penrose pseudoinverse or generalized inverse are also used as synonyms for pseudoinverse. Pseudoinverse is a unique $n×m$ matrix ${A}^{+}$ satisfying all of the following four conditions:

1. $A{A}^{+}A=A$;
2. ${A}^{+}A{A}^{+}={A}^{+}$;
3. ${\left(A{A}^{+}\right)}^{\prime }=A{A}^{+}$ (i.e.,$A{A}^{+}$ is Hermitian);
4. ${\left({A}^{+}A\right)}^{\prime }={A}^{+}A$ (i.e.,${A}^{+}A$ is Hermitian).

Here ${A}^{\prime }$ is the Hermitian transpose (also called conjugate transpose) of a matrix $A$. For matrices whose elements are real numbers instead of complex numbers, ${A}^{\prime }={A}^{T}$, In Matlab notation, Hermitian transpose of A is A’.

Pseudoinverses satisfy many properties. Here are some of them:

• The pseudoinverse exists and is unique for any matrix $A$.
• If $A$ has real entries, then so does ${A}^{+}$; if $A$ has rational entries, then so does ${A}^{+}$.
• If the matrix $A$ is invertible, the pseudoinverse and the inverse coincide, i.e., ${A}^{+}={A}^{-1}$.
• The pseudoinverse of a zero matrix is its transpose.
• The pseudoinverse of the pseudoinverse is the original matrix, i.e., ${\left({A}^{+}\right)}^{+}=A$.
• Pseudoinversion commutes with transposition, conjugation, and taking the conjugate transpose.
• The pseudoinverse of a scalar multiple of $A$ is the reciprocal multiple of ${A}^{+}$, i.e., ${\left(\alpha A\right)}^{+}={\alpha }^{-1}{A}^{+}$ for $\alpha \ne 0$.

A common use of the pseudoinverse is to compute a best fit (least squares) solution to a system of linear equations that lacks a unique solution. In general, a given linear system $AX=b$ may not have a solution. Even if there exists such a solution, it may not be unique. One can however ask for a vector $X$ that brings $AX$ as close as possible to vector $b$, i.e., a vector that minimizes the Euclidean norm

$\parallel AX-b\parallel ,$

where for $v=\left({v}_{1},{v}_{2},\dots ,{v}_{n}\right)$, $\parallel v\parallel ={\left(\sum _{i=1}^{n}|{v}_{i}{|}^{2}\right)}^{1∕2}$. Matlab command for the norm of vector v is norm(v).

If there are several such vectors $X$, one could ask for the one among them with the smallest Euclidean norm. Thus formulated, the problem has a unique solution, given by the pseudoinverse:

$X={A}^{+}b.$

Matlab function that computes the pseudoinverse of a matrix A is pinv(A). If for given A and b, the least square solution X is unique, it can be computed by entering either X=pinv(A)*b or X=A\b. If there exist many least square solutions (for the so called rank-deficient matrices), X=pinv(A)*b computes vector X with the smallest norm norm(X), while X=A\b computes a vector X that has at most $r$ nonzero components. Here, $r$ is the rank of matrix A equal to rank(A).

Sections 12.5.1 and 12.5.2 provide many examples for the use of pinv(A).

12.5. Underdetermined and Overdetermined Linear Systems.

12.5.1. Underdetermined linear systems. Underdetermined linear systems involve more unknowns than equations. Such systems are studied in linear programming, especially when accompanied by additional constraints. A solution, when exists, is never unique. For A=[0 1 1 0; 1 1 0 1; 1 2 1 1] and b=[1;2;1], the corresponding system has three equations and four unknowns:

$\begin{array}{lll}\hfill {x}_{2}+{x}_{3}=1& \phantom{\rule{2em}{0ex}}& \hfill \\ \hfill {x}_{1}+{x}_{2}+{x}_{4}=2& \phantom{\rule{2em}{0ex}}& \hfill \\ \hfill {x}_{1}+2{x}_{2}+{x}_{3}+{x}_{4}=1& \phantom{\rule{2em}{0ex}}& \hfill \end{array}$
>> A=[0 1 1 0; 1 1 0 1; 1 2 1 1]; b=[1;2;1]; reduced_echelon_form =rref([A b])

reduced_echelon_form =
1     0    -1     1     0
0     1     1     0     0
0     0     0     0     1

The row reduced echelon form shows that the system has no exact solutions. In this case, least square solutions $X$ such that

can be obtained by entering A\b or pinv(A)*b (see Section 12.4):

>> format rat; X=A\b

Warning: Rank deficient, rank=2,  tol=2.1756e-15.
X =
1
1/3
0
0

>> format rat; Y=pinv(A)*b

Y =
7/15
2/5
-1/15
7/15

We observe that there are many vectors Z (here, X and Y) that minimize norm(A*Z-b). This is due to the fact that matrix $A$ has only 2 linearly independent columns (so called rank-deficient system, while the full rank would be 3), the least square problem is not unique.

Consider now the following two equations with four unknowns:

$\begin{array}{lll}\hfill 6{x}_{1}+8{x}_{2}+7{x}_{3}+3{x}_{4}=1& \phantom{\rule{2em}{0ex}}& \hfill \\ \hfill 3{x}_{1}+5{x}_{2}+4{x}_{3}+{x}_{4}=2& \phantom{\rule{2em}{0ex}}& \hfill \end{array}$

Enter A=[6 8 7 3;3 5 4 1] and b=[1;2], and calculate rref([A b]):

>> A=[6 8 7 3;3 5 4 1]; b=[1;2]; format rat; rref([A b])

ans =

1    0   1/2   7/6  -11/6
0    1   1/2  -1/2    3/2

Thus, the system has an exact solution that can be obtained by entering X=A\b or W=pinv(A)*b:

>> format rat; X=A\b

X =
0
5/7
0
-11/7

>> format rat; W=pinv(A)*b

W =
-81/137
119/137
19/137
-154/137

Both X and W are exact solutions:

>> A*X

ans =
1
2

>> A*W

ans =
1
2

The general (complete) solution to $AY=b$ is given by

where ${c}_{1}$ and ${c}_{2}$ are any constants (scalars), and ${Z}_{1}$ and ${Z}_{2}$ are fundamental solutions of the system $AZ=0$, as verified by

>>D = null(A,’r’)

D =
-1/2           -7/6
-1/2            1/2
1              0
0              1

>> Z1 = D(:,1); A*Z1

ans =
0
0

>> Z2 = D(:,2); A*Z2

ans =
0
0

The other solution, W=pinv(A)*b, is also of the form ${c}_{1}{Z}_{1}+{c}_{2}{Z}_{2}+X$, for some constants ${c}_{1}$ and ${c}_{2}$. Indeed, the constants ${c}_{1}$ and ${c}_{2}$ are the solutions of the following linear system D*C=W-X with C=[c1;c2], as verified by

>> D=null(A,’r’)

D =
-1/2           -7/6
-1/2            1/2
1              0
0              1

>> format rat; C=D\(W-X)

C =
19/137         %c1 = 19/137
429/959         %c2 = 429/959

12.5.2. Overdetermined linear systems. Overdetermined linear systems involve more equations than unknowns. For example, they are encountered in various kinds of curve fitting to experimental data (see Section 13). Enter the following data into Matlab:

t=[0;0.25;0.5;;0.75;1.00;1.25;1.50;1.75;2.00;2.25;2.50];
y=sin(t)+rand(size(t))/10;

We will model the above randomly corrupted sine function data with a quadratic polynomial function.

$y\left(t\right)={c}_{1}+{c}_{2}t+{c}_{3}{t}^{2}$

In other words, we want to approximate the vector y by a linear combination of three other vectors:

$\left[\begin{array}{c}\hfill y\left(0.00\right)\hfill \\ \hfill y\left(0.25\right)\hfill \\ \hfill y\left(0.50\right)\hfill \\ \hfill y\left(0.75\right)\hfill \\ \hfill y\left(1.00\right)\hfill \\ \hfill y\left(1.25\right)\hfill \\ \hfill y\left(1.50\right)\hfill \\ \hfill y\left(1.75\right)\hfill \\ \hfill y\left(2.00\right)\hfill \\ \hfill y\left(2.25\right)\hfill \\ \hfill y\left(2.50\right)\hfill \end{array}\right]={c}_{1}\left[\begin{array}{c}\hfill 1.00\hfill \\ \hfill 1.00\hfill \\ \hfill 1.00\hfill \\ \hfill 1.00\hfill \\ \hfill 1.00\hfill \\ \hfill 1.00\hfill \\ \hfill 1.00\hfill \\ \hfill 1.00\hfill \\ \hfill 1.00\hfill \\ \hfill 1.00\hfill \\ \hfill 1.00\hfill \end{array}\right]+{c}_{2}\left[\begin{array}{c}\hfill 0.00\hfill \\ \hfill 0.25\hfill \\ \hfill 0.50\hfill \\ \hfill 0.75\hfill \\ \hfill 1.00\hfill \\ \hfill 1.25\hfill \\ \hfill 1.50\hfill \\ \hfill 1.75\hfill \\ \hfill 2.00\hfill \\ \hfill 2.25\hfill \\ \hfill 2.50\hfill \end{array}\right]+{c}_{3}\left[\begin{array}{c}\hfill {\left(0.00\right)}^{2}\hfill \\ \hfill {\left(0.25\right)}^{2}\hfill \\ \hfill {\left(0.50\right)}^{2}\hfill \\ \hfill {\left(0.75\right)}^{2}\hfill \\ \hfill {\left(1.00\right)}^{2}\hfill \\ \hfill {\left(1.25\right)}^{2}\hfill \\ \hfill {\left(1.50\right)}^{2}\hfill \\ \hfill {\left(1.75\right)}^{2}\hfill \\ \hfill {\left(2.00\right)}^{2}\hfill \\ \hfill {\left(2.25\right)}^{2}\hfill \\ \hfill {\left(2.50\right)}^{2}\hfill \end{array}\right].$

The above system has eleven equations and three unknowns: ${c}_{1}$, ${c}_{2}$, and ${c}_{3}$. It is represented 11-by-3 matrix A and b=y. The row reduced echelon form of A shows the above system does not have an exact solution:

>> A = [ones(size(t)) t  t.^2]
A =
1.0000         0         0
1.0000    0.2500    0.0625
1.0000    0.5000    0.2500
1.0000    0.7500    0.5625
1.0000    1.0000    1.0000
1.0000    1.2500    1.5625
1.0000    1.5000    2.2500
1.0000    1.7500    3.0625
1.0000    2.0000    4.0000
1.0000    2.2500    5.0625
1.0000    2.5000    6.2500

>> b=y
b =
0.0863
0.4295
0.5158
0.7344
0.8706
0.9762
1.1714
1.0999
1.0193
0.8071
0.7691

>> rref([A b])
ans =
1     0     0     0
0     1     0     0
0     0     1     0
0     0     0     1
0     0     0     0
0     0     0     0
0     0     0     0
0     0     0     0
0     0     0     0
0     0     0     0
0     0     0     0

Since rank(A)=3, both the backslash operator (i.e., A\b) and pinv(A)*b provide the same (unique) least square solution:

>> X=A\b
X =
0.0207     % c1
1.2570     % c2
-0.4041     % c3

>> Y=pinv(A)*b
Y =
0.0207     % c1
1.2570     % c2
-0.4041     % c3

In other words, the least squares fit to the data is

$y\left(t\right)=0.0207+1.2570t-0.4041{t}^{2}$

The following commands evaluate the model at regularly spaced increments in $t$ (finer than the original spacing) and then plot the result including the original data.

>> T=(0:0.1:3.0)’;
>> Y=[ones(size(T)) T  T.^2]*X;
>> plot(T,Y,’-’,t,y,’o’)

Another way of obtaining the same least square solution is to enter the command
>> polyfit(t,y,2)
ans =
-0.4041    1.2570    0.0207   % c3   c2   c1

For more information on polyfit, see Section 13.1 on curve fitting and interpolation.

Note: Due to randomness of y = sin(t) + rand(size(t))/10, your results will differ slightly from the above interpolating polynomial and the graph.

And here is another attempt at fitting randomly corrupted sine function data into

$y\left(t\right)={c}_{1}+{c}_{2}sin\left(t\right)+{c}_{2}cos\left(t\right)+{c}_{4}exp\left(-t\right).$

A = [ones(size(t)) cos(t) sin(t) exp(-t)]
A =
1.0000    1.0000         0    1.0000
1.0000    0.9689    0.2474    0.7788
1.0000    0.8776    0.4794    0.6065
1.0000    0.7317    0.6816    0.4724
1.0000    0.5403    0.8415    0.3679
1.0000    0.3153    0.9490    0.2865
1.0000    0.0707    0.9975    0.2231
1.0000   -0.1782    0.9840    0.1738
1.0000   -0.4161    0.9093    0.1353
1.0000   -0.6282    0.7781    0.1054
1.0000   -0.8011    0.5985    0.0821

>> b=y
b =
0.0409
0.3069
0.5056
0.7419
0.9126
0.9712
1.0092
1.0137
0.9412
0.8205
0.6493

>> X=A\b
X =
-0.0461    % c1
-0.0336    % c2
1.0582    % c3
0.1777    % c4

>> T=(0:0.1:3.0)’;
>> Y=[ones(size(T)) cos(T) sin(T) exp(-T)]*X;
>> plot(T,Y,’-’,t,y,’o’)

And here is another attempt at fitting randomly corrupted sine function data into

$y\left(t\right)={c}_{1}+{c}_{2}sin\left(t\right)+{c}_{2}cos\left(t\right)+{c}_{4}exp\left(-t\right).$

The next two examples consider the systems with more equations than unknowns that have exact solutions.

Consider the system generated by 5 equations with 4 unknowns.

>> A=magic(5); A = A(:,1:4)
A =
17    24     1     8
23     5     7    14
4     6    13    20
10    12    19    21
11    18    25     2

b=[1;0;-43/648;0;0]
b =
1.0000
0
-0.0664
0
0

By computing the row reduced echelon form of A, we immediately get a solution.

>> rref([A b])
ans =
1.0000         0         0         0   -0.0026
0    1.0000         0         0    0.0434
0         0    1.0000         0   -0.0305
0         0         0    1.0000    0.0040
0         0         0         0         0

>> X=A\b
X =
-0.0026
0.0434
-0.0305
0.0040

>> Y=pinv(A)*b
Y =
-0.0026
0.0434
-0.0305
0.0040

This solution is also confirmed by using the backslash operator, A\b, and pinv(A)*b. The solution is unique since the homogeneous system $AZ=0$ has only the trivial solution, as verified by the command

>> null(A)
ans =
Empty matrix: 4-by-0

Finally, consider the linear system with 8 equations and 6 unknowns as generated by

>> A = magic(8); A = A(:,1:6)
A =
64     2     3    61    60     6
9    55    54    12    13    51
17    47    46    20    21    43
40    26    27    37    36    30
32    34    35    29    28    38
41    23    22    44    45    19
49    15    14    52    53    11
8    58    59     5     4    62

>> b = 260*ones(8,1)
b =
260
260
260
260
260
260
260
260

The scale factor 260 is the 8-by-8 magic sum. With all eight columns, one solution to A*X = b would be a vector of all 1’s. With only six columns, the equations are still consistent, as verified by looking at the row reduced echelon form of A.

rref([A b])
ans =
1     0     0     1     1     0     4
0     1     0     3     4    -3     8
0     0     1    -3    -4     4    -4
0     0     0     0     0     0     0
0     0     0     0     0     0     0
0     0     0     0     0     0     0
0     0     0     0     0     0     0
0     0     0     0     0     0     0

So a solution exists, but it is not all 1’s. Since the matrix is rank deficient (rank(A)=3), there are infinitely many solutions. Two of them are

>> format rat; X=A\b
Warning: Rank deficient, rank=3, tol=1.8829e-13.

X =
4
5
0
0
0
-1

>> format rat; W=pinv(A)*b

W =
15/13
19/13
18/13
18/13
19/13
15/13

The general (complete) solution to $AY=b$ is given by $Y={c}_{1}{Z}_{1}+{c}_{2}{Z}_{2}+{c}_{3}{Z}_{3}+X$ with

where ${c}_{1}$, ${c}_{2}$, and ${c}_{3}$ are any constants (scalars), and ${Z}_{1}$, ${Z}_{2}$, and ${Z}_{3}$ are fundamental solutions of the system $AZ=0$, as verified by

>> format rat; null(A)

ans =
719/2961      -542/2847      1389/2738
373/2112     -1262/1851       475/4663
-1161/6842        93/152        159/326
-2694/3481      -136/1045      -199/1255
316/595        642/2003     -2226/6383
-75/10834      257/3674      -102/173

12.6. Finding Eigenvalues and Eigenvectors.

An eigenvalue and eigenvector of a square $n×n$ matrix $A$ are a scalar $\lambda$ and a nonzero (n-dimensional) column vector $v$ that satisfy

 $Av=\lambda v.$ (1)

Usually, one obtains the solution by solving for the $n$ eigenvalues from the characteristic equation

$det\left(A-\lambda I\right)=0,$

and then solving for the eigenvectors by substituting the corresponding eigenvalues in (1), one at a time. The eigenvectors are determined up to multiplicative scalars. In Matlab, $\lambda$ and $v$ are computed using [V,D]=eig(A) command. For A=magic(3), and requesting the rational format (with fornat rat) for better visibility, we have

>> A=magic(3)

A =
8   1   6
3   5   7
4   9   2

>> format rat; [V,D]=eig(A)

V =
-780/1351      -735/904       -452/1323
-780/1351      1121/2378     -1121/2378
-780/1351       452/1323       735/904

D =
15              0              0
0           4801/980          0
0              0          -4801/980

Matrix $V$ contains the eigenvectors of $A$ as its columns, while matrix $D$ contains the eigenvalues on its diagonal. By comparing

>> A*V(:,1)

ans =
-1351/156
-1351/156
-1351/156

>> 15*V(:,1)

ans =
-1351/156
-1351/156
-1351/156

>> A*V(:,2)

ans =
-3541/889
1351/585
3319/1983

>> (4801/980)*V(:,2)

ans =
-3541/889
1351/585
3319/1983

and similarly for the third column of $V$, we easily recover the eigenvalues and the corresponding eigenvectors of the matrix $A$.

Also by checking

>> format rat; det(V)
ans =
-1121/1189

we observe that ${V}^{-1}$ exists, with the property
>> format rat; V*D*inv(V)

ans =
8    1   6
3    5   7
4    9   2

>> A=magic(3)

A =
8   1   6
3   5   7
4   9   2

Thus, we have $A=VD{V}^{-1}$. Furthermore, ${V}^{-1}AV=D$. Because of this last property, we say that the matrix $A$ is diagonalizable.

The eigenvalues and eigenvectors can be complex valued as the next example shows:

>> A=[0 -6 -1;6 2 -16;-5 20 -10]

A =
0    -6     -1
6     2    -16
-5    20    -10

>> format rat; [V,D]=eig(A)

V =
-383/460     899/4489 - 387/2777i   899/4489 + 387/2777i
-1071/3014  -753/3568 - 519/805i   -753/3568 + 519/805i
-489/1151  -1736/2505             -1736/2505

D =

-2294/747            0                      0
0       -1841/747 + 4277/243i            0
0                 0            -1841/747 - 4277/243i

Not every $n×n$ matrix $A$ has $n$ linearly independent eigenvectors. For example, it is easy to show that matrix $A=\left[\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right]$hasa double eigenvalue $\lambda =0$ and only one eigenvector $v=\left[\begin{array}{c}\hfill 1\hfill \\ \hfill 0\hfill \end{array}\right]$.

One can prove that if $n×n$ matrix $A$ has $n$ distinct eigenvalues ${\lambda }_{1},\phantom{\rule{0.3em}{0ex}}{\lambda }_{2},\phantom{\rule{0.3em}{0ex}}\dots ,\phantom{\rule{0.3em}{0ex}}{\lambda }_{n}$, then $A$ has $n$ linearly independent eigenvectors: ${K}_{1}$, ${K}_{2}$, …, ${K}_{n}$. Furthermore, if $V$ denotes $n×n$ matrix whose columns are eigenvectors ${K}_{1}$, ${K}_{2}$, …, ${K}_{n}$, then

Because of the identity ${V}^{-1}AV=D$, we say that the matrix $A$ is diagonalizable. An important class of diagonalizable matrices are square symmetric matrices whose entries are real. A square matrix $A$ is symmetric if it is equal to its transpose, i.e., $A={A}^{\prime }$. In other words, if for $A={\left({a}_{ij}\right)}_{i,j=1,2,\dots ,n}$ , we have ${a}_{ij}={a}_{ji}$, for $i,j=1,2,\dots ,n$.

A matrix with repeated eigenvalues may still have a full set of linearly independent vectors, as the next example shows:

>> A=[1 -2 2;-2 1 -2;2 -2 1]

A =
1    -2     2
-2     1    -2
2    -2     1

>> format rat; [V,D]=eig(A)

V =
-247/398       1145/2158       780/1351
279/1870      1343/1673      -780/1351
1040/1351      1013/3722       780/1351

D =
-1              0              0
0             -1              0
0              0              5

Matrix $A$ has eigenvalues ${\lambda }_{1}={\lambda }_{2}=-1$ and ${\lambda }_{3}=5$.

A more readable expressions for the eigenvectors can be obtained by entering the command [P,J]=jordan(A) (see Section 12.8.5 for more information on Jordan decomposition):

>> format rat; [P,J]=jordan(A)

P =
2/3      1/3     -1
1/3     -1/3      0
-1/3      1/3      1

J =
-1       0       0
0       5       0
0       0      -1

The columns of $P$ represent the eigenvectors of $A$. Furthermore, after taking into account that eigenvectors are determined up to multiplicative scalars, here are the eigenvalues and the corresponding eigenvectors:

${\lambda }_{1}=-1,\phantom{\rule{1em}{0ex}}{K}_{1}=\left[\begin{array}{c}\hfill 2\hfill \\ \hfill 1\hfill \\ \hfill -1\hfill \end{array}\right];\phantom{\rule{2em}{0ex}}{\lambda }_{2}=5,\phantom{\rule{1em}{0ex}}{K}_{2}=\left[\begin{array}{c}\hfill 1\hfill \\ \hfill -1\hfill \\ \hfill 1\hfill \end{array}\right];\phantom{\rule{2em}{0ex}}{\lambda }_{2}=-1,\phantom{\rule{1em}{0ex}}{K}_{3}=\left[\begin{array}{c}\hfill -1\hfill \\ \hfill 0\hfill \\ \hfill 1\hfill \end{array}\right].$

As verified below, the first two columns of $V$ (i.e., the eigenvectors corresponding to the repeated eigenvalue $-1$) can be expressed as the linear combinations of the eigenvectors ${K}_{1}$ and ${K}_{3}$ :

>> format rat; [3*P(:,1) P(:,3)]\V(:,1)

ans =
279/1870     % c1
624/679      % c2

>> format rat; [3*P(:,1) P(:,3)]\V(:,2)

ans =
1343/1673     % c1
947/881      % c2

Consider the following non diagonalizable matrix:

>> A=[5 4 2 1;0 1 -1 -1;-1 -1 3 0;1 1 -1 2]

A =
5     4     2     1
0     1    -1    -1
-1    -1     3     0
1     1    -1     2

>> [V,D]=eig(A)

V =
0.5774    0.5774    0.5774    0.7071
0.0000    0.0000   -0.5774   -0.7071
-0.5774   -0.5774    0.0000    0.0000
0.5774    0.5774    0.5774   -0.0000

D =
4.0000         0         0         0
0    4.0000         0         0
0         0    2.0000         0
0         0         0    1.0000

We observe that the first two columns of matrix $V$ are identical, thus matrix $A$ does not have 4 linearly independent eigenvectors. Further checking shows that Matlab cannot compute ${V}^{-1}$ within its default error of tolerance:

>> det(V)

ans =
1.3084e-16

>> inv(V)
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.602469e-16.

ans =

1.0e+15 *

1.8014    1.8014    1.8014   -0.0000
-1.8014   -1.8014   -1.8014    0.0000
0.0000   -0.0000    0.0000    0.0000
-0.0000   -0.0000   -0.0000   -0.0000

12.7. Generalized Eigenvectors.

(This section can be skipped in the first reading of the tutorial.)

A non diagonalizable $n×n$ matrix does not have $n$ linearly independent eigenvectors. Generalized eigenvectors are introduced to form a complete basis of such matrix. A generalized eigenvector of matrix $A$ is a nonzero vector $w$, which has associated with it an eigenvalue $\lambda$ having algebraic multiplicity $m\ge 1$, satisfying

 ${\left(A-\lambda I\right)}^{m}w=0,$ (2)

where $I$ denotes $n×n$ identity matrix. For $m=1$, generalized eigenvectors become ordinary eigenvectors.

Assume $\lambda$ is an eigenvalue of multiplicity $m\ge 1$, and there is only one eigenvector ${v}_{1}$. The generalized eigenvectors ${v}_{1},{v}_{2},\dots ,{v}_{m}$ satisfy

and they are linearly independent. Furthermore, for any $k\ge 1$,

 ${\left(A-\lambda I\right)}^{k}{v}_{k}=0.$ (3)

In particular, (3) yields

${\left(A-\lambda I\right)}^{m}{v}_{k}=0,\phantom{\rule{2em}{0ex}}k=1,2,\dots ,m,$

and thus implying (2).

Matlab jordan(A) function can be used to find generalized eigenvectors of matrix $A$ (see also Section 12.8.5 for more information on Jordan decomposition). Consider the following example:

>> A=[5 4 2 1;0 1 -1 -1;-1 -1 3 0;1 1 -1 2]

A =
5     4     2     1
0     1    -1    -1
-1    -1     3     0
1     1    -1     2

>> [P,J]=jordan(A)

P =

1     1    -1     1
0     0     1    -1
-1     0     0     0
1     0     0     1

J =
4     1     0     0
0     4     0     0
0     0     1     0
0     0     0     2

The numbers 1, 2, and 4 are the eigenvalues of matrix $A$. Furthermore 4 is double eigenvalue. Let us denote the columns of $P$ by ${p}_{1}$, ${p}_{2}$, ${p}_{3}$, and ${p}_{4}$, then we have

Indeed,

>> (A-2*eye(4))*P(:,4)
ans =

0
0
0
0

>> (A-1*eye(4))*P(:,3)
ans =

0
0
0
0

>> (A-4*eye(4))*P(:,1)

ans =
0
0
0
0

>> (A-4*eye(4))*P(:,2)

ans =
1
0
-1
1

Thus ${p}_{1}$, ${p}_{3}$, and ${p}_{4}$ are ordinary eigenvectors corresponding to the eigenvalues 4, 1, and 2, respectively. Since $\left(A-4I\right){p}_{2}={p}_{1}$ and $\left(A-4I\right){p}_{1}=0$, we also have ${\left(A-4I\right)}^{2}{p}_{2}=0$ and ${p}_{2}$ is a generalized eigenvector corresponding to the eigenvalue 4.

The matrix in the next example has the eigenvalues 1 (with multiplicity 3), 2 (with multiplicity 2), and 3 with multiplicity 1).

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

A =
1     1     1     2    -1     2
0     1     2    -1    -1     2
0     0     1     2     1     2
0     0     0     2     2    -1
0     0     0     0     2     2
0     0     0     0     0     3

>> format rat; [P,J]=jordan(A)
P =

29/4     -58      116      -29     -145/2       -493/4
7/2    -174/7   2378/49     0       -29         -58
5      -116/7   435/49      0         0         -29/2
3      -58/7    -87/49      0         0           0
2       0       -29/7       0         0           0
1       0         0         0         0           0

J =

3        0         0         0         0           0
0        2         1         0         0           0
0        0         2         0         0           0
0        0         0         1         1           0
0        0         0         0         1           1
0        0         0         0         0           1

Let us denote the columns of $P$ by ${p}_{1}$, ${p}_{2}$, ${p}_{3}$, ${p}_{4}$, ${p}_{5}$, ${p}_{6}$. then we have,

Since $\left(A-2I\right){p}_{2}=0$ and $\left(A-2I\right){p}_{3}={p}_{2}$, we have ${\left(A-2I\right)}^{2}{p}_{3}=0$ and ${p}_{3}$ is a generalized eigenvector corresponding to the eigenvalue 2. Also, since $\left(A-I\right){p}_{4}=0$, $\left(A-I\right){p}_{5}={p}_{4}$, $\left(A-I\right){p}_{6}={p}_{5}$, we have ${\left(A-I\right)}^{3}{p}_{6}=0$ and ${\left(A-I\right)}^{3}{p}_{5}=0$, implying that ${p}_{5}$ and ${p}_{6}$ are generalized eigenvectors corresponding to the eigenvalue 1. On the other hand, ${p}_{1}$, ${p}_{2}$, and ${p}_{4}$ are the ordinary eigenvectors corresponding to the eigenvalues, 3, 2, and 1, respectively.

In the last example, the matrix also has the eigenvalues 1 (with multiplicity 3), 2 (with multiplicity 2), and 3 with multiplicity 1), however here, there are two linearly independent eigenvectors corresponding to the eigenvalue 1.

>> A=[1 0 1 2 3 4;0 1 5 8 -1 2; 0 0 1 7 4 2;0 0 0 2 6 9; 0 0 0 0 2 5; 0 0 0 0 0 3]

A =
1     0     1     2     3     4
0     1     5     8    -1     2
0     0     1     7     4     2
0     0     0     2     6     9
0     0     0     0     2     5
0     0     0     0     0     3

>> format rat; [P,J]=jordan(A)

P =
489/4      -2093/10  -29117/56    -18837/43   22271/56    991/1554
2093/4     -2093/2    -10465/4      -2093       2093         0
295/2         0       -2093/10    -14651/43   22637/404     0
39           0           0        -2093/43   -4216/117     0
5           0           0            0       2093/258     0
1           0           0            0          0         0

J =
3           0           0            0          0         0
0           1           1            0          0         0
0           0           1            0          0         0
0           0           0            2          1         0
0           0           0            0          2         0
0           0           0            0          0         1

Let us denote the columns of $P$ by ${p}_{1}$, ${p}_{2}$, ${p}_{3}$, ${p}_{4}$, ${p}_{5}$, ${p}_{6}$. We have,

Since $\left(A-I\right){p}_{2}=0$ and $\left(A-I\right){p}_{3}={p}_{2}$, we have ${\left(A-I\right)}^{2}{p}_{3}=0$, thus ${p}_{3}$ is a generalized eigenvector corresponding to the eigenvalue 1. Also, since $\left(A-2I\right){p}_{4}=0$ and $\left(A-2I\right){p}_{5}={p}_{4}$, thus we have ${\left(A-2I\right)}^{2}{p}_{5}=0$ and ${p}_{5}$ is a generalized eigenvector corresponding to the eigenvalue 1. On the other hand, ${p}_{1}$, ${p}_{2}$, ${p}_{4}$, and ${p}_{6}$, are the ordinary eigenvectors corresponding to the eigenvalues, 3, 1, 2, and 1, respectively.

12.8. Matrix Factorizations, Matrix Powers and Exponentials.

(This section can be skipped in the first reading of the tutorial.)

12.8.1. Cholesky factorization. The Cholesky factorization expresses a symmetric (or more generally Hermitian) and positive definite matrix $A$ as the product of an upper triangular matrix $R$ and its transpose ${R}^{\prime }$, i.e., $A={R}^{\prime }R$.

>> A=pascal(5)

A =
1     1     1     1     1
1     2     3     4     5
1     3     6    10    15
1     4    10    20    35
1     5    15    35    70

>> R=chol(A)

R =
1     1     1     1     1
0     1     2     3     4
0     0     1     3     6
0     0     0     1     4
0     0     0     0     1

>> R’*R                 % verification

ans =
1     1     1     1     1
1     2     3     4     5
1     3     6    10    15
1     4    10    20    35
1     5    15    35    70

For a symmetric (or more generally Hermitian) and positive definite matrix the Cholesky factorization exists. In this case, it is also unique in the sense that there exists only one upper triangular matrix $R$ with strictly positive diagonal entries such that $A={R}^{\prime }R$.

12.8.2. LU factorization. The LU factorization is a matrix decomposition which writes a matrix $A$ as the product of a lower triangular matrix $L$ and an upper triangular matrix $U$, i.e., $A=LU$. The product sometimes includes a permutation matrix as well. The LU decomposition is a modified form of Gaussian elimination. Not all square matrices have LU factorization, For example, $A=\left[\begin{array}{cc}\hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill \end{array}\right]$cannot be expressed at the product of triangular matrices. Even when such a factorization exists it is not unique (since $A=LU=\left(-L\right)\left(-U\right)$), unless one puts some restrictions on $L$ and $U$ matrices. Often, one such conditions is to require that the diagonal of $L$ (or $U$) consists of ones. Though, the matrix

 $A=\left[\begin{array}{cc}\hfill ϵ\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 1∕ϵ\hfill & \hfill 1\hfill \end{array}\right]\left[\begin{array}{cc}\hfill ϵ\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill -1∕ϵ\hfill \end{array}\right]$ (4)
can be expressed as the product of the lower and upper triangular matrices, however, when $ϵ$ is small, the factors are large and magnify errors, so even the permutations matrix is not necessary, it is desirable. For a given matrix $A$, Matlab command [L,U,P]=lu(A) returns an upper triangular matrix $U$, a lower triangular matrix $L$ with ones on its diagonal, and a permutation matrix $P$, such that $LU=PA$. (A permutation matrix is matrix of zeros and ones that has exactly one entry 1 in each row and column.) Furthermore, such decomposition always exist (even $A$ is not square matrix). On the other hand, the command [L,U]=lu(A) returns an upper triangular matrix $U$, and a permuted lower triangular matrix $L$ such that $LU=A$. LU decomposition is used in numerical analysis to solve systems of linear equations or calculate the determinant. Finally, when $A$ is symmetric (Hermitian) and positive definite matrix, one can choose $L$ in such a way that $L={U}^{\prime }$ and we get back to the Cholesky factorization.

For matrix in (4) with $ϵ=1∕1000$, we have

>> format rat; A=[0.001 1;1 0]
A =
1/1000         1
1              0

>> [L,U,P]=lu(A)
L =
1              0
1/1000         1

U =
1              0
0              1

P =
0              1
1              0

>> format rat; A=[0.001 1;1 0]
A =
1/1000         1
1              0

>> [L,U]=lu(A)
L =
1/1000         1
1              0

U =
1              0
0              1

And here is another example.

>> A=magic(5)
A =
17     24       1       8       15
23      5       7      14       16
4      6      13      20       22
10     12      19      21        3
11     18      25       2        9

>> format rat; [L,U,P]=lu(A)
L =

1          0              0              0              0
17/23        1              0              0              0
11/23      359/467          1              0              0
4/23      118/467       1199/2322         1              0
10/23     226/467       1679/2322        12/13           1

U =
23          5              7             14             16
0        467/23         -96/23         -54/23          73/23
0          0           6787/273      -1350/467       -510/467
0          0              0            845/43        2752/145
0          0              0              0           -200/9

P =
0          1              0              0              0
1          0              0              0              0
0          0              0              0              1
0          0              1              0              0
0          0              0              1              0

12.8.3. QR factorization. An orthogonal matrix $Q$ is a real matrix whose columns all have length one and are perpendicular, equivalently ${Q}^{\prime }Q=I$, where $I$ is the identity matrix. This also implies that ${Q}^{-1}={Q}^{\prime }$. The following is the simplest $2×2$ orthogonal matrix representing two-dimensional coordinate rotations:

$\left[\begin{array}{cc}\hfill cos\left(\theta \right)\hfill & \hfill sin\left(\theta \right)\hfill \\ \hfill -sin\left(\theta \right)\hfill & \hfill cos\left(\theta \right)\hfill \end{array}\right]$

For complex matrices, the corresponding term is unitary. Orthogonal (unitary) matrices preserve length, preserve angles, and do not magnify errors, and thus they are often used in numerical computations.

The qr function performs the orthogonal-triangular decomposition of a matrix. This factorization is useful for both square and rectangular matrices. It expresses the matrix as the product of an orthogonal (unitary) matrix and an upper triangular matrix.

Matlab command [Q,R] = qr(A) produces an upper triangular matrix $R$ of the same dimension as $A$ and a unitary matrix $Q$ so that $A=QR$. If $m×n$ is the size of $A$, then $Q$ is $m×m$ and $R$ is $m×n$. Since we also have $A=\left(-Q\right)\left(-R\right)$, with $-Q$ being orthogonal (unitary) and $-R$ being upper triangular, this factorization is not unique. If $A$ is square and invertible matrix, then this factorization is unique if, for example, we require that the diagonal elements of $R$ are positive (negative).

>> A=[12 -51 4; 6 167 -68;-4 24 -41]
A =
12            -51              4
6            167            -68
-4             24            -41

>> format rat; [Q,R]=qr(A)
Q =
-6/7           69/175         58/175
-3/7         -158/175         -6/175
2/7           -6/35          33/35

R =
-14            -21             14
0           -175             70
0              0            -35

>> Q*R            % verification
ans =
12            -51              4
6            167            -68
-4             24            -41

12.8.4. Schur decomposition. For a matrix $A$, there exist orthogonal (unitary) matrix $Q$ and an upper triangular matrix $U$ such that

$A=QU{Q}^{-1}.$

Note that ${Q}^{-1}={Q}^{\prime }$. Since $U$ is similar to $A$ and $U$ is triangular, the eigenvalues of $A$ are the diagonal entries of $U$. The corresponding Matlab command is [Q,U]=schur(A).

>> A=[6 12 19;-9 -20 -33;4 9 15]
A =
6             12             19
-9            -20            -33
4              9             15

>> format rat; [Q,U]=schur(A)
Q =
-540/1139       926/1393       780/1351
1467/1805       286/3657       780/1351
-489/1444     -2506/3373       780/1351

U =
-1           1351/65       -9520/213
0              1           -523/858
0              0              1

12.8.5. Jordan normal (canonical) form.

(This section requires Symbolic Math Toolbox.)

For a given square matrix $A$ there exists an invertible matrix $P$ (similarity matrix) such that ${P}^{-1}AP=J$ and $J$ is block diagonal matrix

$J=\left[\begin{array}{ccc}\hfill {J}_{1}\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \ddots \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill {J}_{k}\hfill \end{array}\right]$

where each block ${J}_{i}$ ($1\le i\le k$) is a $1×1$ matrix or it is a square matrix of the form of

${J}_{i}=\left[\begin{array}{cccc}\hfill {\lambda }_{i}\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill {\lambda }_{i}\hfill & \hfill \ddots \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill & \hfill 1\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill {\lambda }_{i}\hfill \end{array}\right]$

$J$ is such that the only non-zero entries of $J$ are on the diagonal and the super-diagonal. $J$ is called the Jordan normal form of $A$ and each ${J}_{i}$ is called a Jordan block of $A$. In a given Jordan block, every entry on the super-diagonal is 1. The Jordan decomposition is not unique.

We have the following properties:

• Counting multiplicity, the eigenvalues of $J$ (and therefore $A$) are the diagonal entries.
• Given an eigenvalue ${\lambda }_{i}$, its geometric multiplicity is the dimension of , and it is the number of Jordan blocks corresponding to ${\lambda }_{i}$.
• The sum of the sizes of all Jordan blocks corresponding to an eigenvalue ${\lambda }_{i}$ is its algebraic multiplicity.
• $A$ is diagonalizable if and only if, for every eigenvalue $\lambda$ of $A$, its geometric and algebraic multiplicities coincide.
• $A=PJ{P}^{-1}$ and the columns of $P$ are the generalized eigenvectors of $A$ (see Section 12.7).

The following $11×11$ Jordan block matrix

 $\left[\begin{array}{ccccccccccc}\hfill 2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill i\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill i\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill i\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill i\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 7\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 7\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 7\hfill \end{array}\right]$

has

one $1×1$ block corresponding to eigenvalue 2 (with algebraic multiplicity 1 and geometric multiplicity 1);
one $3×3$ block corresponding to eigenvalue 0 (with algebraic multiplicity 3 and geometric multiplicity 1);
two $2×2$ blocks corresponding to eigenvalue $i$ (with algebraic multiplicity 4 and geometric multiplicity 2);
one $3×3$ block corresponding to eigenvalue 7 (with algebraic multiplicity 3 and geometric multiplicity 1).

Matlab J=jordan(A) computes the Jordan canonical (normal) form of $A$, The matrix must be known exactly. Thus, its elements must be integers or ratios of small integers. Any errors in the input matrix may completely change the Jordan canonical form.

The command [P,J]=jordan(A) computes both $J$, the Jordan canonical form, and the similarity transform, $P$, whose columns are the generalized eigenvectors.

>> A=[5 4 2 1;0 1 -1 -1;-1 -1 3 0;1 1 -1 2]
A =
5     4     2     1
0     1    -1    -1
-1    -1     3     0
1     1    -1     2

>> [P,J]=jordan(A)
P =
1     1    -1     1
0     0     1    -1
-1     0     0     0
1     0     0     1

J =
4     1     0     0
0     4     0     0
0     0     1     0
0     0     0     2

12.8.6. Matrix Powers and Exponentials. If $A$ is a square matrix and $k$ a positive integer, then A^k effectively multiplies $A$ by itself $k-1$ times. If $A$ is a square and invertible matrix and $k$ a positive integer, then A^(-k) effectively multiplies ${A}^{-1}$ (i.e., inv(A), in Matlab notation) by itself $k-1$ times. Fractional powers, like A^(3/4) are also permitted.

The operator .^ produces element-by-element powers. For example, for A=[1 2 3;3 4 5;6 7 8],

>> A^2
ans =
25          31          37
45          57          69
75          96         117

>> A.^2
ans =
1            4          9
9           16         25
36           49         64

The command sqrtm(A) computes A^(1/2) by a more accurate algorithm. The m in sqrtm distinguished this function from sqrt(A) which, like A.^(1/2) does its job element-by-element.

The command expm(A) computes the computes the exponential of a square matrix $A$, i.e., ${e}^{A}$. This operation should be distinguished from exp(A) which computes the element-by-element exponential of $A$.

>> A=[1 2; 0 1]
A =
1              2
0              1

>> expm(A)
ans =
2.7183    5.4366
0    2.7183

Using Symbolic Math Toolbox additional symbolics operations are possible:

>> syms t;                          % defines symbolic variable t
>> f=inline(’expm(t*A)’)            % defines symbolic inline
% function in A and t
f =
Inline function:
f(A,t) = expm(t*A)

>> f(A,t)

ans =

[ exp(t), 2*t*exp(t)]
[      0,     exp(t)]

>> f(A,1)

ans =

2.7183    5.4366
0    2.7183

12.9. Characteristic Polynomial of a Matrix and Cayley-Hamilton Theorem.

Matlab command p = poly(A), where $A$ is $n×n$ matrix returns an $n+1$ element row vector whose elements are the coefficients of the characteristic polynomial of $A$:

$det\left(A-\lambda I\right)=0.$

The coefficients are ordered in descending powers: if a vector $c=\left({c}_{1},{c}_{2},\dots ,{c}_{n+1}\right)$ has $n+1$ components, the polynomial it represents is

${c}_{1}{\lambda }^{n}+\cdots +{c}_{n}\lambda +{c}_{n+1}.$

The command p = poly(r), where $r$ is a vector returns a row vector whose elements are the coefficients of the polynomial whose roots are the elements of $r$.

>> A=[1 2 3;4 5 6;7 8 0]
A =
1     2     3
4     5     6
7     8     0

>> format rat; p=poly(A)
p =
1    -6   -72   -27

>> A=[1 2 3;4 5 6;7 8 0]
A =
1     2     3
4     5     6
7     8     0

>> p=poly(sym(A)) % requires Symbolic Math Toolbox
% computes the characteristic polynomial
% of A in symbolic variable x
p =
x^3 - 6*x^2 - 72*x - 27

>> A=[1 2 3;4 5 6;7 8 0]
A =
1     2     3
4     5     6
7     8     0

>> format rat; p=poly(A)
p =
1    -6   -72   -27

>> pA = polyvalm(p,A)
pA =
*       *       *
*       *       *
*       *       *

The command polyvalm(p,A) evaluates a polynomial in a matrix sense. This is the same as substituting matrix A in the polynomial p. The zero matrix pA verifies Cayley-Hamilton theorem. This theorem states that if $p\left(\lambda \right)={p}_{0}+{p}_{1}\lambda +\dots +{\left(-1\right)}^{n}{\lambda }^{n}$ is the characterstic polynomial of $A$, then

$p\left(A\right)={p}_{0}I+{p}_{1}A+\dots +{\left(-1\right)}^{n}{A}^{n}=0.$

Finally, the last script below computes the roots of the characteristic polynomial p.

>> A=[1 2 3;4 5 6;7 8 0]
A =
1     2     3
4     5     6
7     8     0

>> format rat; p=poly(A)
p =
1    -6   -72   -27

>> format rat; roots(p)    % computes the roots of
ans =                      % the polynomial p, i.e.,
2170/179               % the eigenvalues of A
-648/113
-769/1980

### 13. Curve Fitting and Interpolation

Curve fitting is the process of constructing a curve, or mathematical function, that has the best fit (in a least square sense) to a series of data points, possibly subject to additional constraints. Interpolation is the method of constructing new data points within the range of a discrete set of known data points. Extrapolation is the process of constructing new data points outside a discrete set of known data points.

Matlab has many utilities for curve fitting and interpolation. The example below involves fitting the randomly corrupted sine function using linear, quadratic, and polynomial of degree 4 curves fitting.

>> x = (0: 0.25: 2.5);
>> y = sin(x) + rand(size(x))/5;
>> plot(x,y)

Go to your figure window and click on Tools, and select Basic Fitting. A new window opens with Basic Fitting menu. Check linear, quadratic, and 4th degree polynomial, together with the boxes for Show equations, Plot residuals, and Show norm of residuals. Select also Separate figure.

Note: Due to randomness of y = sin(x) + rand(size(x))/5, your results will differ slightly from the above interpolating polynomials, graphs and norms.

13.1. Hands-on Understanding of Curve Fitting Using Polynomial Functions.

p = polyfit(x,y,n) finds the coefficients of a polynomial p(x) of degree n that fits the data, p(x(i)) to y(i), in a least squares sense. The result p is a row vector of length n+1 containing the polynomial coefficients in descending powers of $x$, i.e, $p=\left[{p}_{1},{p}_{2}\dots ,{p}_{n},{p}_{n+1}\right]$ corresponding to

$p\left(x\right)={p}_{1}{x}^{n}+{p}_{2}{x}^{n-1}+\cdots +{p}_{n}x+{p}_{n+1}$

Next, the command polyval(p,x) evaluates the polynomial at the data points x(i) and generates the values y(i) such that

${y}_{i}={p}_{1}{x}_{i}^{n}+{p}_{2}{x}_{i}^{n-1}+\cdots +{p}_{n}{x}_{i}+{p}_{n+1}.$

The next example involves fitting the error function, erf(x), by a polynomial in x of degree $6$. In mathematics, the error function (also called the Gauss error function) is a special function (non-elementary) of sigmoid shape which occurs in probability, statistics, materials science, and partial differential equations. It is defined as:

$erf\left(x\right)=\frac{2}{\sqrt{\pi }}{\int }_{0}^{x}{e}^{-{t}^{2}}dt.$

>> x = (0: 0.1: 2.5)’;
>> y=erf(x);

First generate x, a (column) vector of equally spaced points in the interval; then evaluate erf(x) at those points.
>> p = polyfit(x,y,6)
p =
0.0084   -0.0983    0.4217   -0.7435    0.1471    1.1064    0.0004

Compute the coefficients in the approximating polynomial of degree 6.

There are seven coefficients and the polynomial is

$0.0084{x}^{6}-0.0983{x}^{5}+0.4217{x}^{4}-0.7435{x}^{3}+0.1471{x}^{2}+1.1064x+0.0004$

>> f = polyval(p,x);

Evaluates the polynomial at the data points.
>> table = [x y f y-f]
table =
0         0    0.0004   -0.0004
0.1000    0.1125    0.1119    0.0006
0.2000    0.2227    0.2223    0.0004
0.3000    0.3286    0.3287   -0.0001
0.4000    0.4284    0.4288   -0.0004
...

2.1000    0.9970    0.9969    0.0001
2.2000    0.9981    0.9982   -0.0001
2.3000    0.9989    0.9991   -0.0003
2.4000    0.9993    0.9995   -0.0002
2.5000    0.9996    0.9994    0.0002

The (partial) table is showing the data, fit, and error.

So, on the interval $\left[0,2.5\right]$, the fit is good to between three and four digits. Beyond this interval, the graph below shows that the polynomial behavior takes over and the approximation quickly deteriorates.

>> x = (0: 0.1: 5)’;
>> y = erf(x);
>> f = polyval(p,x);
>> plot(x,y,’o’,x,f,’-’)
>> axis([0  5  0  2])

13.2. Interpolation.

Given the data points $\left({x}_{i},{y}_{i}\right)$, find ${y}_{j}$ at desired ${x}_{j}\ne {x}_{i}$ from ${y}_{j}=f\left({x}_{j}\right)$, where $f$ is a continuous function to be found from interpolation. The command is ynew = interp1(x,y,xnew, method). The methods are:

’nearest’ Nearest neighbor interpolation

’linear’  Linear interpolation (default)

’spline’  Cubic spline interpolation

’pchip’   Piecewise cubic Hermite interpolation

’cubic’   (Same as ’pchip’)

The default method is linear.

As an example, we generate a coarse sine curve and interpolate over a finer abscissa.

x = 0:10;
y = sin(x);
xi = 0:.25:10;
yi = interp1(x,y,xi,’nearest’);
plot(x,y,’o’,xi,yi)

Uses nearest method.

The figures below show the results of different interpolating methods.

For other interpolating functions/methods check Matlab help on spline or interpft. For two- and three-dimensional analogs of interp1, check Matlab help on interp2 and interp3 commands, respectively.

### 14. Solving Nonlinear Algebraic Equations

The command fzero solves nonlinear equations involving one variable of the form:

$f\left(x\right)=0.$

In the first example, we calculate $\pi$ by finding the zero of the sine function near 3.

>> x0=fzero(’sin’,3)
x0 =
3.1416

>> x0=fzero(’my_f’,10)
x0 =
1.0193

my_f is from the function file defined in Section 5.
>> x = fzero(@cos,[1 2])
x =
1.5708

Finds the zero of cosine between 1 and 2. Use the command help function_handle to check what is the meaning of @ symbol in this command.
>> f = @(x)x.^3-2*x-5;
>> z = fzero(f,2)
z =
2.0946

An alternative way to find the zeros of the function $f\left(x\right)={x}^{3}-2x-5$ is to use an anonymous function @(x)x.^3-2*x-5 and find the zero near $2$.

Because $f\left(x\right)=1\cdot {x}^{3}+0\cdot {x}^{2}+\left(-2\right)\cdot {x}^{1}+\left(-5\right)\cdot {x}^{0}$ is a polynomial, the following command,

>> roots([1 0 -2 -5])
ans =
2.0946
-1.0473 + 1.1359i
-1.0473 - 1.1359i

finds the same zero and the remaining roots of this polynomial.

### 15. Numerical Integration - Quadrature

Numerical integration (often abbreviated to quadrature) is an algorithm for calculating the numerical value of a definite integral. Matlab command y = quad(’your_function’,a,b,tol) evaluates the integral of your_function of one variable from $a$ to $b$, using adoptive Simpson’s rule. The optional parameter tol specifies absolute tolerance. The default value is $1{0}^{-6}$,

If one uses the statement [y,evals] = quad(’your_function’,a,b,tol), the procedure also returns, evals, the number of function evaluations used to find the numerical value of the integral.

For example, in order to compute

$\underset{0}{\overset{2}{\int }}\frac{1}{{x}^{3}-2x-5}\phantom{\rule{0.3em}{0ex}}dx$

type the following:

y =
-0.4605

>> f = inline(’1./(x.^3-2*x-5)’);
y =
-0.4605
evals =
41

An alternative way to find the integral of the same function with the number of function evaluations provided as an outcome.
integral =
0.0917
evals =
125

The definite integral from $0$ to $1$ of my_f, the function file defined in Section 5.

$\underset{0.3}{\overset{0.7}{\int }}exp\left(-{x}^{2}\right)\phantom{\rule{0.3em}{0ex}}dx=\frac{\sqrt{\pi }}{2}\left[erf\left(0.7\right)-erf\left(0.3\right)\right],$

where $erf\left(x\right)$ is the internal Matlab function (see also Section 13.1),

>> format long
>> value = (sqrt(pi)/2)*(erf(0.7)-erf(0.3));
>> table = [value y yl (value-y)/value (value-yl)/value evals evalsl]’
table =
0.309447785425779
0.309447785513074
0.309447785425838
-0.000000000282099
-0.000000000000191
13.000000000000000
18.000000000000000

With the default accuracy quad gives only $\left(0.00028×1{0}^{-6}\right)%$ error, while quadl gives $\left(0.00000019×1{0}^{-6}\right)%$ error. The algorithm quad uses $13$ and quadl uses $18$ functions evaluations to compute the integral.

For double integrals,

$\underset{{y}_{min}}{\overset{{y}_{max}}{\int }}\left(\underset{{x}_{min}}{\overset{{x}_{max}}{\int }}f\left(x,y\right)\phantom{\rule{0.3em}{0ex}}dx\right)dy,$

Matlab provides q = dblquad(’your_function’,xmin,xmax,ymin,ymax,tol,@method), where @method can be @quad (default) or @quadl. You need [] as a placeholder for the default value of tol, when the optional argument @method is used. For example, for the double integral,

$\underset{0}{\overset{\pi }{\int }}\left[\underset{\pi }{\overset{2\pi }{\int }}\left(ysinx+xcosy\right)\phantom{\rule{0.3em}{0ex}}dx\right]dy=-{\pi }^{2},$

>> format long
>> value =-pi^2;
>> table = [value q ql (value-q)/value (value-ql)/value]’
table =
-9.869604401089358
-9.869604377254573
-9.869604289913852
0.000000002414969
0.000000011264434

Note that the lower order method (quad) performs better in this case than the higher order method (quadl).

### 16. Solving Ordinary Differential Equations

Matlab has many algorithms for solving systems of ordinary differential equations (ODE). Two procedures, ode23 and ode45, described below are implementations of 2nd/3rd-order and 4th/5th-order Runge-Kutta methods, respectively.

For the system of ODE in the vector form,

$\stackrel{̇}{x}=f\left(x,t\right),$

where

$\begin{array}{c}\hfill x={\left[{x}_{1}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{x}_{2}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\dots \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{x}_{n}\right]}^{T}=\left[\begin{array}{c}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \\ \hfill ⋮\hfill \\ \hfill {x}_{n}\hfill \end{array}\right],\end{array}\phantom{\rule{2em}{0ex}}\begin{array}{c}\hfill \stackrel{̇}{x}={\left[\frac{d{x}_{1}}{dt}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\frac{d{x}_{2}}{dt}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\dots \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\frac{d{x}_{n}}{dt}\right]}^{T}=\left[\begin{array}{c}\hfill \frac{d{x}_{1}}{dt}\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{d{x}_{2}}{dt}\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill ⋮\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{d{x}_{n}}{dt}\hfill \end{array}\right],\end{array}$

and

$\begin{array}{c}\hfill f={\left[{f}_{1}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{f}_{2}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\dots \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{f}_{n}\right]}^{T}=\left[\begin{array}{c}\hfill {f}_{1}\left({x}_{1},{x}_{2},\dots ,{x}_{n},t\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill {f}_{2}\left({x}_{1},{x}_{2},\dots ,{x}_{n},t\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill ⋮\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill {f}_{n}\left({x}_{1},{x}_{2},\dots ,{x}_{n},t\right)\hfill \end{array}\right],\end{array}$

one has to write a function that computes ${f}_{1},{f}_{2}\dots ,{f}_{n}$, given the input $\left(x,t\right)$, where $x$ is the column vector $x={\left[{x}_{1}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{x}_{2}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\dots \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{x}_{n}\right]}^{T}$. Your function must return the state derivative $\stackrel{̇}{x}$ as a column vector.

Here is the syntax for ode23 (for ode45, just replace ode23 with ode45):

(time, solution) = ode23(’your_function’, tspan, x0),

where tspan$=\left[{t}_{0},{t}_{final}\right]$ and x0 is the initial condition(s). For a system of $n$ equations, the output solution contains $n$ columns. Note which column corresponds to which variable in order to extract the correct column.

Example 1. For the first order differential equation,

 (5)

create the following function file:

function xdot = ode1(t,x);
% ODE1: computes xdot.
xdot = x.^2-t.^2./(t.^4+1);

and save it as ode1.m in the current working directory of Matlab.

Example 2a. The one-dimensional motion of an object projected vertically downward (toward the surface of the Earth) in the media offering resistance proportional to the square root of velocity can be modeled by the following differential equation:

 $m\frac{dv}{dt}=mg-k\sqrt{v},\phantom{\rule{1em}{0ex}}v\left(0\right)={v}_{0}\ge 0,$ (6)

where $v$ is velocity, $m$ is the mass of the object, $g=9.81\phantom{\rule{0.3em}{0ex}}\left(m∕se{c}^{2}\right)$ is the acceleration due to gravity, $k$ is the proportionality constant, and ${v}_{0}$ is the initial velocity. The positive direction of the velocity is in the downward direction.

In this example we are not interested in the positions of the moving object, though you can imagine that the object was thrown vertically downward from a very high building. Computation of the positions is done in Example 2b.

In the calculations below, we choose $m=1$, the drag constant $k=9$, and experiment with different (positive) initial velocities.

Create the following function file:

function dvdt = ode2(t,v);
% ODE2A: computes dvdt.
m =1; g = 9.81; k = 9;
dvdt = g-(k/m)*sqrt(v);

and save it as ode2a.m in the current working directory of Matlab.

The dashed curve in Figure 8 represents a function decreasing in time and thus, it does not describe the intended physical situation. Indeed, the terminal velocity for the problem described by (6) is given by

${V}_{terminal}={\left(\frac{mg}{k}\right)}^{2}={\left(\frac{9.81}{9}\right)}^{2}=1.1881,$

corresponding to one of the equilibrium solutions of (6): $v\left(t\right)={\left(mg∕k\right)}^{2}=1.1881$. Therefore, for initial velocities greater than $1.1881$, as is the case of the initial velocity ${v}_{0}=1.5$, the typical behavior of such solutions is indicated by the dashed curve in Figure 8.

Practice exercise. Write you own function file to solve the problem of an object projected vertically upward in the media offering resistance proportional to the square root of velocity (with $k=1$, $g=9.81$, and $m=1$). As before, choose the positive direction of the velocity in the downward direction and the initial velocity ${v}_{0}=-10$. Additionally, find the time ${t}_{0}$ such that $v\left({t}_{0}\right)=0$.

Note: The direction of the drag force always opposes the motion of the object, and thus it is different for upward/downward motions.

Example 2b.

In this example we consider the same physical problem as in Example 2a and, in addition to the velocity, we also compute the positions of the object of mass $m=1$ kg that is dropped from the height of 100 meters.

 $m\frac{{d}^{2}x}{d{t}^{2}}=mg-k\sqrt{\frac{dx}{dt}},\phantom{\rule{1em}{0ex}}x\left(0\right)=0,\phantom{\rule{1em}{0ex}}\frac{dx}{dt}\left(0\right)=0,$ (7)

We assume that $k=1$ and the positive direction of the velocity and the position is in the downward direction.

Equation (7) can be written as the system of two differential equations for ${z}_{1}=x$ and ${z}_{2}=dx∕dt=v$ as follows:

$\begin{array}{c}\hfill \frac{dz}{dt}=\left[\begin{array}{c}\hfill \frac{d{z}_{1}}{dt}\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{d{z}_{2}}{dt}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {z}_{2}\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill g-\frac{k}{m}\sqrt{{z}_{2}}\hfill \end{array}\right],\phantom{\rule{2em}{0ex}}\left[\begin{array}{c}\hfill {z}_{1}\left(0\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill {z}_{2}\left(0\right)\hfill \end{array}\right]=\left[\begin{array}{c}\hfill 0\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill 0\hfill \end{array}\right].\end{array}$

Create the following function file:

function dzdt = ode2b(t,z);
% ODE2B: computes dzdt.
g = 9.81; k = 1; m = 1;
dzdt = [z(2); g-(k/m)*sqrt(abs(z(2)))];

and save it as ode2b.m in the current working directory of Matlab.

Please note that in the commands x=z(:,1) and v=z(:,2), $x$ denotes the 1st column (i.e., the position) and $v$ denotes the 2nd column (i.e, the velocity) of the vector $z=\left[{z}_{1},{z}_{2}\right]$, respectively.

Next, we would like to compute the time needed for the object to reach the ground level. Let us look at the output:

>> [t,x]
ans =
0         0
0.0000    0.0000
0.0000    0.0000
0.0002    0.0000
0.0007    0.0000
0.0019    0.0000
0.0048    0.0001
0.0113    0.0006
0.0225    0.0024
0.0397    0.0075
0.0676    0.0214
0.1111    0.0572
0.1772    0.1432
0.2750    0.3391

0.4168    0.7626
0.6183    1.6368
0.8994    3.3665
1.2857    6.6575
1.8092   12.6973
2.5101   23.4179
3.4389   41.8660
4.4389   66.6306
5.4389   96.0028
6.4389  129.5746
7.4389  166.9978
8.4389  207.9697
9.4389  252.2230
10.0000  278.3993

We see that this event occurs between $5.4389$ and $6.4389$ seconds. In order to determine this time more precisely we use the interpolation from Section 13.2 to build a continuous function that represents the trajectory of the object (we use here $s$ variable to denote time since $t$ is already used as a variable). This can be done with the command interp1(t,x,s). Now, the function , where x_dist=interp1(t,x,s), is the height of the object as a function of time. Furthermore, with fzero command (nonlinear equation solver) from Section 14, we can compute the zero of $f\left(s\right)$, which will represent the time needed for the object to reach the ground level.

And here is the actual Matlab session to do the above calculations:

>> f = inline(’100-interp1(t,x,s)’)
f =
Inline function:
f(s,t,x) = 100-interp1(t,x,s)
>> fzero(@(s)f(s,t,x),5)
ans =
5.5580

This time can be compared to the time needed by the same object to reach the ground if there is no drag force:

$t=\sqrt{\frac{2h}{g}}=\sqrt{\frac{2\cdot 100}{9.81}}=4.5152.$

Practice exercise. Repeat Example 2b (the same mass and height, and with $k=1$) when the drag force is proportional to $v$ and to ${v}^{2}$.

Example 3.

The motion of the nonlinear pendulum is

 $\stackrel{̈}{\theta }+{\omega }^{2}sin\theta =0.$ (8)

We consider the initial conditions

$\theta \left(0\right)=0,\phantom{\rule{2em}{0ex}}\stackrel{̇}{\theta }\left(0\right)=1.$

Equation (8) can be written as the system of two differential equations for $z=\left[{z}_{1},{z}_{2}\right]$, with ${z}_{1}=\theta$ and ${z}_{2}=\stackrel{̇}{\theta }$, as follows:

$\begin{array}{c}\hfill \frac{dz}{dt}=\left[\begin{array}{c}\hfill \frac{d{z}_{1}}{dt}\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{d{z}_{2}}{dt}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {z}_{2}\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill -{\omega }^{2}sin\left({z}_{1}\right)\hfill \end{array}\right],\phantom{\rule{2em}{0ex}}\left[\begin{array}{c}\hfill {z}_{1}\left(0\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill {z}_{2}\left(0\right)\hfill \end{array}\right]=\left[\begin{array}{c}\hfill 0\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill 1\hfill \end{array}\right].\end{array}$

Create the following function file:

function dzdt = ode3(t,z);
% ODE3: computes dzdt.
w=1;
dzdt = [z(2); -w^2*sin(z(1))];

and save it as ode3.m in the current working directory of Matlab.

>> z0=[0,1];
>> tspan =[0,20];
>> [t,z] = ode45(’ode3’,tspan,z0);
>> x=z(:,1); v=z(:,2);
>> plot(t,x,t,v,’--’)
>> figure(2)
>> plot(x,v)

Please note that in the commands x=z(:,1); v=z(:,2); $x$ denotes the 1st column (i.e., the position $\theta$ ) and $v$ denotes the 2nd column (i.e, the velocity $\stackrel{̇}{\theta }$) of the vector $z=\left[{z}_{1},{z}_{2}\right]$, respectively.

16.1. Solving Linear Systems of First Order Differential Equations.

Consider a system of homogeneous linear first order differential equations

 $\stackrel{̇}{X}=AX,$ (9)

where

$\begin{array}{c}\hfill X={\left[{x}_{1}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{x}_{2}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\dots \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{x}_{n}\right]}^{T}=\left[\begin{array}{c}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \\ \hfill ⋮\hfill \\ \hfill {x}_{n}\hfill \end{array}\right],\end{array}\phantom{\rule{2em}{0ex}}\begin{array}{c}\hfill \stackrel{̇}{X}={\left[\frac{d{x}_{1}}{dt}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\frac{d{x}_{2}}{dt}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\dots \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\frac{d{x}_{n}}{dt}\right]}^{T}=\left[\begin{array}{c}\hfill \frac{d{x}_{1}}{dt}\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{d{x}_{2}}{dt}\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill ⋮\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{d{x}_{n}}{dt}\hfill \end{array}\right],\end{array}$

and $A$ is $n×n$ matrix with real entries.

16.1.1. Distinct real eigenvalues. In the case $A$ has $n$ real distinct eigenvalues ${\lambda }_{1},{\lambda }_{2},\dots ,{\lambda }_{n}$ and the corresponding (linearly independent) eigenvectors ${K}_{1},{K}_{2},\dots ,{K}_{n}$, the general solution of (9) can be written in the form

$X\left(t\right)={c}_{1}{X}_{1}\left(t\right)+{c}_{2}{X}_{2}\left(t\right)+\cdots +{c}_{n}{X}_{n}\left(t\right),$

where

${X}_{1}\left(t\right)=exp\left({\lambda }_{1}t\right){K}_{1},\phantom{\rule{1em}{0ex}}{X}_{2}\left(t\right)=exp\left({\lambda }_{2}t\right){K}_{2},\phantom{\rule{1em}{0ex}}\dots ,\phantom{\rule{1em}{0ex}}{X}_{n}\left(t\right)=exp\left({\lambda }_{n}t\right){K}_{n}\phantom{\rule{0.3em}{0ex}},$

and constants ${c}_{1},{c}_{2},\dots ,{c}_{n}$ are determined from the initial condition at $t={t}_{0}$ on $X$, i.e., $X\left({t}_{0}\right)={X}_{0}$, where vector ${X}_{0}$ is provided.

>> A=[-4 1 1;1 5 -1;0 1 -3]

A =
-4        1        1
1        5       -1
0        1       -3

>> format rat; [V,D]=eig(A)
V =

-129/1048       201/203        985/1393
-129/131       -201/2030         *
-129/1048       201/2030       985/1393

D =

5              0              0
0             -4              0
0              0             -3

The columns of matrix $V$ provide the eigenvectors of $A$, where * in the last column of $V$ signifies $0$. Eigenvectors are determined up to multiplicative scalars and since $129∕131=8\cdot \left(129∕1048\right)$ and $201∕203=10\cdot \left(201∕2030\right)$, the general solution is

$X\left(t\right)={c}_{1}exp\left(5t\right)\left[\begin{array}{c}\hfill 1\hfill \\ \hfill 8\hfill \\ \hfill 1\hfill \end{array}\right]+{c}_{2}exp\left(-4t\right)\left[\begin{array}{c}\hfill 10\hfill \\ \hfill -1\hfill \\ \hfill 1\hfill \end{array}\right]+{c}_{3}exp\left(-3t\right)\left[\begin{array}{c}\hfill 1\hfill \\ \hfill 0\hfill \\ \hfill 1\hfill \end{array}\right].$

For the initial condition ${X}_{0}={\left(1,2,3\right)}^{T}$, the constants ${c}_{1}$, ${c}_{2}$, and ${c}_{3}$ are given by

>> B=[1 10 1;8 -1 0;1 1 1]
B =
1     10      1
8     -1      0
1      1      1

>> X0=[1;2;3]
X0 =
1
2
3

>> format rat; C=B\X0
C =
2/9     % c1
-2/9     % c2
3       % c3

and the solution is

The general solution to (9) has also the following form,

$X\left(t\right)=exp\left(tA\right)C,$

where $C={\left({c}_{1},{c}_{2},\dots ,{c}_{n}\right)}^{T}$ and $exp\left(tA\right)$ is the exponential of the matrix $tA$. In particular, for

$\stackrel{̇}{X}=AX,\phantom{\rule{2em}{0ex}}X\left(0\right)={X}_{0}\phantom{\rule{0.3em}{0ex}},$

the solution is given by

$X\left(t\right)=exp\left(tA\right){X}_{0}\phantom{\rule{0.3em}{0ex}}.$

Using the above example, the following Matlab session (that uses Symbolic Math Toolbox)

>> A=[-4 1 1;1 5 -1;0 1 -3]
A =
-4        1        1
1        5       -1
0        1       -3

>> X0=[1;2;3]
X0 =
1
2
3

>> syms t    % defines symbolic variable t

>> f=inline(’expm(t*A)’)

f =
Inline function:
f(A,t) = expm(t*A)

>> X=f(A,t)*X0
X =
3/exp(3*t) - 20/(9*exp(4*t)) + (2*exp(5*t))/9
2/(9*exp(4*t)) + (16*exp(5*t))/9
3/exp(3*t) - 2/(9*exp(4*t)) + (2*exp(5*t))/9

provides the solution already obtained in (10).

16.1.2. Repeated real eigenvalues. If $n×n$ matrix $A$ in the system (9)has only $k linearly independent eigenvectors ${K}_{1}$, ${K}_{2}$,…${K}_{k}$ then, we have only $k$ linearly independent solutions of the form $exp\left({\lambda }_{i}t\right){K}_{i}$. $i=1,\dots k$. To find additional solutions we pick an eigenvalue $\lambda$ of $A$ and find all vectors $v$ for which ${\left(A-\lambda I\right)}^{2}v=0$, but $\left(A-\lambda I\right)v\ne 0$. For such vectors $v$,

$exp\left(At\right)v=exp\left(\lambda t\right)exp\left(A-\lambda I\right)tv=exp\left(\lambda t\right)\left[v+t\left(A-\lambda I\right)v\right]$

is an additional solution of (9).

If we still don’t have enough solutions, then we find all vectors $v$ for which ${\left(A-\lambda I\right)}^{3}v=0$, but ${\left(A-\lambda I\right)}^{2}v\ne 0$. For such vectors $v$,

$exp\left(At\right)v=exp\left(\lambda t\right)exp\left(A-\lambda I\right)tv=exp\left(\lambda t\right)\left[v+t\left(A-\lambda I\right)v+\frac{{t}^{2}}{2!}{\left(A-\lambda I\right)}^{2}v\right]$

is an additional solution of (9). We keep proceeding in this manner until we obtain $n$ linearly independent solutions. The result from linear algebra guarantees that this algorithm works, Indeed, if $n×n$ matrix $A$ has $k$ distinct eigenvalues ${\lambda }_{1},{\lambda }_{2}\dots ,{\lambda }_{k}$ with multiplicity ${n}_{1}\dots ,{n}_{k}$, then for each $j=1,\dots k$, there exists ${d}_{j}\le {n}_{j}$ such that the equation ${\left(A-\lambda I\right)}^{{d}_{j}}v=0$ has at least ${n}_{j}$ linearly independent solutions. Thus, for each eigenvalue ${\lambda }_{j}$ of $A$, we can compute ${n}_{j}$ linearly solutions of (9). All these solutions have the form

$X\left(t\right)=exp\left({\lambda }_{j}t\right)\left[v+t\left(A-\lambda I\right)v+\dots +\frac{{t}^{{d}_{j}-1}}{\left({d}_{j}-1\right)!}{\left(A-\lambda I\right)}^{{d}_{j}-1}v\right].$

Furthermore, ${n}_{1}+{n}_{2}+\dots {n}_{k}=n$ and such obtained solutions must be linearly independent. We know from Section 12.7 that vectors $v$ satisfying ${\left(A-\lambda I\right)}^{{d}_{j}}v=0$ are the generalized eigenvectors of $A$.

Consider the system (9) with $A$ given by

$A=\left[\begin{array}{cccc}\hfill 1\hfill & \hfill -2\hfill & \hfill 3\hfill & \hfill -1\hfill \\ \hfill 1\hfill & \hfill -5\hfill & \hfill 7\hfill & \hfill -2\hfill \\ \hfill 2\hfill & \hfill -9\hfill & \hfill 10\hfill & \hfill -2\hfill \\ \hfill 2\hfill & \hfill -8\hfill & \hfill 6\hfill & \hfill 1\hfill \end{array}\right]$

The command [P,J]=jordan(A) computes both $J$, the Jordan canonical form, and the similarity transform, $P$, whose columns are the generalized eigenvectors.

>> A=[1 -2 3 -1; 1 -5 7 -2; 2 -9 10 -2;2 -8 6 1]
A =
1   -2   3   -1
1   -5   7   -2
2   -9  10   -2
2   -8   6    1

>> format rat; [P,J]=jordan(A)
P =
4   -1   3   -3
4   -2   5   -4
4   -3   6   -4
4   -4   6   -4

J =
1    0   0    0
0    2   1    0
0    0   2    1
0    0   0    2

Let us denote the columns of $P$ by ${v}_{1}$, ${v}_{2}$, ${v}_{3}$, ${v}_{4}$. We have,

$\left(A-I\right){v}_{1}=0,\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\left(A-2I\right){v}_{2}=0,\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\left(A-2I\right){v}_{3}={v}_{2},\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\left(A-2I\right){v}_{4}={v}_{3},\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{\left(A-2I\right)}^{2}{v}_{3}=0,\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{\left(A-2I\right)}^{3}{v}_{4}=0,$

as also verified by the following Matlab commands:

>> (A-eye(4))*P(:,1)

ans =
0
0
0
0

>> (A-2*eye(4))*P(:,2)

ans =
0
0
0
0

>> (A-2*eye(4))*P(:,3)

ans =
-1
-2
-3
-4

>> (A-2*eye(4))*P(:,4)

ans =
3
5
6
6

>> (A-2*eye(4))^2*P(:,3)

ans =
0
0
0
0

>> (A-2*eye(4))^3*P(:,4)

ans =
0
0
0
0

The solution of (9) corresponding to $\lambda =1$ (eigenvalue of multiplicity 1) is

${X}_{1}\left(t\right)=exp\left(t\right){v}_{1}=exp\left(t\right)\left[\begin{array}{c}\hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \end{array}\right].$

Three linearly independent solutions of (9) corresponding to $\lambda =2$ (eigenvalue of multiplicity 3) are

${X}_{2}\left(t\right)=exp\left(2t\right){v}_{1}=exp\left(2t\right)\left[\begin{array}{c}\hfill -1\hfill \\ \hfill -2\hfill \\ \hfill -3\hfill \\ \hfill -4\hfill \end{array}\right],\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{X}_{3}\left(t\right)=exp\left(2t\right)\left[I+t\left(A-2I\right)\right]{v}_{3}=exp\left(2t\right)\left[{v}_{3}+t{v}_{2}\right]=exp\left(2t\right)\left[\begin{array}{c}\hfill 3-t\hfill \\ \hfill 5-2t\hfill \\ \hfill 6-3t\hfill \\ \hfill 6-4t\hfill \end{array}\right],\phantom{\rule{1em}{0ex}}$

and

${X}_{4}\left(t\right)=exp\left(2t\right)\left[I+t\left(A-2I\right)+\frac{{t}^{2}}{2!}{\left(A-2I\right)}^{2}\right]{v}_{4}=exp\left(2t\right)\left[{v}_{4}+t{v}_{3}+\frac{{t}^{2}}{2!}{v}_{2}\right]=exp\left(2t\right)\left[\begin{array}{c}\hfill 3t-{t}^{2}∕2-3\hfill \\ \hfill 5t-{t}^{2}-4\hfill \\ \hfill 6t-\left(3∕2\right){t}^{2}-4\hfill \\ \hfill 6t-2{t}^{2}-4\hfill \end{array}\right]$

and the general solution is

$X\left(t\right)={c}_{1}{X}_{1}\left(t\right)+{c}_{2}{X}_{2}\left(t\right)+{c}_{3}{X}_{3}\left(t\right)+{c}_{4}{X}_{4}\left(t\right)$

For the initial condition ${X}_{0}={\left(1,0,0,0\right)}^{T}$, the constants ${c}_{1}$, ${c}_{2}$, ${c}_{3}$, and ${c}_{4}$ are solutions of the system:

$\left[\begin{array}{c}\hfill 1\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right]={X}_{0}={c}_{1}{v}_{1}+{c}_{2}{v}_{2}+{c}_{3}{v}_{3}+{c}_{4}{v}_{4}=\left[\begin{array}{cccc}\hfill 4\hfill & \hfill -1\hfill & \hfill 3\hfill & \hfill -3\hfill \\ \hfill 4\hfill & \hfill -2\hfill & \hfill 5\hfill & \hfill -4\hfill \\ \hfill 4\hfill & \hfill -3\hfill & \hfill 6\hfill & \hfill -4\hfill \\ \hfill 4\hfill & \hfill -4\hfill & \hfill 6\hfill & \hfill -4\hfill \end{array}\right]\left[\begin{array}{c}\hfill {c}_{1}\hfill \\ \hfill {c}_{2}\hfill \\ \hfill {c}_{3}\hfill \\ \hfill {c}_{4}\hfill \end{array}\right]$

Using Matlab, we obtain

>> X0=[1;0;0;0]

X0 =
1
0
0
0

>> format rat; C=P\X0

C =
1      % c1
0      % c2
0      % c3
1      % c4

and the solution of (9) with the initial condition ${X}_{0}={\left(1,0,0,0\right)}^{T}$ is

 $X\left(t\right)={X}_{1}\left(t\right)+{X}_{4}\left(t\right)=\left[\begin{array}{c}\hfill 4exp\left(t\right)+\left(3t-{t}^{2}∕2-3\right)exp\left(2t\right)\hfill \\ \hfill 4exp\left(t\right)+\left(5t-{t}^{2}-4\right)exp\left(2t\right)\hfill \\ \hfill 4exp\left(t\right)+\left(6t-\left(3∕2\right){t}^{2}-4\right)exp\left(2t\right)\hfill \\ \hfill 4exp\left(t\right)+\left(6t-2{t}^{2}-4\right)exp\left(2t\right)\hfill \end{array}\right].$ (11)

Finally, the following Matlab session (that uses Symbolic Math Toolbox)

>> A=[1 -2 3 -1; 1 -5 7 -2; 2 -9 10 -2;2 -8 6 1]
A =
1   -2   3   -1
1   -5   7   -2
2   -9  10   -2
2   -8   6    1

>> X0=[1;0;0;0]
X0 =
1
0
0
0

>> syms t   % defines symbolic variable t

>> f=inline(’expm(t*A)’)

f =
Inline function:
f(A,t) = expm(t*A)

>> X=f(A,t)*X0

X =
4*exp(t) - 3*exp(2*t) + 3*t*exp(2*t) - (t^2*exp(2*t))/2
4*exp(t) - 4*exp(2*t) + 5*t*exp(2*t) - t^2*exp(2*t)
4*exp(t) - 4*exp(2*t) + 6*t*exp(2*t) - (3*t^2*exp(2*t))/2
4*exp(t) - 4*exp(2*t) + 6*t*exp(2*t) - 2*t^2*exp(2*t)

verifies the solution already obtained in (11).

16.1.3. Complex eigenvalues. If $\lambda =\alpha +i\beta$ is a complex eigenvalue of multiplicity 1 of the coefficient matrix $A$ (with real entries) and $K$ is the corresponding eigenvector, then

are complex-valued solutions of (9), where $\stackrel{̄}{\lambda }=\alpha -i\beta$ and $\overline{K}$ is the conjugate vector of $K$. Furthermore, if

then

$\begin{array}{llll}\hfill {X}_{1}\left(t\right)& =\left[{B}_{1}cos\left(\beta t\right)-{B}_{2}sin\left(\beta t\right)\right]exp\left(\alpha t\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill {X}_{2}\left(t\right)& =\left[{B}_{2}cos\left(\beta t\right)+{B}_{1}sin\left(\beta t\right)\right]exp\left(\alpha t\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$

are linearly independent real-valued solutions of (9) corresponding to the eigenvalue $\lambda$ of multiplicity 1.

Consider (9) with $A=\left[\begin{array}{cc}\hfill 4\hfill & \hfill -5\hfill \\ \hfill 5\hfill & \hfill -4\hfill \end{array}\right]$. After using [P,J]=jordan(A) command,

>> A=[4 -5;5 -4]

A =
4     -5
5     -4

>> format rat; [P,J]=jordan(A)

P =
1/2 + 2/3i     1/2 - 2/3i
0   + 5/6i      0  - 5/6i

J =
0 - 3i           0
0            0 + 3i

>> B1=(1/2)*(P(:,1)+conj(P(:,1)))

B1 =
1/2
0

>> B2=(i/2)*(conj(P(:,1))-P(:,1))

B2 =
2/3
5/6

we obtain

$\begin{array}{llll}\hfill {X}_{1}\left(t\right)& ={B}_{1}cos\left(3t\right)+{B}_{2}sin\left(3t\right)=\left[\begin{array}{c}\hfill \frac{1}{2}cos\left(3t\right)+\frac{2}{3}sin\left(3t\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{5}{6}sin\left(3t\right)\hfill \end{array}\right]\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill {X}_{2}\left(t\right)& ={B}_{2}cos\left(3t\right)-{B}_{1}sin\left(3t\right)=\left[\begin{array}{c}\hfill \frac{2}{3}cos\left(3t\right)-\frac{1}{2}sin\left(3t\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{5}{6}cos\left(3t\right)\hfill \end{array}\right]\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$ The general solution is a linear combination of ${X}_{1}$ and ${X}_{2}$ in the form
$X\left(t\right)={c}_{1}\left[\begin{array}{c}\hfill \frac{1}{2}cos\left(3t\right)+\frac{2}{3}sin\left(3t\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{5}{6}sin\left(3t\right)\hfill \end{array}\right]+{c}_{2}\left[\begin{array}{c}\hfill \frac{2}{3}cos\left(3t\right)-\frac{1}{2}sin\left(3t\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{5}{6}cos\left(3t\right)\hfill \end{array}\right]$

With ${X}_{0}={\left(1,0\right)}^{T}$, ${c}_{1}=2$ and ${c}_{2}=0$ and the solution is

$X\left(t\right)=\left[\begin{array}{c}\hfill cos\left(3t\right)+\frac{4}{3}sin\left(3t\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \frac{5}{3}sin\left(3t\right)\hfill \end{array}\right].$

One can easily plot phase portrait of the solutions using Matlab expm(A) and for loop commands:

>> X=[]; Y=[]; Z=[];
>> for t=0:.01:10
X=[X expm(t*A)*[1;0]]; Y=[Y expm(t*A)*[1;3]]; Z=[Z expm(t*A)*[10;-2]];
end
>> plot(X(1,:),X(2,:),’r’,Y(1,:),Y(2,:),’b’,Z(1,:),Z(2,:),’g’)

And here is another example of use of expm(A) and three dimensional phase plot (plot3) commands:

16.1.4. Nonhomogeneous systems. The solution of the nohomogeneous system

 $\stackrel{̇}{X}=AX+f,\phantom{\rule{2em}{0ex}}X\left({t}_{0}\right)={X}_{0},$ (12)
is given by the variation of parameters formula,
 $X\left(t\right)=exp\left[\left(t-{t}_{0}\right)A\right]{X}_{0}+\underset{{t}_{0}}{\overset{t}{\int }}exp\left[\left(t-s\right)A\right]f\left(s\right)\phantom{\rule{0.3em}{0ex}}ds.\phantom{\rule{2em}{0ex}}t\ge {t}_{0}.$ (13)

Here, $f\left(t\right)$ is a given vector valued function of $t$.

Although the method of variation of parameters is often a time consuming process if done by hand, with Matlab, it can be easily performed.

Consider the nonhomogeneous system,

 $\stackrel{̇}{X}=\left[\begin{array}{cccc}\hfill 2\hfill & \hfill -2\hfill & \hfill 2\hfill & \hfill 1\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill -1\hfill & \hfill 3\hfill & \hfill 0\hfill & \hfill 3\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 4\hfill & \hfill -2\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 2\hfill & \hfill -1\hfill \end{array}\right]X+\left[\begin{array}{c}\hfill texp\left(t\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill exp\left(-t\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill exp\left(2t\right)\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill 1\hfill \end{array}\right],\phantom{\rule{2em}{0ex}}X\left(0\right)=\left[\begin{array}{c}\hfill 1\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill 0\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill 1\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill 0\hfill \end{array}\right].$ (14)

First, we find the general solution of (14).

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

A =
2     -2      2      1
-1      3      0      3
0      0      4     -2
0      0      2     -1

>> format rat; [P,J]=jordan(A)

P =
1     -8/3    4      2/3
2/3    4/3   4/3    -2/3
-1/6     0    8/3      0
-1/3     0    4/3      0

J =
0      0      0       0
0      1      0       0
0      0      3       0
0      0      0       4

The eigenvalues are 0, 1, 3, and 4, with the corresponding eigenvectors

A particular solution is given by

${X}_{p}\left(t\right)=exp\left(tA\right)\int exp\left(-tA\right)f\left(t\right)\phantom{\rule{0.3em}{0ex}}dt.$

The corresponding Matlab script is,

>> syms t
>> Xp = expm(t*A)*int(expm(-t*A),t);    % The output is a very long
% expression and thus initially
% is suppressed.

>> Xp = simplify(Xp)                    % simplify command makes Xp
% much simpler

Xp =
(t^2*exp(t))/3 - 1/(5*exp(t)) - 5*exp(2*t) - exp(t)/27 - 4*t - (t*exp(t))/9 - 59/12
exp(t)/27 - 3/(10*exp(t)) - 2*exp(2*t) - (8*t)/3 + (t^2*exp(t))/6 + (t*exp(t))/9 - 95/36
(2*t)/3 - (3*exp(2*t))/2 + 2/9
(4*t)/3 - exp(2*t) + 1/9

Thus, the general solution of (14) is

$X\left(t\right)={c}_{1}\left[\begin{array}{c}\hfill 1\hfill \\ \hfill 2∕3\hfill \\ \hfill -1∕6\hfill \\ \hfill -1∕3\hfill \end{array}\right]+{c}_{2}exp\left(t\right)\left[\begin{array}{c}\hfill -8∕3\hfill \\ \hfill 4∕3\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right]+{c}_{3}exp\left(3t\right)\left[\begin{array}{c}\hfill 4\hfill \\ \hfill 4∕3\hfill \\ \hfill 8∕3\hfill \\ \hfill 4∕3\hfill \end{array}\right]+{c}_{4}exp\left(4t\right)\left[\begin{array}{c}\hfill 2∕3\hfill \\ \hfill -2∕3\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right]+{X}_{p}\left(t\right),$

where

${X}_{p}\left(t\right)=\left[\begin{array}{c}\hfill \left({t}^{2}∕3\right)exp\left(t\right)-\left(1∕5\right)exp\left(-t\right)-5exp\left(2t\right)-\left(1∕27\right)exp\left(t\right)-4t-\left(t∕9\right)exp\left(t\right)-59∕12\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \left(1∕27\right)exp\left(t\right)-\left(3∕10\right)exp\left(-t\right)-2exp\left(2t\right)-\left(8∕3\right)t+\left({t}^{2}∕6\right)exp\left(t\right)+\left(t∕9\right)exp\left(t\right)-95∕36\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \left(2∕3\right)t-\left(3∕2\right)exp\left(2t\right)+2∕9\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \left(4∕3\right)t-exp\left(2t\right)+1∕9\hfill \end{array}\right]$

Finally, the following Matlab script computes the solution to the initial value problem (14):

>> syms s t
>> X=expm(t*A)*[1;0;1;0]+expm(t*A)*int(expm(-s*A)*[s*exp(s);exp(-s);exp(2*s);1],0,t);

>> X=simplify(X)

X =

(8*exp(3*t))/3 - 1/(5*exp(t)) - 5*exp(2*t) - 4*t - (457*exp(4*t))/540
+ (332*exp(t))/27 + (t^2*exp(t))/3 - (t*exp(t))/9 - 95/12

(8*exp(3*t))/9 - 3/(10*exp(t)) - 2*exp(2*t) - (8*t)/3 + (457*exp(4*t))/540
+ (335*exp(t))/54 + (t^2*exp(t))/6 + (t*exp(t))/9 -167/36

(2*t)/3 - (3*exp(2*t))/2 + (16*exp(3*t))/9 + 13/18

(4*t)/3 - exp(2*t) + (8*exp(3*t))/9 + 10/9

The solution is

$X\left(t\right)=\left[\begin{array}{c}\hfill \left(8∕3\right){e}^{3t}-\left(1∕5\right){e}^{-t}-5{e}^{2t}-4t-\left(457∕540\right){e}^{4t}+\left(332∕27\right){e}^{t}+\left({t}^{2}∕3\right){e}^{t}-\left(t∕9\right){e}^{t}-95∕12\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \left(8∕3\right){e}^{3t}-\left(3∕10\right){e}^{-t}-2{e}^{2t}-\left(8∕3\right)t+\left(457∕540\right){e}^{4t}+\left(335∕54\right){e}^{t}+\left({t}^{2}∕6\right){e}^{t}+\left(t∕9\right){e}^{t}-167∕36\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \left(2∕3\right)t-\left(3∕2\right){e}^{2t}+\left(16∕9\right){e}^{3t}+13∕18\hfill \\ & & & & & & & & & & & & & & & & & & & \\ \hfill \left(4∕3\right)t-{e}^{2t}+\left(8∕9\right){e}^{3t}+10∕9\hfill \end{array}\right].$

16.2. Linearization of a Double Pendulum.

Consider the double-pendulum system consisting of a pendulum attached to a pendulum shown in Figure 13. The system oscillates in a vertical plane under influence of gravity, the mass of each rod is negligible, and there are no dumping forces acting on the system. The displacement ${\theta }_{1}$ is measured (in radians) from vertical line extending downward from the pivot of the system. The displacement ${\theta }_{2}$ is measured from a vertical line extending downward from the center of mass ${m}_{1}$. The positive direction is to the right, the negative direction is to the left. The system of nonlinear differential equations describing the motion is:

 $\begin{array}{c}\hfill \left({m}_{1}+{m}_{2}\right){L}_{1}^{2}{\stackrel{̈}{\theta }}_{1}+{m}_{2}{L}_{1}{L}_{2}{\stackrel{̈}{\theta }}_{2}cos\left({\theta }_{1}-{\theta }_{2}\right)++{m}_{2}{L}_{1}{L}_{2}{\left({\stackrel{̇}{\theta }}_{2}\right)}^{2}sin\left({\theta }_{1}-{\theta }_{2}\right)+\left({m}_{1}+{m}_{2}\right){L}_{1}gsin{\theta }_{1}=0;\\ \hfill {m}_{2}{L}_{2}^{2}{\stackrel{̈}{\theta }}_{2}+{m}_{2}{L}_{1}{L}_{2}{\stackrel{̈}{\theta }}_{1}cos\left({\theta }_{1}-{\theta }_{2}\right)-{m}_{2}{L}_{1}{L}_{2}{\left({\stackrel{̇}{\theta }}_{1}\right)}^{2}sin\left({\theta }_{1}-{\theta }_{2}\right)+{m}_{2}{L}_{2}gsin{\theta }_{2}=0.\end{array}$ (15)

If the displacements ${\theta }_{1}\left(t\right)$ and ${\theta }_{2}\left(t\right)$ are assumed to be small, then $cos\left({\theta }_{1}-{\theta }_{2}\right)\approx 1$, $sin\left({\theta }_{1}-{\theta }_{2}\right)\approx 0$, $sin\left({\theta }_{1}\right)\approx {\theta }_{1}$, $sin\left({\theta }_{2}\right)\approx {\theta }_{2}$, and the following linearization of (15) is obtained:

 $\begin{array}{c}\hfill \left({m}_{1}+{m}_{2}\right){L}_{1}^{2}{\stackrel{̈}{\theta }}_{1}+{m}_{2}{L}_{1}{L}_{2}{\stackrel{̈}{\theta }}_{2}+\left({m}_{1}+{m}_{2}\right){L}_{1}g{\theta }_{1}=0;\\ \hfill {m}_{2}{L}_{2}^{2}{\stackrel{̈}{\theta }}_{2}+{m}_{2}{L}_{1}{L}_{2}{\stackrel{̈}{\theta }}_{1}+{m}_{2}{L}_{2}gsin{\theta }_{2}=0.\end{array}$ (16)

Consider (16) with ${m}_{1}=3$, ${m}_{2}=1$, ${L}_{1}={L}_{2}=16$, $g=32$, and the initial conditions: ${\theta }_{1}\left(0\right)=1$, ${\theta }_{2}\left(0\right)=-1$, $\stackrel{̇}{{\theta }_{1}}\left(0\right)=0$, and $\stackrel{̇}{{\theta }_{2}}\left(0\right)=0$. We obtain the following system:

 $\begin{array}{c}\hfill 4{\stackrel{̈}{\theta }}_{1}+{\stackrel{̈}{\theta }}_{2}+8{\theta }_{1}=0;\\ \hfill {\stackrel{̈}{\theta }}_{1}+{\stackrel{̈}{\theta }}_{2}+2{\theta }_{2}=0.\end{array}$ (17)

Finally, with the $z={\left({z}_{1},{z}_{2},{z}_{3},{z}_{4}\right)}^{T}$ and ${z}_{1}={\theta }_{1}$, ${z}_{2}={\stackrel{̇}{\theta }}_{1}$, ${z}_{3}={\theta }_{2}$, and ${z}_{4}={\stackrel{̇}{\theta }}_{2}$, the system of second order differential equations becomes the system of first order differential equations:

 (18)

We compute the eigenvalues of $A$ by computing the roots of the characteristic polynomial $det\left(A-\lambda I\right)=0$.

>> A = [0 1 0 0;-8/3 0 2/3 0;0 0 0 1;8/3 0 -8/3 0]

A =
0       1       0         0
-8/3     0      2/3        0
0       0       0         1
8/3     0     -8/3        0

>> factor(poly(sym(A)))

ans =

((x^2 + 4)*(3*x^2 + 4))/3

The eigenvalues of matrix $A$ are $±2i$ and $±\left(2∕\sqrt{3}\right)i$. We can find the solution $z\left(t\right)$ of (18) by computing $z\left(t\right)=exp\left(tA\right){z}_{0}$. Here is the corresponding Matlab script:

>> z0=[1;0;-1;0]

z0 =
1
0
-1
0

>> syms t; f=inline(’expm(t*A)’)

f =

Inline function:
f(A,t) = expm(t*A)

>> z=f(A,t)*z0

z =
1/(8*exp((2*3^(1/2)*i*t)/3)) + exp((2*3^(1/2)*i*t)/3)/8 + 3/(8*exp(2*i*t)) + (3*exp(2*i*t))/8

(3*i*exp(2*i*t))/4 - (3*i)/(4*exp(2*i*t)) - (3^(1/2)*i)/(12*exp((2*3^(1/2)*i*t)/3)) +
(3^(1/2)*i*exp((2*3^(1/2)*i*t)/3))/12

1/(4*exp((2*3^(1/2)*i*t)/3)) + exp((2*3^(1/2)*i*t)/3)/4 - 3/(4*exp(2*i*t)) - (3*exp(2*i*t))/4

(3*i)/(2*exp(2*i*t)) - (3*i*exp(2*i*t))/2 - (3^(1/2)*i)/(6*exp((2*3^(1/2)*it)/3)) +
(3^(1/2)*i*exp((2*3^(1/2)*i*t)/3))/6

The first and third rows of z provide expressions for ${\theta }_{1}\left(t\right)$ and ${\theta }_{2}\left(t\right)$, respectively:

 $\begin{array}{c}\hfill {\theta }_{1}\left(t\right)=\frac{1}{4}cos\left(\frac{2}{\sqrt{3}}t\right)+\frac{3}{4}cos\left(2t\right)\\ \hfill {\theta }_{2}\left(t\right)=\frac{1}{2}cos\left(\frac{2}{\sqrt{3}}t\right)-\frac{3}{2}cos\left(2t\right)\end{array}$ (19)

From (19) , we observe that the motion of the (linearized) pendulums is not periodic since $cos\left(2t∕\sqrt{3}\right)$ has period $\sqrt{3}\pi$, $cos\left(2t\right)$ has period $\pi$, and the ratio of these periods is $\sqrt{3}$, which is not a rational number.

The following Matlab script computes the plots of ${\theta }_{1}\left(t\right)$, ${\theta }_{2}\left(t\right)$, and the plot of ${\theta }_{2}$ versus ${\theta }_{1}$ (i.e., the Lissajous curve of compound pendulum).

>> A = [0 1 0 0;-8/3 0 2/3 0;0 0 0 1;8/3 0 -8/3 0]

A =
0       1       0         0
-8/3     0      2/3        0
0       0       0         1
8/3     0     -8/3        0

>> X=[];
for t=0:.1:30
X=[X expm(t*A)*[1;0;-1;0]];
end

>> t=0:.1:30;
plot(t,X(1,:))           % plot of theta1
figure(2)
>> plot(t,X(3,:))        % plot of theta2
figure(3)
>> plot(X(1,:),X(3,:))   % plot of X(3,:) versus X(1,:)

And finally, the link below provides animation of a double compound nonlinear pendulum showing its chaotic behaviour:

### 17. Use of Global Variables

Occasionally, it is desirable to declare variables globally. Such variables will be accessible to all or some functions without passing these variables in the input list. The command global u v w declares variables u, v, and w globally. This statement must go before any executable in the functions and scripts that will use them. Global declarations of variables should be used with a caution.

Here is one example of a global declaration use. In Example 2a of Section 16, the one-dimensional motion of an object projected vertically downward (toward the surface of the Earth) in the media offering resistance proportional to the square root of velocity was modeled by the equation:

$m\frac{dv}{dt}=mg-k\sqrt{v},\phantom{\rule{1em}{0ex}}v\left(0\right)={v}_{0}\ge 0,$

where $v$ is velocity, $m$ is the mass of the object, $g=9.81\phantom{\rule{0.3em}{0ex}}\left(m∕se{c}^{2}\right)$ is the acceleration due to gravity, $k$ is the proportionality constant, and ${v}_{0}$ is the initial velocity. The positive direction of the velocity is in the downward direction.

The use of global variables in the following function (on the left) and script files

function dvdt = ode2a_global(t,v);
% ODE2A_GLOBAL: computes dvdt.
global m_value g_value k_value
dvdt = g_value-(k_value)/(m_value)*sqrt(v);

% script file to solve 1st order ode.
global m_value g_value k_value
m_value=3; g_value=9.81; k_value=6;
v0 = 2;
tspan = [0,30];
[t,v] = ode23(’ode2a_global’,tspan,v0);
plot(t,v)

allow for changing the values of m_value, g_value, and k_value in the script file without changing those variables in the function file. It’s important to note that these global declarations will be available only to the above function and script files.

Save the above function and script file in ode2a_global.m and in ode2a_global_script.m files, respectively. Remember to place them in the working directory of Matlab and try them out. Alternatively, you can also download them from:

and

### 18. Publishing Documents and Opening Web Pages in Matlab

Matlab (version 7, or greater) includes easy to use automatic report generator. It publishes scripts in various formats: HTML/XML, LATEX, (MS Word) DOC, and (MS PowerPoint) PPT formats. The last two formats are available only on PC Windows systems with MS Office installed.

Script my_scrip.m, placed in the working Matlab directory, can be published with the following command:

>> publish(’my_script’)

publish(’my_script’) command creates my_script.html file in html sub-directory of your current working directory. If the script creates plots/graphics, the plots and graphic images will be also included in html sub-directory and automatically linked to my_script.html file.

You can view this newly created my_script.html file in Matlab’s internal browser with the command:

>> open html/my_script.html

or using the command:

>> web(’html/ex1_script.html’)

By the way, the command web(’url’) can be used to display any url in Matlab’s internal browser, for example,

>> web(’www.csun.edu’)

Furthermore, any other browser can also be used to view this file. While viewing my_script.html file in Matlab’s internal or in any external browser, one can also print the file using Print option from the File menu.

The command publish(’my_script’) is a default for publish(’my_script’,’html’) that generates HTML format. The other formats are created by the following commands:

• publish(’my_script’,’latex’)
LATEX source file my_script.tex is generated. If the script creates plots/graphics, the plots and graphic images will be included in .eps (encapsulated Postscript format) files.
• publish(’my_script’,’xml’)
XML file my_script.xml is generated and can be transformed with XSLT or other tools.
• publish(’my_script’,’doc’)
MS Word file my_script.doc is generated. Available only on PC Windows systems with MS Office installed.
• publish(’my_script’,’ppt’)
PowerPoint file my_script.ppt is generated. Available only on PC Windows systems with MS Office installed.

Consider the script from Example 2 of Section 9, where the the Euler method is used to approximate the solution to a given differential equation. The differences between the exact and approximate solutions and the corresponding plots are also provided.

Here is the complete Matlab script:

h=0.1;                             % grid size
x=0:h:2;                           % specify the partition of the interval [0,2]
clear y;                           % clear possibly existing variable y
y(1)=1;                            % initial value: x1=0, thus y(1) corresponds to y(x1)=1
f=inline(’3*exp(-4*x)-2*y’);       % define your function
size(x);                           % size of x to be used in determining the size of vector i
for i=1:20 y(i+1)=y(i)+h*f(x(i),y(i));end   % compute the approximation
Y_exact=2.5*exp(-2*x)-1.5*exp(-4*x);        % the exact solution
error = abs(Y_exact-y)’                     % the differences between exact
% and approximate solutions
plot(x,y,’--’,x,Y_exact)                    % plot of the exact and approximate
% solutions

Save the script in ex1_script.m file (remember to place it in the working directory of Matlab) and try it out. Alternatively, you can also download it from:

Once the Matlab script is running and placed in the current working directory of Matlab, you can publish it using the following Matlab command:

>> publish(’ex1_script’)

ex1_script.html file will be created in html sub-directory of your current working directory. If the script creates plots/graphics (as this is the case here), the plots and graphic images will be also included in html sub-directory and automatically linked to ex1_script.html file.

You can view this newly created ex1_script.html file in Matlab’s internal browser with the command:

>> open html/ex1_script.html

or using the command:

>> web(’html/ex1_script.html’)

And here is the link to ex1_script.html file:

Any other browser can also be used to view this file. While viewing ex1_script.html file in Matlab’s internal or in any external browser, you can also print the file using Print option from the File menu.

It is possible to add further formatting features to ex1_script.m for better readability. Furthermore, one can also use LATEX markup commands to obtain nicely formatted mathematical expressions. And here is an example of such an improved script:

%% A simple example of publishing reports
%% The Euler approximation of the differential equation on [0,2]:
% $$y’+2y=3\exp(-4x),\quad y(0)=1.$$
%% Define the grid size, interval’s partition, and the initial value
h=0.1;                           % grid size
x=0:h:2;                         % specify the partition of the interval [0,2]
clear y;                         % clear possibly existing variable y
y(1)=1;                          % initial value: x1=0, thus y(1) corresponds to y(x1)=1
%% Define the function
f=inline(’3*exp(-4*x)-2*y’);
%% Compute the Euler approximation
size(x);                         % size of x to be used in determining the size of vector i
for i=1:20 y(i+1)=y(i)+h*f(x(i),y(i));end
%% The exact solution
Y_exact=2.5*exp(-2*x)-1.5*exp(-4*x);
%% The differences between exact and approximate solutions
error = abs(Y_exact-y)’
%% Plot the exact and approximate solutions
plot(x,y,’--’,x,Y_exact)

As before, save the script in ex1_script_improved.m file (remember to place it in the working directory of Matlab) and try it out. Alternatively, you can also download it from:

Once the Matlab script is running and placed in the current working directory of Matlab, you can publish it using the following Matlab command:

>> publish(’ex1_script_improved’)

ex1_script_improved.html file will be created in html sub-directory of your current working directory.

You can view this newly created ex1_script_improved.html file in Matlab’s internal browser with the command:

>> open html/ex1_script_improved.html

And here is the link to ex1_script_improved.html file:

After typesetting the LATEX file publish(’ex1_script_improved’,’latex’), its pdf output can be found here:

The last example publishes the script in Example 3 of Section 16. The script uses an additional function file ode3.m that can be downloaded from

and placed in the working Matlab directory.

Here is the corresponding Matlab script:

%% Nonlinear pendulum
%% The differential equation:
% $$\ddot{\theta}+\omega^2\sin\theta=0,\qquad \theta(0)=0,\qquad \dot{\theta}(0)=1.$$
%% The corresponding system of the first order differential equations:
% $$\frac{dz_1}{dt}=z_2,\,\,\frac{dz_2}{dt}=-\omega^2\sin(z_1),\,\, z_1(0)=0,\,\, z_2(0)=1.$$
%% Input initial conditions:
z0=[0,1];
%% Define the interval on which solution is computed:
tspan =[0,20];
%% Solve the system using |ode45| procedure:
[t,z] = ode45(’ode3’,tspan,z0);
%% Extract the positions and velocities:
x=z(:,1); v=z(:,2);
%% Plots of the positions and velocities as functions of time:
% *Note:* _The dashed curve indicates velocities_
plot(t,x,t,v,’--’)
%% Plot of the phase portrait (velocity as the function of position):
figure(2)
plot(x,v)

Save the script in ex2_script.m file (remember to place it in the working directory of Matlab) and try it out. Alternatively, you can also download it from:

Publish it using the Matlab command:

>> publish(’ex2_script’)

and view it with the command:

>> open html/ex2_script.html

And here is the link to ex2_script.html file:

Finally, the output generated by publish(’ex2_script’,’latex’) and converted topdf file can be found here:

### 19. How to Download the Matlab Files Used in this Tutorial

The Matlab files created and used in this tutorial can be downloaded from:

Department of Mathematics, CSUN