[PPT] - Introduction to Matlab with mathematical and engineering PowerPoint Presentation

SLIDE 1

Introduction to Matlab with mathematical and engineering applications

Hanumant Singh Shekhawat

Indian Institute of Technology Guwahati India

April 11,2015

(IITG, India) April 11,2015 1 / 30

SLIDE 2

Objective

Learning basic Matlab with some elementary engineering and mathematical applications

(IITG, India) April 11,2015 2 / 30

SLIDE 3

Outline

1 Round off error 2 Matrices 3 Basic Numerical Mathematics 4 Ordinary differential equations

(IITG, India) April 11,2015 3 / 30

SLIDE 4

Round off error

Define variable x = 1 − 1/3 − 2/3.
Theoretically x = 0
check using x == 0
What do you get ? and why ?

(IITG, India) April 11,2015 4 / 30

SLIDE 5

Matrices

Different data types

(IITG, India) April 11,2015 5 / 30

SLIDE 6

Matrices

Linear equations

Consider the following linear equation

3x1 + 4x2 = 1 x2 = 2 x1 + x3 = 0

The above is equivalent to

Ax = B where 3 × 3 matrix A =   3 4 1 1 1   , B =

1

2 T and x =

x1

x2 x3 T

x = A−1B
The matrix notation is useful in solving large number of

equations.

(IITG, India) April 11,2015 6 / 30

SLIDE 7

Matrices

Solving linear equation

Define a 3 × 3 matrix A =

  3 4 1 1 1   in Matlab as A =

3

4 0; 1 0; 1 1

Solve for x in Ax = B where vector B =
1

2 T (Matlab: B =

1

2 ′ ).

x = A−1B (in Matlab: x = inv(A) ∗ B)

(IITG, India) April 11,2015 7 / 30

SLIDE 8

Matrices

Solving linear equation in a faster way

Define a 3 × 3 matrix A =

  3 4 1 1 1   in Matlab as A =

3

4 0; 1 0; 1 1

Solve for x in Ax = B where vector B =
1

2 T (Matlab: B =

1

2 ′ )

Matlab: x = A\B
The time difference is visible in case of a large size matrix.
N = 1500;
A = rand(N, N);
B = rand(N, 1);
tic; inv(A) ∗ B; toc
tic; A\B; toc

(IITG, India) April 11,2015 8 / 30

SLIDE 9

Matrices

Least square solution of a linear equation

Define a 3 × 2 matrix A =

  3 4 1 1   in Matlab as A =

3

4; 1; 1

Solve for x in Ax = B where vector B =
1

2 T (Matlab: B =

1

2 ′ )

Theory: No solution
Matlab: x = A\B. You will get an answer !!
Matlab gives the least square solution
x minimizes y − B2 = 3

1(yi − Bi)2 where y = Ax

(IITG, India) April 11,2015 9 / 30

SLIDE 10

Matrices

Many solutions of a linear equation

Define a 2 × 3 matrix A =

3 4 1 1

in Matlab as

A =

3

4 0; 1 1

Solve for x in Ax = B where vector B =
1

2 T (Matlab: B =

1

2 ′ ).

Matlab: x = A\B
Obtain y = null(A);
Check x + λy is also a solution of Ax = B.
Which solution to use in case of many solutions?

(IITG, India) April 11,2015 10 / 30

SLIDE 11

Matrices

Exercise

Let continuous function f : [−1 1] → R be unknown. However, we have r measurements of it, say at x1, · · · , xr. Find an approximation fr (it is a continuous function) such that f − fr minimized at the measurement point i.e. minimize r

i=1(f (xi) − fr(xi))2.

Out of many methods we choose the 3rd order polynomial

approximation i.e. fr(x) = a0 + a1x + a2x2 + a3x3

let the measurement is given by

xi = i ∗ 0.01; f (xi) = cos(xi) + sin(xi)3; for i = 0, 1, · · · 100.

