>> picocrse %%%%%%% Introduction MATLAB is a specialized software system for numerical computations with matrices. We will see that very sophisticated functions and computation routines can be built from the basic matrix operations within MATLAB. This 2-hour "Pico"-Course consists of a series of brief lessons intended to introduce the syntax of MATLAB and especially to present some commonly used functions and techniques. In addition, some very specific topics useful for 6.003 are discussed in more detail. We *strongly advise* that the student becomes accustomed to plotting capabilities, Help features, and defining one's own specialized functions using M-files in MATLAB. With this goal in mind, we close our presentation by using the power of MATLAB to analyze the step response for the parallel RLC circuit. The overall presentation consists of five lessons: Lesson 1: Generating Matrices Lesson 2: Basic Operations: Matrix and Element-wise Lesson 3: Working with Complex Numbers and Polynomials Lesson 4: Plotting and Printing Lesson 5: Illustrative Example -- Parallel RLC Circuit (Using MATLAB Help features and defining one's own Scripts/Functions using M-files) Select a lesson [1,2,3,4,5]: 1 %%%%%%% Lesson 1: Generating Matrices % % 1) Basic syntax for defining matrices, vectors, and scalars % 2) Some special functions: "colon" operator,"linspace","eye","zeros","ones" % 3) Extracting elements, vectors, and blocks from existing matrices % % At every keypress, one or more MATLAB commands is executed and an % explanation will follow the format below: (Actual MATLAB syntax) %(brief explanation of the command) (MATLAB's response to the issued command) pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%(ready for the next keypress) A = [1 7; -4 8] %2x2 matrix A = 1 7 -4 8 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% b = [9; 3] %2x1 column vector b = 9 3 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c1 = [9 3] %1x2 row vector c1 = 9 3 c2 = b' %1x2 row vector -- via transpose c2 = 9 3 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% alpha = 1 %scalar alpha = 1 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% H = [A b] %2x3 matrix -- via block form (DIMENSIONS!) H = 1 7 9 -4 8 3 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% keyboard %To continue, type: return K>> K>> return return t1 = [0:4] %1x5 row vector -- via "colon" t1 = 0 1 2 3 4 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% t2 = linspace(0,4,5) %1x5 row vector -- via "linspace" t2 = 0 1 2 3 4 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% t3 = [-1:0.2:0]' %6x1 column vector -- via "colon" t3 = -1.0000 -0.8000 -0.6000 -0.4000 -0.2000 0 t4 = linspace(-1,0,6)' %6x1 column vector -- via linspace t4 = -1.0000 -0.8000 -0.6000 -0.4000 -0.2000 0 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% t5 = [0:-0.2:-1] %And again, in descending order t5 = 0 -0.2000 -0.4000 -0.6000 -0.8000 -1.0000 t6 = linspace(0,-1,6) t6 = 0 -0.2000 -0.4000 -0.6000 -0.8000 -1.0000 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% I1 = [1 0 0; 0 1 0 ; ... I1 = 1 0 0 0 1 0 0 0 1 I2 = eye(3) %3x3 Identity Matrix -- via "eye" I2 = 1 0 0 0 1 0 0 0 1 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Z1 = [0 0 0; 0 0 0] %2x3 zero matrix Z1 = 0 0 0 0 0 0 Z2 = zeros(2,3) %2x3 zero matrix -- via "zeros" Z2 = 0 0 0 0 0 0 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% O1 = [1 1; 1 1; 1 1] %3x2 ones matrix O1 = 1 1 1 1 1 1 O2 = ones(3,2) %3x2 ones matrix -- via "ones" O2 = 1 1 1 1 1 1 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% txt = 'string' %Can also generate matrices/vectors of characters txt = string txt = ['s' 't' 'r'; 'i', 'n', 'g'] %Character-by-character example, txt = str ing pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [numrows,numcols] = size(H) %Getting dimensions of existing variables numrows = 2 numcols = 3 numrows = size(H,1) %Getting number of rows only numrows = 2 numcols = size(H,2) %Getting number of columns only numcols = 3 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% A %Extracting elements from existing matrices A = 1 7 -4 8 alpha = A(1,1) %scalar, upper-left element of A alpha = 1 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% H %can use the "colon" operator H = 1 7 9 -4 8 3 A = H(1:2,1:2) %2x2 matrix block, upper-left block of A A = 1 7 -4 8 b = H(1:2,3) %2x1 column vector, top of 3rd column of A b = 9 3 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% A = H(:,1:2) %2x2 matrix block, all rows of left twp columns of H A = 1 7 -4 8 b = H(:,3) %2x2 matrix block, all rows of third column of H b = 9 3 F = H(:,[1 3]) %2x2 matrix block, 1st & 3rd columns of A F = 1 9 -4 3 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% H %can use the "size" function H = 1 7 9 -4 8 3 A = H(1:size(H,1),1:2) %2x2 matrix block, all rows of left twp columns of H A = 1 7 -4 8 b = H(1:size(H,1),size(H,2)) %2x2 matrix block, all rows of third column of H b = 9 3 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% echo off >> picocrse %%%%%%% Introduction MATLAB is a specialized software system for numerical computations with matrices. We will see that very sophisticated functions and computation routines can be built from the basic matrix operations within MATLAB. This 2-hour "Pico"-Course consists of a series of brief lessons intended to introduce the syntax of MATLAB and especially to present some commonly used functions and techniques. In addition, some very specific topics useful for 6.003 are discussed in more detail. We *strongly advise* that the student becomes accustomed to plotting capabilities, Help features, and defining one's own specialized functions using M-files in MATLAB. With this goal in mind, we close our presentation by using the power of MATLAB to analyze the step response for the parallel RLC circuit. The overall presentation consists of five lessons: Lesson 1: Generating Matrices Lesson 2: Basic Operations: Matrix and Element-wise Lesson 3: Working with Complex Numbers and Polynomials Lesson 4: Plotting and Printing Lesson 5: Illustrative Example -- Parallel RLC Circuit (Using MATLAB Help features and defining one's own Scripts/Functions using M-files) Select a lesson [1,2,3,4,5]: 2 %%%%%%% Lesson 2: Basic Operations: Matrix and Element-wise % % 1) Matrix arithmetic and special functions for matrices % addition, subtraction, multiplication, inverse ("division"), exponent % 2) Element-wise arithmetic and special elememt-wise functions % multiplication, division, exponent, trig and log/exp functions % % At every keypress, one or more MATLAB commands is executed and an % explanation will follow the format below: (Actual MATLAB syntax) %(brief explanation of the command) (MATLAB's response to the issued command) pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%(ready for the next keypress) A1 = [1 7; -4 8]; %Generate "workspace" for lesson A2 = [-1 0; 3 2]; H = [-2 5 -7; 3 1 2]; b2 = [9; 3]; b3 = [b2; 1]; c = [2 0]; alpha = 1; pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% B1 = A1 + A2 %Basic Matrix Operations -- addition B1 = 0 7 -1 10 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% B2 = A2 - A1 %Basic Matrix Operations -- subtraction B2 = -2 -7 7 -6 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C1 = A1*H %Basic Matrix Calculations -- multiplication C1 = 19 12 7 32 -12 44 %Pay attention to dimensions!!!! keyboard %To continue lesson, type: return K>> K>> K>> K>> return return C2 = H'*A1 %Inner dimensions agree now...phew! 8-) C2 = -14 10 1 43 -15 -33 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% v1 = A1*b2 %Example with vector v1 = 30 -12 v2 = alpha*A1*b2 %Example with scalar and vector v2 = 30 -12 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% A2inv = inv(A2) %Basic Matrix Operations -- inverse A2inv = -1.0000 0 1.5000 0.5000 x = inv(A1)*b2 %Solve linear equation Ax = b for variable x x = 1.4167 1.0833 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Basic Matrix Operations -- ?????division?????? x = A1\b2 %Equivalent to x = inv(A1)*b2 x = 1.4167 1.0833 y = A1\A2 %Equivalent to x = inv(A1)*A2 y = -0.8056 -0.3889 -0.0278 0.0556 y = A1/A2 %Equivalent to x = A1*inv(A2) y = 9.5000 3.5000 16.0000 4.0000 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Asquared = A1^2 %Basic Matric Operations -- exponent Asquared = -27 63 -36 36 Asquared = A1*A1 %Equivalently, Asquared = -27 63 -36 36 Aroot = A1^0.5 %Matrix square root -- via exponent operation Aroot = 1.5275 1.5275 + 0.0000i -0.8729 - 0.0000i 3.0551 + 0.0000i Aroot = sqrtm(A1) %Matrix square root -- via "sqrtm" function Aroot = 1.5275 1.5275 -0.8729 3.0551 %Notice complex numbers! (See lesson 3) pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Can combine operations in many ways T = alpha*(inv(A1)*v1+c')'*H %using standard "order-of-operations" T = -13.0000 58.0000 -71.0000 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear %Clear "workspace" -- take all variables out of memory pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% A1 = [1 7; -4 8]; %regenerate "workspace" and continue lesson A2 = [-1 0; 3 2]; H = [-2 5 -7; 3 1 2]; b2 = [9; 3]; b3 = [b2; 1]; c = [2 0]; alpha = 1; pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C1 = A1.*A2 %Basic Element-wise Operations -- multiplication C1 = -1 0 -12 16 v1 = b2.*c' %Example with vectors v1 = 18 0 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Basic Element-wise Operations -- division M1 = A1./A2 %Equivalent to each element m = a1 divided by a2 Warning: Divide by zero. > In /afs/athena.mit.edu/course/6/6.003/picocourse/s97less2.m at line 64 In /afs/athena.mit.edu/course/6/6.003/picocourse/picocrse.m at line 15 M1 = -1.0000 Inf -1.3333 4.0000 M1 = A1.\A2 %Equivalent to each element m = a2 divided by a1 M1 = -1.0000 0 -0.7500 0.2500 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Basic Element-wise Operations -- exponent Melsqua = A1.^2 %Square each element Melsqua = 1 49 16 64 Melroot = A1.^0.5 %Take square root of each element -- via exponent Melroot = 1.0000 2.6458 0.0000 + 2.0000i 2.8284 Melroot = sqrt(A1) %Take square root of each element -- via "sqrt" Melroot = 1.0000 2.6458 0 + 2.0000i 2.8284 Wow = 10.^A1 %Elements are 10 raised to the power a1 Wow = 1.0e+08 * 0.0000 0.1000 0.0000 1.0000 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Note that Wow's elements can be displayed meaningfully format short e %Change format of MATLAB output Wow = 10.^A1 %Elements are 10 raised to the power a1 Wow = 1.0000e+01 1.0000e+07 1.0000e-04 1.0000e+08 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% format short %Change format back to MATLAB default setting %Many familiar element-wise functions yrad = [-pi/2 0 pi/2] %Note predefined variable pi yrad = -1.5708 0 1.5708 ydeg = degrees(yrad) %Readily convert from radians to degrees ydeg = -90 0 90 yrad = radians(ydeg) %And vice-versa...use radians for MATLAB trig functions yrad = -1.5708 0 1.5708 Sinofy = sin(yrad) %Trigonometric functions and their inverses Sinofy = -1 0 1 ydeg0 = degrees(atan(tan(radians(ydeg)))) %Nest functions easily! ydeg0 = -90 0 90 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %More element-wise functions y = [-1 0 1] %RECALL: Could have used y = [-1:1] = linspace(-1,1,3) y = -1 0 1 Expofy = exp(y) %Logarithmic/exponential functions Expofy = 0.3679 1.0000 2.7183 Logofy = log(Expofy) %Natural logarithm (base e) Logofy = -1 0 1 What = log10(10.^y) %Logarithm base 10 What = -1.0000 0 1.0000 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% echo off >> picocrse %%%%%%% Introduction MATLAB is a specialized software system for numerical computations with matrices. We will see that very sophisticated functions and computation routines can be built from the basic matrix operations within MATLAB. This 2-hour "Pico"-Course consists of a series of brief lessons intended to introduce the syntax of MATLAB and especially to present some commonly used functions and techniques. In addition, some very specific topics useful for 6.003 are discussed in more detail. We *strongly advise* that the student becomes accustomed to plotting capabilities, Help features, and defining one's own specialized functions using M-files in MATLAB. With this goal in mind, we close our presentation by using the power of MATLAB to analyze the step response for the parallel RLC circuit. The overall presentation consists of five lessons: Lesson 1: Generating Matrices Lesson 2: Basic Operations: Matrix and Element-wise Lesson 3: Working with Complex Numbers and Polynomials Lesson 4: Plotting and Printing Lesson 5: Illustrative Example -- Parallel RLC Circuit (Using MATLAB Help features and defining one's own Scripts/Functions using M-files) Select a lesson [1,2,3,4,5]: 3 %%%%%%% Lesson 3: Working with Complex Numbers and Polynomials % % 1) Complex Numbers: operations and special related functions % 2) Polynomials: MATLAB representation and special related functions % % At every keypress, one or more MATLAB commands is executed and an % explanation will follow the format below: (Actual MATLAB syntax) %(brief explanation of the command) (MATLAB's response to the issued command) pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%(ready for the next keypress) %%%%%%%%%%%%%%%% Complex numbers %MATLAB interprets the letter i to represent sqrt(-1) sqrt(-1) %PROVIDED THE USER HASN"T DEFINED IT AS SOMETHING ELSE ans = 0 + 1.0000i s1 = 3 + 4i %Defining complex numbers s1 = 3.0000 + 4.0000i s2 = 5 - 12j %Can use j also -- MATLAB caters to electrical engineers 8-) s2 = 5.0000 -12.0000i pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %MATLAB does "everything" with complex numbers as with real!!! A = [s1 j; j s2] %Generate matrix of complex entries A = 3.0000 + 4.0000i 0 + 1.0000i 0 + 1.0000i 5.0000 -12.0000i ca = s1 + s2 %Addition ca = 8.0000 - 8.0000i cm = s1*s2 %Multiplication cm = 63.0000 -16.0000i A2 = A.^2 %Element-wise Exponent -- each a2 = a*a A2 = 1.0e+02 * -0.0700 + 0.2400i -0.0100 -0.0100 -1.1900 - 1.2000i Ainv = inv(A) %Matrix inverse Ainv = 0.1176 - 0.1581i 0.0037 - 0.0147i 0.0037 - 0.0147i 0.0294 + 0.0699i pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %MATLAB also provides special functions for complex numbers a = real([s1 s2]) %extract real part -- element-wise a = 3 5 b = imag([s1 s2]) %extract imaginary part -- element-wise b = 4 -12 r = abs([s1 s2]) %calculate magnitude -- element-wise r = 5 13 alpha = angle([s1; s2]) %calculate phase angle in rad -- element-wise alpha = 0.9273 -1.1760 beta = degrees(angle([s1; s2])) %nest functions if want in degrees!! beta = 53.1301 -67.3801 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !clear %%%%%%%%%%%%%%% Polynomials %MATLAB represents polynomials by a vector of its coefficients %in order of descending powers p1 = [1 4 4] %Represents the polynomial s^2 + 4s + 4 p1 = 1 4 4 p2 = [2 0 0 -1] %Represents the polynomial 2s^3 +0s^2 + 0s - 1 p2 = 2 0 0 -1 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %MATLAB also provides special functions for polynomials r1 = roots(p1) %two roots of second-order polynomial -- via "roots" r1 = -2 -2 r2 = roots(p2) %three roots of third-order polynomial -- via "roots" r2 = -0.3969 + 0.6874i -0.3969 - 0.6874i 0.7937 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% p3 = poly(r2) %Coefficients of polynomial with given set of roots p3 = 1.0000 0.0000 0 -0.5000 a = polyval(p3,5) %Evaluate polynomial p3 at s = 5 a = 124.5000 c = polyval(p1,r1(1,1)) %Evaluating at root gives zero...whew! c = 0 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Coefficients of product of two polynomials p4 = conv([1 1],[1 1]) %(s+1)*(s+1) = s^2 + 2s + 1 p4 = 1 2 1 poly2str(p4,'s') %display polynomial as equation in variable s ans = s^2 + 2 s + 1 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% poly2str(conv([1 1],poly([j; -j])),'z') %Can combine operations in many ways ans = z^3 + z^2 + z + 1 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %We use polynomials primarily to express system functions num = [1 0]; den = [1 1 5]; printsys(num,den,'s') %Show system function in "pretty" form num/den = s ------------ s^2 + s + 5 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% echo off >> picocrse %%%%%%% Introduction MATLAB is a specialized software system for numerical computations with matrices. We will see that very sophisticated functions and computation routines can be built from the basic matrix operations within MATLAB. This 2-hour "Pico"-Course consists of a series of brief lessons intended to introduce the syntax of MATLAB and especially to present some commonly used functions and techniques. In addition, some very specific topics useful for 6.003 are discussed in more detail. We *strongly advise* that the student becomes accustomed to plotting capabilities, Help features, and defining one's own specialized functions using M-files in MATLAB. With this goal in mind, we close our presentation by using the power of MATLAB to analyze the step response for the parallel RLC circuit. The overall presentation consists of five lessons: Lesson 1: Generating Matrices Lesson 2: Basic Operations: Matrix and Element-wise Lesson 3: Working with Complex Numbers and Polynomials Lesson 4: Plotting and Printing Lesson 5: Illustrative Example -- Parallel RLC Circuit (Using MATLAB Help features and defining one's own Scripts/Functions using M-files) Select a lesson [1,2,3,4,5]: 4 %%%%%%% Lesson 4: Plotting and Printing % % 1) Generating and annotating simple plots of signals % 2) Multiple plots on a page, multiple curves on a set of axes, and advanced % annotation of plots % 3) Other useful types of plots % % At every keypress, one or more MATLAB commands is executed and an % explanation will follow the format below: (Actual MATLAB syntax) %(brief explanation of the command) (MATLAB's response to the issued command) pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%(ready for the next keypress) t1 = linspace(0,2,11); %Generate signal to be plotted y1 = sin(2*pi*t1); %Sinusoid with period 1 sec plot(t1,y1) %Open figure window and generate plot %NOTICE: plot does not appear sinusoidal...? pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% t2 = linspace(0,2,51); %Regenerate signal using more data points y2 = sin(2*pi*t2 + pi/2); %Sinusoid with period 1 sec shifted by 90 deg figure(1),plot(t2,y2) %Make figure window #1 active (rather than open %a second figure window) and plot pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Annotate plot in current figure window title('Sinusoid y2(t) shifted by 90 deg') %Add title xlabel('t (sec)') %Label x-axis ylabel('Amplitude') %label y-axis pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(1) %Make figure window #1 active %Can put many plots in same figure window subplot(2,1,1) %Generate 2x1 figure window, activate plot 1 plot(t1,y1) %Generate plot 1 title('Sinusoid y1(t) with too few datapoints') %Add title xlabel('t (sec)') %Label x-axis ylabel('Amplitude') %label y-axis subplot(2,1,2) %Activate plot 2 plot(t2,y2) %Generate plot 2 title('Sinusoid y2(t) shifted by 90 deg') %Add title xlabel('t (sec)') %Label x-axis ylabel('Amplitude') %label y-axis pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(1) %Make figure window #1 active plot(t1,y1,'y-',t2,y2,'r--') %Generate plots using different line styles title('Sinusoid y1(t) and y2(t)') %Add title xlabel('t (sec)') %Label x-axis ylabel('Amplitude') %label y-axis pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Can customize figure further text(0.4,0.8,'y1(t)') %Label curve for y1(t) text(0.15,-0.5,'y2(t)') %Label curve for y2(t) axis([0,2,-1.2,1.2]) %Can change range of figure axes grid on %Can add a grid... pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% grid off %...and remove it if clutters the plot pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% y3 = exp(-4*t2); %Generate a third signal--decaying exponential %Plot new curve without erasing previous plots hold on, plot(t2,y3,'g:'), hold off title('Sinusoid y1(t) and y2(t), and Exponential y3(t)') text(0.6,0.3,'y3(t)') pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %To send active figure to printer, simply type % print print -deps plot1.ps %Prints active figure to postscript file, can %then check it before printing on paper 8-) keyboard %To continue, type: return K>> K>> t1 = linspace(0,2,11); %Generate signal to be plotted t1 = linspace(0,2,11); %Generate signal to be plotted K>> y1 = sin(2*pi*t1); %Sinusoid with period 1 sec y1 = sin(2*pi*t1); %Sinusoid with period 1 sec K>> plot(t1,y1) plot(t1,y1) K>> clf clf K>> plot(t1,y1) plot(t1,y1) K>> print plot0.ps print plot0.ps K>> t2 = linspace(0,2,51); %Regenerate signal using more data points t2 = linspace(0,2,51); %Regenerate signal using more data points K>> y2 = sin(2*pi*t2 + pi/2); %Sinusoid with period 1 sec shifted by 90 deg y2 = sin(2*pi*t2 + pi/2); %Sinusoid with period 1 sec shifted by 90 deg K>> figure(1),plot(t2,y2) %Make figure window #1 active (rather than open figure(1),plot(t2,y2) %Make figure window #1 active (rather than open K>> %Annotate plot in current figure window %Annotate plot in current figure window K>> title('Sinusoid y2(t) shifted by 90 deg') %Add title title('Sinusoid y2(t) shifted by 90 deg') %Add title K>> xlabel('t (sec)') %Label x-axis xlabel('t (sec)') %Label x-axis K>> ylabel('Amplitude') ylabel('Amplitude') K>> print plot01.ps print plot01.ps K>> %Make figure window #1 active %Make figure window #1 active K>> %Can put many plots in same figure window %Can put many plots in same figure window K>> subplot(2,1,1) %Generate 2x1 figure window, activate plot 1 subplot(2,1,1) %Generate 2x1 figure window, activate plot 1 K>> plot(t1,y1) %Generate plot 1 plot(t1,y1) %Generate plot 1 K>> title('Sinusoid y1(t) with too few datapoints') %Add title title('Sinusoid y1(t) with too few datapoints') %Add title K>> xlabel('t (sec)') %Label x-axis xlabel('t (sec)') %Label x-axis K>> ylabel('Amplitude') %label y-axis ylabel('Amplitude') %label y-axis K>> subplot(2,1,2) %Activate plot 2 subplot(2,1,2) %Activate plot 2 K>> plot(t2,y2) %Generate plot 2 plot(t2,y2) %Generate plot 2 K>> title('Sinusoid y2(t) shifted by 90 deg') %Add title title('Sinusoid y2(t) shifted by 90 deg') %Add title K>> xlabel('t (sec)') %Label x-axis xlabel('t (sec)') %Label x-axis K>> ylabel('Amplitude') ylabel('Amplitude') K>> print plot02.ps print plot02.ps K>> return return pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close(1) %Close figure window #1 pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Several other useful types of plots subplot(2,2,1),stem(t2,y3) %Plots discrete "stems" title('Discrete Stem Plot: y = exp(-4t)') subplot(2,2,2),loglog(t2,y3) %Log scale on both axes title('Log Scale on Both Axes: y = exp(-4t)') grid on subplot(2,2,3),semilogx(t2,y3) %Log scale on x-axis only title('Log Scale on x-axis Only: y = exp(-4t)') grid on subplot(2,2,4),semilogy(t2,y3) %Log scale on y-axis only title('Log Scale on y-axis Only: y = exp(-4t)') grid on %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% NOTE: MAXIMIZE THE FIGURE WINDOW TO IMPROVE VIEW %%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% echo off >> %Several other useful types of plots >> subplot(2,2,1),stem(t2,y3) %Plots discrete "stems" ??? Undefined function or variable 't2'. >> title('Discrete Stem Plot: y = exp(-4t)') >> subplot(2,2,2),loglog(t2,y3) %Log scale on both axes ??? Undefined function or variable 't2'. >> title('Log Scale on Both Axes: y = exp(-4t)') >> grid on >> subplot(2,2,3),semilogx(t2,y3) %Log scale on x-axis only ??? Undefined function or variable 't2'. >> title('Log Scale on x-axis Only: y = exp(-4t)') >> grid on >> subplot(2,2,4),semilogy(t2,y3) %Log scale on y-axis only ??? Undefined function or variable 't2'. >> title('Log Scale on y-axis Only: y = exp(-4t)') >> grid on >> t2 = linspace(0,2,51); %Regenerate signal using more data points >> y2 = sin(2*pi*t2 + pi/2); %Sinusoid with period 1 sec shifted by 90 deg >> %Several other useful types of plots >> subplot(2,2,1),stem(t2,y3) %Plots discrete "stems" ??? Undefined function or variable 'y3'. >> title('Discrete Stem Plot: y = exp(-4t)') >> subplot(2,2,2),loglog(t2,y3) %Log scale on both axes ??? Undefined function or variable 'y3'. >> title('Log Scale on Both Axes: y = exp(-4t)') >> grid on >> subplot(2,2,3),semilogx(t2,y3) %Log scale on x-axis only ??? Undefined function or variable 'y3'. >> title('Log Scale on x-axis Only: y = exp(-4t)') >> grid on >> subplot(2,2,4),semilogy(t2,y3) %Log scale on y-axis only ??? Undefined function or variable 'y3'. >> title('Log Scale on y-axis Only: y = exp(-4t)') >> grid on >> picocrse %%%%%%% Introduction MATLAB is a specialized software system for numerical computations with matrices. We will see that very sophisticated functions and computation routines can be built from the basic matrix operations within MATLAB. This 2-hour "Pico"-Course consists of a series of brief lessons intended to introduce the syntax of MATLAB and especially to present some commonly used functions and techniques. In addition, some very specific topics useful for 6.003 are discussed in more detail. We *strongly advise* that the student becomes accustomed to plotting capabilities, Help features, and defining one's own specialized functions using M-files in MATLAB. With this goal in mind, we close our presentation by using the power of MATLAB to analyze the step response for the parallel RLC circuit. The overall presentation consists of five lessons: Lesson 1: Generating Matrices Lesson 2: Basic Operations: Matrix and Element-wise Lesson 3: Working with Complex Numbers and Polynomials Lesson 4: Plotting and Printing Lesson 5: Illustrative Example -- Parallel RLC Circuit (Using MATLAB Help features and defining one's own Scripts/Functions using M-files) Select a lesson [1,2,3,4,5]: 5 %%%%%%% Lesson 5: Illustrative Example -- Parallel RLC Circuit % % 1) Problem Statement (mimic types of statements you will get in the 6.003 HW) % 2) Basic MATLAB Solution: M-files are rlcstep.m (function) and % rlc1.m (script) % % (Consult the Lesson 5 supplementary handout which describes the analytical % solutions that this solution relies upon) % % 3) Advanced MATLAB solution: M-file is rlc2.m (script) % % (Uses advanced and useful MATLAB commands so solution does not rely on % analytical solution!) %%%%%%% Parallel RLC Circuit: Problem Statement (mimic HW problems) Consider a parallel RLC circuit driven by a current source. Determine the step response for the following sets of component values, and use MATLAB to plot each response from t = [0,T]: Case (i) : R = 1 Ohm, C = 1 Farad, L = 6 Henrys, T = 20 sec Case (ii) : R = 1 Ohm, C = 1 Farad, L = 1 Henry , T = 20 sec Case (iii): R = 1 Ohm, C = 1 Farad, L = 4 Henrys, T = 20 sec Basic Note: You may find the following commands useful: plot, title, xlabel, ylabel Advanced Note: You may find the following commands useful: residue, step % Outline of Solution: We will first demonstrate how to create a % specialized function and script (using M-files) % involving the already familiar commands given in the % Basic Note. This solution method requires that we % first obtain the analytical solution of the step % response in our circuit. % % Then, we will make use of MATLAB's Help features to % learn about the currently unfamiliar, yet very powerful % commands given in the Advanced Note. We will % create another script (M-file) that uses these advanced % commands, and see that we will have saved ourselves % a significant amount of work because we will not % require grinding through the analytical solution! pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !clear %%%%% Basic Solution % % The approach we take here is to define a MATLAB function which takes as % inputs values for R, L, C, and T and gives the corresponding step response % as the output. This function we will call rlcstep.m The equations in % this file come from the analytical solution (consult the Lesson 5 % supplementary handout). % Then, we write a MATLAB script rlc1.m which calls on function % rlcstep.m for each of the three cases, and includes the necessary plot % commands. % % NOTE: While studying this solution, we will also be introduced to MATLAB's % C language-based "control flow" operations (e.g., if-elseif-else % conditional statements, while loops, for loops, etc.). pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !clear type rlcstep.m %View the MATLAB function M-file function V = rlcstep(R,L,C,T) % function V = rlcstep(R,L,C,T) % % % RLCSTEP Function RLCSTEP returns the step response of a parallel RLC % circuit for the given input parameters. Input variables R, L, % and C denote the values of the circuit components: resistor (in % Ohms), inductor (in Henrys), and capacitor (in Farads), % respectively. Output variable V is a row vector representing % the step response v(t) from t = [0,T]. % % This function first checks whether the system poles (natural % frequencies) are real & distinct, complex conjugates, or real % & equal. Depending on which of these cases is true for the % given component values, the correct "analytical" solution is % applied (consult the Lesson 5 supplementary handout). At the % end of this example, we will solve the same problem much more % quickly, without relying on our analytical solution, by using % some of MATLAB's more sophisticated functions. t = linspace(0,T,201); %define the time vector dis = (R*C)^(-2) - 4/(L*C); %calculate the "discriminant" for the poles p1 = -1/(2*R*C) + sqrt(dis)/2; %compute the values for the two poles p2 = -1/(2*R*C) - sqrt(dis)/2; if dis > 0 %Case (i): poles are real, distinct K1 = 1/((p1-p2)*C); V = K1*(exp(p1*t) - exp(p2*t)); elseif dis < 0 %Case (ii): poles are complex conjugates K1 = 1/((p1-p2)*C); V = -2*imag(K1)*exp(real(p1)*t) .* sin(imag(p1)*t); else %Case (iii): poles are real, equal V = (1/C)*t .* exp(p1*t); end pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %NOTICE: 1) Function header line (distinguishes function from script) % 2) Help section (the many rows of commented lines at top of file) % 3) Use of an if-elseif-else statement % 4) Included in-line comments for future reference pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !clear %The comments section after the function header is the help rlcstep %help information! function V = rlcstep(R,L,C,T) RLCSTEP Function RLCSTEP returns the step response of a parallel RLC circuit for the given input parameters. Input variables R, L, and C denote the values of the circuit components: resistor (in Ohms), inductor (in Henrys), and capacitor (in Farads), respectively. Output variable V is a row vector representing the step response v(t) from t = [0,T]. This function first checks whether the system poles (natural frequencies) are real & distinct, complex conjugates, or real & equal. Depending on which of these cases is true for the given component values, the correct "analytical" solution is applied (consult the Lesson 5 supplementary handout). At the end of this example, we will solve the same problem much more quickly, without relying on our analytical solution, by using some of MATLAB's more sophisticated functions. pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !clear rlc1 %Run MATLAB script file % RLC1 MATLAB script file which generates the basic solution to % Lesson 5 of the 6.003 MATLAB "Pico"-Course. This script calls % on function rlcstep.m for each of the three cases, and % includes the necessary plot commands. % Case (i) : R = 1 Ohm, C = 1 Farad, L = 6 Henrys, T = 20 sec % Case (ii) : R = 1 Ohm, C = 1 Farad, L = 1 Henry , T = 20 sec % Case (iii): R = 1 Ohm, C = 1 Farad, L = 4 Henrys, T = 20 sec % Define component values T = 20; R = 1; C = 1; L = [6; 1; 4]; words = ['poles real,distinct'; ... % Cycle through the cases, generating all three plots on the same figure for case = 1:3 v(case,:) = rlcstep(R,L(case,1),C,T); %function call subplot(3,1,case),plot( linspace(0,T,size(v,2)), v(case,:) ); xlabel('t (sec)'); ylabel('v(t)'); text(10,0.4,words(case,:)); end v(case,:) = rlcstep(R,L(case,1),C,T); %function call subplot(3,1,case),plot( linspace(0,T,size(v,2)), v(case,:) ); xlabel('t (sec)'); ylabel('v(t)'); text(10,0.4,words(case,:)); end v(case,:) = rlcstep(R,L(case,1),C,T); %function call subplot(3,1,case),plot( linspace(0,T,size(v,2)), v(case,:) ); xlabel('t (sec)'); ylabel('v(t)'); text(10,0.4,words(case,:)); end subplot(3,1,1) title('Parallel RLC Circuit: Step Response for Three Cases'); pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear !clear %%%%% Advanced Solution % % The approach we take here is to rely on the powerful functions given in % the advanced note of the problem statement. The first step is to use % MATLAB's Help features to learn exactly how these functions work. Then, % we write a MATLAB script rlc2.m which calls on these functions for % each of the three cases. These functions will allow us to determine % the form of the step response, in addition to plotting them, without % relying on the analytical solutions! pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !clear help residue %Learn about function residue RESIDUE Partial-fraction expansion (residues). [R,P,K] = RESIDUE(B,A) finds the residues, poles and direct term of a partial fraction expansion of the ratio of two polynomials B(s)/A(s). If there are no multiple roots, B(s) R(1) R(2) R(n) ---- = -------- + -------- + ... + -------- + K(s) A(s) s - P(1) s - P(2) s - P(n) Vectors B and A specify the coefficients of the numerator and denominator polynomials in descending powers of s. The residues are returned in the column vector R, the pole locations in column vector P, and the direct terms in row vector K. The number of poles is n = length(A)-1 = length(R) = length(P). The direct term coefficient vector is empty if length(B) < length(A), otherwise length(K) = length(B)-length(A)+1. If P(j) = ... = P(j+m-1) is a pole of multplicity m, then the expansion includes terms of the form R(j) R(j+1) R(j+m-1) -------- + ------------ + ... + ------------ s - P(j) (s - P(j))^2 (s - P(j))^m [B,A] = RESIDUE(R,P,K), with 3 input arguments and 2 output arguments, converts the partial fraction expansion back to the polynomials with coefficients in B and A. Warning: Numerically, the partial fraction expansion of a ratio of polynomials represents an ill-posed problem. If the denominator polynomial, A(s), is near a polynomial with multiple roots, then small changes in the data, including roundoff errors, can make arbitrarily large changes in the resulting poles and residues. Problem formulations making use of state-space or zero-pole representations are preferable. See also POLY, ROOTS, DECONV. pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !clear help step %Learn about function step STEP Step response of LTI systems. STEP(SYS) plots the step response of each input channel of the LTI system SYS (created with either TF, ZPK, or SS). The time range and number of points are chosen automatically. STEP(SYS,TFINAL) simulates the step response from t = 0 to the final time t = TFINAL. For discrete-time systems with unspecified sampling time, TFINAL is interpreted as the number of samples. STEP(SYS,T) uses the user-supplied time vector T for simulation. For discrete-time systems, T should be of the form Ti:Ts:Tf where Ts is the sample time of the system. For continuous systems, T should be of the form Ti:dt:Tf where dt will become the sample time of a discrete approximation to the continuous system. The step input is always assumed to start at t = 0 (regardless of Ti). STEP(SYS1,SYS2,...,T) plots the step response of multiple LTI systems SYS1,SYS2,... on a single plot. The time vector T is optional. You can also specify a color, line style, and marker for each system, as in step(sys1,'r',sys2,'y--',sys3,'gx'). When invoked with left-hand arguments, [Y,T] = STEP(SYS,...) returns the output response Y and the time vector T used for simulation. No plot is drawn on the screen. If SYS has NU inputs and NY outputs, and LT=length(T), the array Y is LT-by-NY-by-NU and Y(:,:,j) gives the step response of the j-th input channel. For state-space systems, [Y,T,X] = STEP(SYS,...) also returns the state trajectory X which is an LT-by-NX-by-NU array if SYS has NX states. See also INITIAL, IMPULSE, LSIM. Overloaded methods help lti/step.m pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !clear rlc2 %Run MATLAB script file % RLC2 MATLAB script file which generates the advanced solution to % Lesson 5 of the 6.003 MATLAB "Pico"-Course. This script calls % on functions residue.m and step.m for each of the % three cases, and includes the necessary plot commands. % Case (i) : R = 1 Ohm, C = 1 Farad, L = 6 Henrys, T = 20 sec % Case (ii) : R = 1 Ohm, C = 1 Farad, L = 1 Henry , T = 20 sec % Case (iii): R = 1 Ohm, C = 1 Farad, L = 4 Henrys, T = 20 sec % Define component values T = 20; R = 1; C = 1; L = [6; 1; 4]; words = ['poles real,distinct'; ... t = linspace(0,T,201); % Cycle through the cases, generating all three plots on the same figure for case = 1:3 Hnum = [1/C 0]; %Define system function Hden = [1 1/(R*C) 1/(L(case,1)*C)]; %To continue lesson, type: return keyboard %Use printsys to verify system function K>> K>> print plotles5.ps print plotles5.ps K>> return return [K,P] = residue(Hnum,Hden) %Determine Partial Fraction Expansion K = 1.3660 -0.3660 P = -0.7887 -0.2113 v(case,:) = step(Hnum,Hden,t)'; %function call subplot(3,1,case), hold on plot( t, v(case,:), 'rx'); xlabel('t (sec)'); ylabel('v(t)'); text(10,0.4,words(case,:)); end Hnum = [1/C 0]; %Define system function Hden = [1 1/(R*C) 1/(L(case,1)*C)]; %To continue lesson, type: return keyboard %Use printsys to verify system function K>> print plotles51.ps print plotles51.ps K>> return return [K,P] = residue(Hnum,Hden) %Determine Partial Fraction Expansion K = 0.5000 + 0.2887i 0.5000 - 0.2887i P = -0.5000 + 0.8660i -0.5000 - 0.8660i v(case,:) = step(Hnum,Hden,t)'; %function call subplot(3,1,case), hold on plot( t, v(case,:), 'rx'); xlabel('t (sec)'); ylabel('v(t)'); text(10,0.4,words(case,:)); end Hnum = [1/C 0]; %Define system function Hden = [1 1/(R*C) 1/(L(case,1)*C)]; %To continue lesson, type: return keyboard %Use printsys to verify system function K>> return return [K,P] = residue(Hnum,Hden) %Determine Partial Fraction Expansion K = 1.0000 -0.5000 P = -0.5000 -0.5000 v(case,:) = step(Hnum,Hden,t)'; %function call subplot(3,1,case), hold on plot( t, v(case,:), 'rx'); xlabel('t (sec)'); ylabel('v(t)'); text(10,0.4,words(case,:)); end subplot(3,1,1) title('Parallel RLC Circuit: Step Response for Three Cases'); pause %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% echo off >> diary off