MEEN 364 Vijay 1/16/2008
Tutorial on Matlab This tutorial is intended to be a review of some concepts in Matlab about matrix algebra, and ODE’s. 1. Matrix Algebra I.1 Multiplication of matrices
Multiplication of matrices is defined in a way that reflects composition of the underlying linear transformations and allows compact representation of systems of simultaneous linear equations. The matrix product C = AB is defined when the column dimension of A is equal to the row dimension of B, or when one of them is a scalar. If A is m-by- p and B is p-by-n, their product C is m-by-n. The product can actually be defined using MATLAB's for loops, colon notation, and vector dot products. for i = 1:m for j = 1:n C(i, j) = A(i,:)*B(:,j); A(i,:)*B(:,j); end end MATLAB uses a single asterisk to denote matrix multiplication. The next two examples illustrate the fact that matrix multiplication is not commutative; AB is usually not equal to BA.
⎡1 1 1⎤ ⎢ ⎥ A= 1 2 3 ⎢ ⎥ ⎢⎣1 3 6⎥⎦ ⎡8 1 6⎤ ⎢ ⎥ B= 3 5 7 ⎢ ⎥ ⎢⎣ 4 9 2⎥⎦ X = A*B
⎡15 15 15 ⎤ ⎢ ⎥ X = 26 38 26 ⎢ ⎥ ⎢⎣ 41 70 39⎥⎦ Y = B*A
1
MEEN 364 Vijay 1/16/2008
⎡15 28 47⎤ ⎢ ⎥ Y = 15 34 60 ⎢ ⎥ ⎢⎣15 28 43⎥⎦ Rectangular matrix multiplications must satisfy the dimension compatibility conditions. X = A*C
⎡17 19 ⎤ ⎢ ⎥ X= 31 41 ⎢ ⎥ ⎢⎣51 70⎥⎦ Y = C*A gives Error using ==> * Inner matrix dimensions must agree. Anything can be multiplied by a scalar. I.2 Solving Linear Equations
One of the most important problems in technical computing is the solution of simultaneous linear equations. In matrix notation, this problem can be stated as follows. Given two matrices A and B, does there exist a unique matrix X so that AX = B or XA = B? It is instructive to consider a 1-by-1 example. Does the equation have a unique solution? The answer, of course, is yes. The equation has the unique solution x = 3. The solution is easily obtained by division.
The solution is not ordinarily obtained by computing the inverse of 7, that is 7 1 -1 -1 = 0.142857..., and then multiplying 7 by 21. This would be more work and, if 7 is represented to a finite number of digits, less accurate. Similar considerations apply to sets of linear equations with more than one unknown.
MATLAB solves such equations without computing the inverse of the matrix. Although it is not standard mathematical notation, MATLAB uses the division terminology familiar in the scalar case to describe the solution of a general system of simultaneous equations. The two division symbols, slash, /, and backslash, \, are used for the two situations where the unknown matrix appears on the left or right of the coefficient matrix. X = A\B X = B/A
Denotes the solution to the matrix equation AX = B. Denotes the solution to the matrix equation XA = B. 2
MEEN 364 Vijay 1/16/2008 You can think of "dividing" both sides of the equation AX = B or XA = B by A. The coefficient matrix A is always in the "denominator." The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows. The solution X then has the same number of columns as B and its row dimension is equal to the column dimension of A. For X = B/A, the roles of rows and columns are interchanged. In practice, linear equations of the form AX = B occur more frequently than those of the form XA = B. Consequently, backslash is used far more frequently than slash. The remainder of this section concentrates on the backslash operator; the corresponding properties of the slash operator can be inferred from the identity (B/A)' = (A'\B') The coefficient matrix A need not be square. If A is m-by-n, there are three cases. m=n
Square system. Seek an exact solution.
m>n
Over determined system. Find a least squares solution.
m
Underdetermined system. Find a basic solution with at most m nonzero components.
I.3 Eigenvalues and Eigenvalue Decomposition
An eigenvalue and eigenvector of a square matrix A are a scalar satisfy
With the eigenvalues on the diagonal of a diagonal matrix eigenvectors forming the columns of a matrix V , we have
and a vector v that
and the corresponding
If V is nonsingular, this becomes the eigenvalue decomposition
Lets look at an example fro matrix A
⎡ 0 − 6 −1 ⎤ ⎢ ⎥ A= 6 2 − 16 ⎢ ⎥ ⎢⎣− 5 20 − 10⎥⎦ The statement lambda = eig (A) produces a column vector containing the eigenvalues. For this matrix, the eigenvalues are complex.
3
MEEN 364 Vijay 1/16/2008
− 3.0710 ⎡ ⎤ ⎢ ⎥ Lambda = − 2.4645 + 17.0068i ⎢ ⎥ ⎢⎣ − 2.4645 − 17.6008i ⎥⎦ The real part of each of the eigenvalues is negative, so approaches zero as t increases. The nonzero imaginary part of two of the eigenvalues, , contributes the oscillatory component, , to the solution of the differential equation. With two output arguments, eig computes the eigenvectors and stores the eigenvalues in a diagonal matrix. [V, D] = eig (A)
⎡− 0.8326 0.2003 − 0.1394i 0.2003 + 0.1394i ⎤ ⎢ ⎥ V = − 0.3553 − 0.2110 − 0.6447i − 0.2001 + 0.6447i ⎢ ⎥ ⎢⎣ − 0.4248 ⎥⎦ − 0.6930 − 0.6930 0 0 ⎡− 3.0710 ⎤ ⎢ ⎥ − D= 0 2.4645 + 17.6008i 0 ⎢ ⎥ ⎢⎣ 0 − 2.4645 − 17.6008i ⎥⎦ 0
II. ODEs
This section describes the process for solving ODE problems using one of the MATLAB ODE solvers. It uses the Van Der Pol equation as an example. Solving ODE’s of the first order using matlab is pretty straightforward. This section describes how to solve higher order ODE’s. 1 Rewrite the Problem as a System of First-Order ODEs.
ODEs often involve a number of dependent variables, as well as derivatives of order higher than one. To use the MATLAB ODE solvers, you must rewrite such equations as an equivalent system of first-order differential equations of the form
You can write any ordinary differential equation
as a system of first-order equations by making the substitutions
4
MEEN 364 Vijay 1/16/2008
The result is an equivalent system of
first-order ODEs.
For example, you can rewrite the van der Pol equation (second-order)
where is a scalar parameter, by making the substitution system of first-order ODEs is
. The resulting
2 Code the System of First-Order ODEs in MATLAB.
Once you represent the equation as a system of first-order ODEs, you can code it as a function that a MATLAB ODE solver can use. The function must be of the form dydt = odefun(t,y) Although t and y must be the function's first two arguments, the function does not need to use them. The output dydt, a column vector, is the derivative of y. The code below represents the van der Pol system in a MATLAB function, vdp1. The vdp1 function assumes that element vector.
.
and
become elements y(1) and y(2) of a two-
function dydt = vdp1(t,y) dydt = [y(2); (1-y(1)^2)*y(2)-y(1)]; Note that, although vdp1 must accept the arguments t and y, it does not use t in its computations. 3 Apply a Solver to the Problem.
Decide which solver you want to use to solve the problem. Then call the solver and pass it the function you created to describe the first-order system of ODEs, the time interval on which you want to solve the problem, and an initial condition vector. For the van der Pol system, you can use ode45 on time interval [0 20] with initial values y (1) = 2 and y (2) = 0.
5
MEEN 364 Vijay 1/16/2008 [t, y] = ode45 (@vdp1,[0 20],[2; 0]); This example uses @ to pass vdp1 as a function handle to ode45. The resulting output is a column vector of time point’s t and a solution array y. Each row in y corresponds to a time returned in the corresponding row of t. The first column of y corresponds to the second column to
, and
.
4 View the Solver Output.
You can simply use the plot command to view the solver output. plot (t, y(:,1),'-',t,y(:,2),'--') title ('Solution of van der Pol Equation, \mu = 1'); xlabel('time t'); ylabel('solution y'); legend('y_1','y_2')
Let us now consider solving more complicated ODE’s using Matlab. These are differential equations whose coefficients are not constant. They generally represent nonlinear systems. Example A Bar supported by a wire and a horizontal plane
6
MEEN 364 Vijay 1/16/2008 Consider the system shown below. It consists of a bar BC of length l and mass m that is supported by a mass less wire from point A to end B. End C of the bar rests on a horizontal plane that lies a distance d below point A. The wire from point A to B has length L and. Choose ‘l’, ‘L’, ‘m’ and ‘d’ as 1m, 1m, 5 Kg and 0.5m respectively. The initial values of ‘ θ’and ‘φ’ are taken to be 0.4857 radians and 0.5233 radians respectively. Let us plot the variation of ‘ θ’and ‘φ’ with respect to time. A L
φ B d
C
θ
Solution The model for the above system is
⎡ ml 2 ⎤ −l l 0 ⎧ ⎫ .. + 0 sin( θ φ ) cos θ 2 . ⎢ 12 ⎥ ⎧ ⎫ ⎪ ⎪ ml 2 2 ⎢ ml ⎥ ⎪ θ .. ⎪ ⎪ − w + sin θ θ ⎪ 2 ⎢ cosθ ⎪ − sin φ − 1 ⎥ ⎪ φ ⎪ = ⎪ 0 2 2 ⎢2 ⎥ ⎨ ⎬ ⎨mL cos φ φ . − ml cosθ θ . ⎬ ⎪ ⎢ ml sin θ − mL sin φ ⎥ ⎪T c ⎪ ⎪ 2 cos φ 0 ⎪ ⎢ 2 ⎥ ⎪⎩ N ⎪⎭ ⎪ . 2 . 2 ⎪ ⎢ l cosθ L cos φ 0 0 ⎥⎦ ⎩ L sin φ φ + l sin θ θ ⎪⎭ ⎣ Where T c represents the tension in the chord and N the normal reaction exerted on the bar at the ground. Other parameters are clear from the diagram. .
First we solve for θ and
θ. Then we substitute them into the following constraint .
equations to obtain the values of φ and φ. The constraint equations are the geometrical relations, which can be derived from the above figure to be .
.
L cos φ φ + l cosθ θ = 0 d = L sin φ + l sin θ
The above model involves second derivatives of the variables and must be converted to an equation involving first derivatives so that MATLAB can be used to solve the differential equation. So with
7
MEEN 364 Vijay 1/16/2008
θ = y (1); .
θ = y (2); φ = y (3); .
φ = y (4);
The model becomes
⎡1 ⎢ ⎢0 ⎢0 ⎢ ⎢0 ⎢ ⎢0 ⎢ ⎢⎣0
0 2 ml
ml
12 0
2 ml
0 0 1
cosθ 0
sin θ 2 l cosθ
0 0
y (2) ⎫ ⎤ θ ⎧⎪ ⎪ ⎥⎧ . ⎫ ⎪ 0 −l l ⎪ 0 sin(θ + φ ) cosθ ⎥ ⎪ θ ⎪ 2 2 ⎪ y (4) ⎥ ⎪ ⎪ ⎪⎪ 2 0 0 0 ⎪⎪ . ⎥ ⎪⎪ φ ⎪⎪ ⎪ − + ml . w sin θ θ = ⎬ − sin φ − 1 ⎥ ⎨ φ ⎬ ⎨ 0 2 ⎪ ⎪ ⎪ ⎪ ⎥ . 2 . 2 ml ⎪ ⎥⎪ ⎪ ⎪ − mL sin φ cos φ 0 ⎥ ⎪T c ⎪ ⎪mL cos φ φ − 2 cosθ θ ⎪ ⎪ N ⎪ ⎪ ⎪ . 2 . 2 L cos φ 0 0 ⎥⎦ ⎩ ⎭ ⎪ L ⎩ sin φ φ + l sin θ θ ⎪⎭
0
0
0
[T, Y] = ODEsolver ('F', TSPAN, Y0) with TSPAN = [T0 TFINAL] integrates the system of differential equations y' = F (t,y) from time T0 to TFINAL with initial conditions Y0. 'F' is a string containing the name of an ODE file. Function F(T,Y) must return a column vector. Each row in solution array Y corresponds to a time returned in column vector T. To obtain solutions at specific times T0, T1, ..., TFINAL (all increasing or all decreasing), use TSPAN = [T0 T1 ... TFINAL]. Some ODEsolvers can solve problems M(t,y)*y' = F(t,y) with a mass matrix that is non-singular. ( eg.ODE113)
8
MEEN 364 Vijay 1/16/2008
MATLAB Code The code, compatible with all version of Matlab after version 5, is included below. To store the right hand side vector of the original model, a separate function file must be written as shown below. Note that the name of the function is ‘FF1’, so this file must be saved as ‘FF1.m’. The following function contains the right hand side of the differential equation of the form M(t,y)*y'=F(t,y) i.e. it contains F(t,y).it is also stored in a separate file named, FF1.m. function yp=FF1(t,y) l=1; L=1; m=5; g=9.81; w=m*g; yp=zeros(6,1); yp(1)=y(2); yp(2)=0; yp(3)=y(4); yp(4)=-w+(m*l/2)*sin(y(1))*(y(2)^2); yp(5)=m*L*cos(y(3))*(y(4)^2)-(m*l/2)*cos(y(1))*(y(2)^2); yp(6)=L*sin(y(3))*(y(4)^2)+l*sin(y(1))*(y(2)^2); Similarly, to store the coefficient matrix, a separate function file is written which is stored as ‘MM1.m’.
The following function contains the mass separately stored in a file named, MM1.m
matrix.
It
is
function n = MM1(t,y) l=1; L=1; m=5; g=9.81; w=m*g; n1=[1 0 0 0 0 0]; n2=[0 (m*l^2)/12 0 0 (-l/2)*sin(y(1)+y(3)) (l/2)*cos(y(1))]; n3=[0 0 1 0 0 0]; n4=[0 (m*l/2)*cos(y(1)) 0 0 -sin(y(3)) -1]; n5=[0 (m*l/2)*sin(y(1)) 0 -m*L*sin(y(3)) cos(y(3)) 0]; n6=[0 l*cos(y(1)) 0 L*cos(y(3)) 0 0]; n=[n1;n2;n3;n4;n5;n6];
9
MEEN 364 Vijay 1/16/2008 To obtain the response, the main file should call the function ‘indmot_ode1.m’, which has the switch/case programming which in turn calls the corresponding functions depending on the value of the flag. For the main file to recognize the coefficient matrix, the MATLAB command ODESET is used to set the mass to ‘M (t, y)’. l=1; L=1; d=0.5; tspan=[0 10] options=odeset('mass',@MM1) y0=[0.5233;0;1.0467;0;0;0] [t,y]=ode113(@FF1,tspan,y0,options) %plot(t,y(:,1)) %grid theta=y(:,1); thetadot=y(:,2); % solving the constraint equation to get the value of phi. for i = 1:size(theta,1) phi(i)=asin((d-l*sin(theta(i)))/L); end % solving the constraint equation to get the phidot. for i = 1:size(theta,1) phidot(i)=(l*(thetadot(i))*cos(theta(i)))/(L*cos(phi(i))); end
value
of
t1=t'; plot(t1,phi) .
To plot the value of θ with respect to time, change the variable in line 8 of the main code from ‘y(:,1)’ to ‘y(:,2)’. This step is required because we initially assigned y(1) to ‘ θ’ and .
.
y(2) to θ . Once the values of ‘ θ’ and θ are known, they are substituted in the constraint .
equation to get the values of ‘ φ’ and φ . The plots are attached below.
10
MEEN 364 Vijay 1/16/2008
11
MEEN 364 Vijay 1/16/2008
12
MEEN 364 Vijay 1/16/2008 III. FFT One-dimensional fast Fourier transform. Y = fft (X) returns the discrete Fourier transform (DFT) of vector X, computed with a fast Fourier transform (FFT) algorithm = fft (X, n) returns the n-point DFT. If the length of X is less than n, X is padded with trailing zeros to length n. If the length of X is greater than n, the sequence X is truncated. When X is a matrix, the length of the columns are adjusted in
the same manner. Examples A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal. Consider data sampled at 1000 Hz. Form a signal containing 50 Hz and 120 Hz and corrupt it with some zero-mean random noise: t = 0:0.001:0.6; x = sin(2*pi*50*t)+sin(2*pi*120*t); y = x + 2*randn(size(t)); plot(y(1:50)) title('Signal Corrupted with Zero-Mean Random Noise') xlabel('time (seconds)')
It is difficult to identify the frequency components by looking at the original signal. Converting to the frequency domain, the discrete Fourier transform of the noisy signal y is found by taking the 512-point fast Fourier transform (FFT): Y = fft(y,512); The power spectrum, a measurement of the power at various frequencies, is Pyy = Y.* conj(Y) / 512;
13
MEEN 364 Vijay 1/16/2008 Graph the first 257 points (the other 255 points are redundant) on a meaningful frequency axis. f = 1000*(0:256)/512; plot(f,Pyy(1:257)) title('Frequency content of y') xlabel('frequency (Hz)')
14