In reality, we just have data points f (xi) instead of an equation.
Use f (xi) = fr(xi) to generate 101 equation. Use matrices here
therwise you will waste a lot of time and paper.
Find the least square solution a =
a0

a1 a2 a3

Plot f (xi) and fr(xi).

(IITG, India) April 11,2015 11 / 30

SLIDE 12

Matrices

Advance Exercises

Repeat the last exercise with the following approximations:

Higher order polynomial approximation ? Do you observe

singularity ?

Legendre polynomials approximation

fr(x) = a0 + a1x + a2 3x2 − 1 2 + a3 5x3 − 3x 2 Do you observe any advantage over polynomials?

Fourier approximation

fr(x) = a0 + a1cos(2πx) + a2cos(4πx) + a3cos(6πx)

Which one gives a better approximation ?

(IITG, India) April 11,2015 12 / 30

SLIDE 13

Matrices

Eigenvalue decomposition

Define a 3 × 3 matrix A =

  10 12 1 12 17 1 1  

Calculate A3.
Now represent A as A = MDM−1 (the eigenvalue

decomposition) where D is a diagonal matrix (non-zero entries

n the diagonal only). Here use [M, D] = eig(A) in Matlab.
A3 = MD3M−1. Which one is easy to compute?
Check the above in Matlab.
Take any colom of M (say x) and calculate y = Ax.
Divide first entry of y with first entry of x. Repeat for second

and third.

Did you get same result ?
Coloms of M are eigenvectors and diagonal entries of D are

called eigenvalues.

(IITG, India) April 11,2015 13 / 30

SLIDE 14

Matrices

Singular value decomposition (SVD)

Matrix A can be written as

A = UDV T where D has entries (known as singular values) only along the diagonal and UUT = VV T = I. I is the identity matrix.

Matlab
A= rand(3,2);
[U, D, V ] = svd(A);
UUT = I and VV T = I (size of I may change)
Get size of all matrices.
Get eigenvalue decomposition of AAT and ATA.
Do you see U and V there ?
Any relationship between D here and the eigenvalues of AAT

and ATA.

(IITG, India) April 11,2015 14 / 30

SLIDE 15

Matrices

Image

Download the image from http://img1.jurko.net/wall/uploads/wallpaper_13409.jpg and save as ’end.jpeg’.

(IITG, India) April 11,2015 15 / 30

SLIDE 16

Matrices

Image

Download the image from http://img1.jurko.net/wall/uploads/wallpaper_13409.jpg and save as ’end.jpeg’. Use the following commands to read the image in Matlab as a matrix

X = imread(’end.jpeg’,’jpg’);
X = rgb2gray(X);
image(X)
colormap(gray(256))

X is a matrix of size 600 × 800

(IITG, India) April 11,2015 16 / 30

SLIDE 17

Matrices

Image compression

Compression means reducing the memory requirement to store

the data.

Image is a matrix M = UDV T (SVD).
One method of image compression is to replace smaller singular

values with zero.

˜

M = U ˜ DV T where ˜ D is the same matrix as D except that it contains only the r largest singular values (the other singular values are replaced by zero).

M − ˜

MF = D − ˜

DF. Here . just means root mean square
f all entries of the matrix.
Try this for a random matrix A = rand(4, 4).
Why this is a compression ?

(IITG, India) April 11,2015 17 / 30

SLIDE 18

Matrices

Exercise (Image compression)

Download the image from http:

//img1.jurko.net/wall/uploads/wallpaper_13409.jpg and save as ’end.jpeg’.

Matlab: X = imread(’end.jpeg’,’jpg’);
Matlab: X = rgb2gray(X);
Matlab: image(X)
Matlab: colormap(gray(256))
Compute the SVD of X, using the command: [U,D,V] = svd(X);
Retain only k = 50 singular values and make the rest zero.
Try with k equal to 100, 150, 200.
What do you observe ?
Calculate compression ratio ?
Plot the singular values (divided by the largest) on the semi-log

scale.

(IITG, India) April 11,2015 18 / 30

SLIDE 19

Matrices

Advance Exercises

How many singular values are sufficient for an approximation ?
Divide the image in four equal pieces and do an approximation
f all four images. Now join them together. Is this method gives

you good result ?

(IITG, India) April 11,2015 19 / 30

SLIDE 20

Basic Numerical Mathematics

Fixed points

Define variable x0 = 100.
Find x1 using iterations xn+1 = cos(xn)
Repeat 30 times to get x30.
Did you get 0.73908 · · · ? why ?
Repeat with x0 = 0
Repeat with x0 = 9e6
0.73908 · · · is a fixed point of iterations xn+1 = cos(xn) and is a

solution of x = cos(x). Why ?

No dynamics at the fixed point !
Can you write a program without worrying of number of

repetitions to get a fixed point?

In general, all solutions of x = f (x) are called fixed points.

(IITG, India) April 11,2015 20 / 30

SLIDE 21

Basic Numerical Mathematics

Fixed points, contd

Define variable x0 = 0.673.
Find x1 using iterations xn+1 = x2

n

Repeat 20 times to get x20.
Did you get 0 ? why ?
Repeat with x = 1
Repeat with x = 1.000001 and x = 0.999999
Slight change in 1 and we are away from 1.
1 is an unstable fixed point of iterations xn+1 = x2

n.

Repeat with x = 1e − 6 and x = −1e − 6
Slight change in 0 and we are on 0 again.
0 is a stable fixed point of iterations xn+1 = x2

n.

(IITG, India) April 11,2015 21 / 30

SLIDE 22

Basic Numerical Mathematics

Discrete time system

A discrete time system (zero-input) is defined by an iteration or

difference equation xk+1 = f (xk), k = 1, 2, · · · where f : Rn → Rn.

x is called fixed point of discrete time system xk+1 = f (xk) if

x = f (x)

Loosely speaking stability of a fixed point x means that a

variation in x will lead iteration xk+1 = f (xk) to x.

(IITG, India) April 11,2015 22 / 30

SLIDE 23

Basic Numerical Mathematics

Stability of discrete time system

A discrete time system (zero-input) is defined by an iteration or

difference equation xk+1 = f (xk), k = 1, 2, · · · where f : Rn → Rn.

Calculate Jacobin/Gradient J of f (x) = [f1(x), · · · , fN(x)]T at

the fixed point x∗ = [x∗

1, · · · , x∗ N]T as

J(x∗) =   

∂f1 ∂x1(x∗)

· · ·

∂f1 ∂xN (x∗)

. . . ... . . .

∂fN ∂x1 (x∗)

· · ·

∂fN ∂xN (x∗)

  

x∗ is stable if all eigenvalues λ of J(x∗) satisfy |λ| < 1
x∗ is unstable if all eigenvalues λ of J(x∗) satisfy |λ| > 1
No conclusions on stability if |λ| = 1 for some eigenvalue λ of

J(x∗).

(IITG, India) April 11,2015 23 / 30

SLIDE 24

Basic Numerical Mathematics

Examples of stability

Jacobin/Gradient J of f (x) = [f1(x), · · · , fN(x)]T at the fixed

point x∗ = [x∗

1, · · · , x∗ N]T is given by

J(x∗) =   

∂f1 ∂x1(x∗)

· · ·

∂f1 ∂xN (x∗)

. . . ... . . .

∂fN ∂x1 (x∗)

· · ·

∂fN ∂xN (x∗)

  

Now you can proof stability in all previous example.
xk+1 = cos(xk). Here J(x∗) = − sin(x∗).
For x∗ = 0.73908 we have that |λ| = | − sin(0.73908)| < 1 (as

0.73908 < π/2). This implies x∗ = 0.73908 is a stable fixed point.

xk+1 = x2
k. Here J(x∗) = 2x∗
For x∗ = 0 we have that |λ| = 0 < 1. This implies x∗ = 0 is a

stable fixed point.

For x∗ = 1 we have that |λ| = 2 > 1. This implies x∗ = 1 is an

unstable fixed point.

(IITG, India) April 11,2015 24 / 30

SLIDE 25

Basic Numerical Mathematics

Exercise

Jacobin/Gradient J(x∗) of f (x) = [f1(x), · · · , fN(x)]T at the

fixed point x∗ = [x∗

1, · · · , x∗ N]T is given by

J(x∗) =   

∂f1 ∂x1(x∗)

· · ·

∂f1 ∂xN (x∗)

. . . ... . . .

∂fN ∂x1 (x∗)

· · ·

∂fN ∂xN (x∗)

  

let f (x) = M(B − g(x)) where M =

2/3 1/3 1/3 2/3

,

B =

−1

1 T and g(x) = 1

9

e−x1

e−x2T

Find a fixed point of iteration xk+1 = f (xk) and show that the

fixed point is stable. Use x0 =

1

1 T.

Use J(x∗) = 1

9M

−e−x∗

1

−e−x∗

2

. Verify !

(IITG, India) April 11,2015 25 / 30

SLIDE 26

Ordinary differential equations

Consider ˙

x = f (x).

There is no dynamics at fixed point x∗, so x∗(t) do not very

with time. This means ˙ x∗ = f (x∗) = 0.

Fixed points are also called equilibrium points, constant

solutions, working points, steady solutions, stagnation points.

For a small time step h, the numerical solution can be obtained

using different approximation ˙ x ≈ x(t + h) − x(t) h ˙ x ≈ x(t) − x(t − h) h ˙ x ≈ x(t + h) − x(t − h) 2h

(IITG, India) April 11,2015 26 / 30

SLIDE 27

Ordinary differential equations

Solving differential equation dx

dt = rx(1 − x α), x(0) = 1 in Matlab.

>> t0=0; te=10; % Initial and the final time for simulation >> x0=1; % Initial value >> r=1; alpha=10; % parameter of the differential equation >> [t,x] = ode45(@mydiff,[t0 te],x0,[],r,alpha); % [] for some

ptions, see Matlab help

>> plot(t,x); >> title(’Solution’) mydiff is the function defined in file mydiff.m function dxdt = mydiff(t,x,r,alpha); dxdt=rx(1-x/alpha); Try with different initial values.

(IITG, India) April 11,2015 27 / 30

SLIDE 28

Ordinary differential equations

Exercise (Van der Pol oscillator)

Solve for z(t) where ¨ z + r(z2 − 1)˙ z + z = 0 and r = 1.5.

define y = ˙

z then ˙ z ˙ y

=
y

−r(z2 − 1)y + z

x =
z

y T.

function dxdt=mydiffVan(t,x,r)

z=x(1); y=x(2); dxdt=[ y;

r ∗ (zˆ2 - 1) ∗ y - z ];
Use [t,x] = ode45(@mydiffVan,[t0 te],x0,[],r);
plot(t,x)
Use different initial values (start with x =
0.5

T)

(IITG, India) April 11,2015 28 / 30

SLIDE 29

Applications

All this is applied in

Linear algebra
Chaos theory
Control systems
Signal processing
Tensors
Data compression
And many more .....

(IITG, India) April 11,2015 29 / 30

SLIDE 30

Main References

1 Strang, Gilbert. Linear algebra and its applications. Belmont,

CA: Thomson, Brooks/Cole, 2006.

Good for Matrices and the related stuff

2 Strogatz, Steven. Nonlinear dynamics and chaos: with

applications to physics, biology, chemistry, and engineering. Westview press, 2014.

Good for discrete systems, differential equations, fixed point etc.
There are various definitions of stability. Please see this book

for details.

3 Epperson, James. An introduction to numerical methods and

analysis. John Wiley & Sons, 2013.

(IITG, India) April 11,2015 30 / 30

SLIDE 31

Contact

Please contact me at h.s.shekhawat [AT] iitg.ernet.in for any error in the Matlab code and in the presentation.

(IITG, India) April 11,2015 31 / 30