Join StudySoup for FREE

Get Full Access to
UTA - MAE 3360 - Class Notes - Week 14

Description

Reviews

Application of Math Principles To Engineering Problems 8th edition By Professor David A. Hullender January 1, 20171 Preface The objectives of this notebook are to introduce solution techniques for solving typical equations encountered in the analysis and simulation of engineering systems. When obtaining solutions to equations, it is very beneficial to focus on the common sense details associated with the performance of the actual engineering system for which the equations are written. Consequently, throughout this text, equations for actual engineering systems are utilized so that it is possible to apply common sense to the prediction and approximation of the solution to the equations before actually solving the equations. Knowing the approximate solution provides confirmation to the actual solution in the end. From the very beginning of the notebook, computer algorithms in MATLAB will be utilized to obtain and plot solutions for typical equations encountered when solving engineering problems. Having experience with MATLAB prior to reading this notebook is not necessary. Specific details and examples of using MATLAB for solving engineering equations and for simulating dynamic systems are included in order for the reader not experienced with MATLAB to be able to utilize the software effectively without having to be creative and start from scratch. At the end of several of the sections, a list of previous quiz, homework, and exam problems are listed that utilize some of the concepts introduced in that particular section. The solutions to almost all of the problems are provided; however, it is very important to test your understanding of the concepts by trying to work each problem without looking at the solution. Then and only then will you know your weaknesses; it is better to determine your weaknesses prior to an exam instead of during an exam! Reviewing the solution to a problem before trying to work it without looking at the solution is a total waste of time; in my opinion, this is the number one mistake students make when preparing for exams. The second greatest mistake is preparing for an exam with other students; you are not learning what you don’t know and, unfortunately, when you take an exam, that is when you will find out what you don’t know. Table of Contents Chapter 1 Basic Math Principles 1.1 Complex Numbers 5 1.2 Notation for Time Derivatives and Use of the ‘D’ Operator 6 1.2.1 Notation for Time Derivatives 6 1.2.2 Use of the ‘D’ Operator to convert Differential Eqn’s to Algebraic Eqn’s 6 1.3 Eigenvalues of a System 7 1.3.1 General Properties of Eigenvalues 7 1.3.2 Using MATLAB to Compute Eigenvalues 8 1.4 Solving Simultaneous Algebraic Equations 9 1.4.1 Manually Solving Simultaneous Equations 9 1.4.2 Using Symbolic Math in MATLAB to Solve Simultaneous Equations 112 Chapter 2. Solving Differential Equations 2.1 Analytical Solutions of Differential Equations Using the Laplace Transform 15 2.1.1 Linearizing Nonlinear Differential Equations 15 2.1.2 The definition of the Laplace Transform 20 2.1.3 Properties of the Laplace Transform 20 2.1.4 Final Value Theorem 21 2.1.5 Initial value theorem 21 2.1.6 Converting Differential Eqn’s to Algebraic Eqn’s Using Laplace Transform 22 2.1.7 Common Inputs and Their Laplace transforms 23 2.1.7.1 Step Functions and Constants 23 2.1.7.2 Pulse and Impulse Functions 24 2.1.7.3 Exponential Functions 25 2.1.7.4 Periodic Functions 25 2.1.8 Inverse Laplace Transform 25 2.1.8.1 Inverse Laplace Transform Using Partial Fractions 25 Obtaining Partial Fractions Using MATLAB 28 2.1.8.2 Inverse Laplace Transform Using the Residue Theorem 28 Shortcut for Complex Poles 29 2.1.8.3 Using MATLAB to Plot the Inverse Laplace Transform without Getting an Equation for the Inverse Laplace Transform 32 2.1.9 Modes of a System 34 2.2 Numerical Solutions to Differential Equations Using MATLAB 36 2.2.1 General Considerations 36 2.2.2 Using the ‘impulse’ command in MATLAB 36 2.2.3 Using the ‘step’ command in MATLAB 37 Obtaining DC Gain using ‘dcgain’ command 37 2.2.4 Using the Command ‘lsim’ 38 Pulse Series input 39 2.2.5 Using the Command ‘initial’ 41 2.3 Using MATLAB to Get the System Response Using Numerical Integration 41 2.3.1 Example Problems using ‘ode45’ in MATLAB 42 Deterministic Inputs 42 Random Inputs 43 2.3.2 Example Simulation Using Options and Event Functions with ‘ode45’ 44 2.4 Using MATLAB to Get the System Steady State Response for Periodic Inputs 46 2.4.1 Frequency Response Analysis 46 2.4.2 Examples of Frequency Response Function Plots 47 2.5 Using MATLAB to Get Lower Order Approximations for Transfer Functions 51 2.5.1 Examples of Modal Approximation 51 2.5.2 Transfer Function Approximations using Inverse Freq. Domain Analysis 54 2.6 Converting Linear Continuous Differential Eqn’s to Discrete Time Eqn’s 65 2.6.1 Converting 1st Order Differential Eqn’s Using Inverse Laplace transform 65 2.6.2 Conversion Using the Continuous-to-discrete MATLAB command, ‘c2d’ 663 Chapter 3. Expressing Differential Equations in State Variable Format 3.1 Representing an nth Order System by n First Order Differential Eqn’s. 68 3.2 Simulation Diagram Approach When There are Input Derivatives 69 3.3 General Matrix Format for State Variable Equations 70 3.4 State Variable Equations in Matrix Format if there are Input Derivatives 71 3.4.1 Only 1st Order Input Derivatives 71 3.4.2 Format Options for Input Derivatives Higher than 1st Order 72 3.4.2.1 Simulation Diagram Approach 72 3.4.2.2 Phase Variable Approach 73 3.4.2.3 MATLAB Approach 73 3.4.2.4 Example Demonstrating the Different Methods 74 3.4.2.5 Example of Phase Variable Approach when the Numerator and Denominator Orders are Equal 75 3.5 Maintaining Access to Natural Variables When Creating State Variable Equations from Block Diagrams 76 3.6 SIMULINK 77 Chapter 4 Example Problems 81 Chapter 5 Solutions to Example Problems in Chapter 4 103 Chapter 6 Previous Quiz Problems and Solutions 134 Chapter 7 Previous Homework and Solutions 151 Chapter 8 Previous Exams and Solutions 261 Chapter 9 MATLAB M-files 3234 Chapter 1 Basic Math Principles 1.1 Complex numbers A complex number can be expressed in exponential format, i.e. y a jb Me j = + = 2 2 M a b = + φ φ = + angle a jb ( ) Note, the angle is in one of four quadrants. Simple use of the inverse tangent on a calculator may give the incorrect answer. In MATLAB, the command ‘angle’ with a complex number gives the correct answer for any of the four quadrants. Consider the examples below. imaginary part b M Example: y j ej = − − = 4 3 5 φ 3.785 real part a M 2 2 = + = 3 4 5 − 1 φ angle j π = − − = + = ( 4 3) tan (3/ 4) 3.785 imaginary part 3.785 rad -4 -3 real part Example: The solution y(t) to a differential equation is found to be the following:5 ��(��) = 3 + 25����−4��sin(3�� + ∅) where �� = |−2 + ��3 −4 + ��3| ������ ∅ = ���������� (−2 + ��3 −4 + ��3) What are M and ∅? �� =|−2 + ��3| |−4 + ��3|=√(−2)2 + 32 √(−4)2 + 32=√13 √25= 0.7211 ∅ = ����������(−2 + ��3) − ����������(−4 + ��3) = 2.159 − 2.498 = −0.339 ������ 1.2 Notation for Time Derivatives and Use of the ‘D’ Operator 1.2.1 Notation for Time Derivatives Throughout this book, time derivatives of functions will be denoted with ‘dots’ over the function. For example, for the first three time derivatives of function y(t) ��̇ ≡�������� ��̈≡��2�� ����2��⃛ ≡��3�� ����3 1.2.2 Use of the ‘D’ Operator to convert linear Differential Equations to Algebraic Equations The ‘D’ operator is used to express a linear differential equation with an input and with no initial conditions into an algebraic equation. In ‘D’ operator format, the first derivative with respect to time of a function, f(t), is denoted by Df(D) where f(D) denotes the function f(t) in ‘D’ operator format. The 2nd derivative of f(t) in D operator format is denoted by D2f(D), etc. Consider the following 3rd order linear ordinary differential equation for the output y(t) of a system with input u(t): ����⃛+16��̈+ ������̇ + �������� = ����̈+ +������̇ + �������� Initial conditions ��(0−), ��̇(0−), ��̈(0−), ��(0−), ������ ��̇(0−) are in general needed as well as the input function u(t) to obtain the solution for y(t). If all of these initial conditions are zero, then the ‘D’ operator can be used to express the differential equation as an algebraic equation. Note, if any of the initial conditions are not zero, then the Laplace operator ‘s’ will be required instead of ‘D’. In “D’ operator format, the differential equation above is written as follows: (������ + �������� + ������ + ������)��(��) = (������ + ������ + ������)��(��) Solving for y(D) gives ��(��) = [������ + ������ + ������ ������ + �������� + ������ + ������]��(��) The function of ‘D’ in the brackets is called the transfer function relating the input u to the output y. Note that the numerator and denominator polynomial coefficients in the transfer function are simply the coefficients in the differential equation. The denominator of the transfer function is called the characteristic polynomial; the poles of the transfer function (roots of the characteristic polynomial) are called the eigenvalues of the system. Knowing the transfer function and eigenvalues of a system are very beneficial to analyzing time response dynamics for different system inputs.6 In general if u(t) is a constant, then the solution of the differential equation for y(t) will also become a constant as time goes to infinity; when y(t) becomes a constant, ��̇(��) will become zero as will ��̈(��), ��⃛(��), ��̇(��), and ��̈(��). Thus, in the end (as time goes to infinity), the original differential equation above simplifies to 0 + 0 + 0 + 100�� = 0 + 0 + 150�� This equation can be solved for the limiting value (also called the final value) of y(t), i.e. �� = 150 100�� where u is the constant value for the input. For example, if u =4, then the final value of y is 6. 1.3 Eigenvalues of a System 1.3.1 General Properties of Eigenvalues The eigenvalues of a system define the time response characteristics of a system resulting from an external input or from nonzero initial conditions. The eigenvalues are obtained by factoring the characteristic polynomial which is the denominator of the transfer function of the system. Eigenvalues will either be real or complex. The real parts of all of the eigenvalues of the system must be negative otherwise the system will be unstable; thus, a single eigenvalue, real or complex, with a positive real part will cause the time response of the system to diverge to infinity. A complex pair of eigenvalues, −�� ± ����, corresponds to the following polynomial formats for a polynomial: (�� + �� + ����)(�� + �� − ����) = [(�� + ��)2 + ��2] = ��2 + 2���� + ��2 + ��2 = ��2 + 2�������� + ����2 where -r is the real part of the eigenvalue, �� is the imaginary part also called the damped natural frequency, �� is the damping ratio, and ���� is the undamped natural frequency. Note that �� =������. Also note that eigenvalues have time constants, ��; �� = |1��|. Time constants determine how long it will take for the time response of the system to essentially end; it takes approximately 5 times the largest eigenvalue time constant for the response to be within 1% of its final value. So, when performing simulations for problems in which the solution becomes zero or a constant, a rough estimate of the final time should be five times the largest time constant, i.e. �������������������� ���������� �������� ≈ 5�������� Example: A system has the following eigenvalues: −5, −2, −4 ± ��3 The time constants of this system are 1/5=0.2 sec., 1/2=0.5 sec., 1/4=0.25 sec. The largest time constant is 0.5; thus, it will take about 5*0.5=2.5 seconds for this system to return to equilibrium following an input disturbance or nonzero initial condition. Usually, ‘returning to equilibrium’ means that the solution to the equations become constant. For the complex pair of eigenvalues −4 ± ��3, we can write (�� + 4 + ��3)(�� + 4 − ��3) = (�� + 4)2 + 32 = ��2 + 8�� + 25 Thus, ���� = 3������ ������, ���� = √25 = 5������ ������, ������ �� = 0.8 Eigenvalue problems in this notebook: 4.16, 4.13, 4.16, 4.22, 4.25, 4.26, 4.27, 4.28, 4.29, 4.30, 4.36, 4.37, 4.40, 4.43, 5.6, 5.13, 5.16, 5.22, 5.25, 5.26, 5.27, 5.28, 5.29, 5.30, 5.33, 5.36, 5.37, 5.43, Q7, Q14, Q18, H3, H5, H9, H12, H13, H20, H21, H24, H25, H28, H34, H35, H43, E1C, E2A, E2C, E3B, E3D, E4B,E5A, E5C, E6A, E6C7 1.3.2 Using MATLAB to Compute Eigenvalues All MATLAB commands are defined using lower case letters. The best way to learn to use any command is to type ‘help command-name’. Using ‘help’ works for all of the MATLAB commands. If you don’t know the name of a command, type ‘help’ to get a list of the different categories of commands. For example, typing ‘help’ lists ‘elfun’ for elementary functions. To see a list of these elementary functions type ‘help elfun’. As mentioned in section 1.2.2, the transfer function defines how the output of a system responds to the input to that system. Also, the denominator of a transfer function is the characteristic polynomial and the roots of this polynomial are called the eigenvalues of the system. If the numerator and denominator polynomials of a transfer function are known, then the transfer function can be entered into MATLAB. Consider the transfer found in section 1.2.2: ��(��) = [������+������+������ ������+��������+������+������] ��(��) The procedure for entering the coefficients of polynomials and creating transfer functions is shown below. Note, putting a semi-colon at the end of a statement stops MATLAB from echoing the result; this is really important if the result consists of hundreds of numbers. In the MATLAB commands below, the polynomials are given names ‘num’ and ‘den’; the transfer function is given the name ‘G’: >> num=[6 28 150]; Entering the numerator polynomial into MATLAB >> den=[2 16 74 100]; Entering the denominator polynomial into MATLAB >> G=tf(num,den) Entering the transfer function into MATLAB �� = [������ + ������ + ������ ������ + �������� + ������ + ������] The ‘damp’ command gives the eigenvalues, damping ratios, and undamped natural frequencies. >> damp(G) Eigenvalues Damping Time Undamped Natural Ratio Constant, sec. Frequency, rad/sec -2 not applicable 0.5 not applicable -3 ± j4 0.6 0.333 5 The time constants are computed using one over the absolute value of the real parts of the eigenvalues. For this example, the time constants are 1/2 and 1/3. Time constants are positive numbers. Note, MATLAB will incorrectly give values for the damping ratio and frequency for real eigenvalues. Real eigenvalues do not have a damping ratio and frequency. I have replaced the numbers with ‘not applicable’ in the MATLAB response to the ‘damp’ command. The eigenvalues can also be obtained from the denominator polynomial using the command ‘roots’, i.e. >> roots(den) 8 1.4 Solving Simultaneous Algebraic Equations The concept of solving simultaneous algebraic equations is included due to the fact that it will likely be the greatest source of solution challenges unless a simple appropriate approach is used and unless thorough attention to detail is given at each step. When solving engineering problems, we often start with a combination of differential equations and algebraic equations that represent a simultaneous set of equations defining the dynamics of the system to be simulated. When appropriate, the differential equations can first be converted to algebraic equations using either the ‘D’ operator or the Laplace transform so that all of the equations are then algebraic. 1.4.1 Manually Solving Simultaneous Equations There are numerous ways to solve simultaneous algebraic equations. For engineering problems, there is normally only a single variable representing the output of interest. Consequently, it is suggested that a simple step by step approach be used for simplifying and reducing the simultaneous equations down to a single equation with the variable representing the output of interest being the only remaining unknown. Consider the following three simultaneous algebraic equations with unknowns x, y, and z. (1) 2�� + 3�� − 10�� = 4 (2) − 5�� + 2�� + 8�� = 1 (3) 3�� − 4�� + 2�� = 0 These equations are considered simultaneous because you can’t solve for either x, y or z without using all three equations. Unless the solutions for all three unknowns are needed, an approach giving solutions for all of the unknowns, such as the matrix inverse approach, may be computationally burdensome and not necessary. Suppose x is the output of interest meaning we only want a solution for x. The suggested approach is to solve for one of the other unknowns in one of the equations and substitute that solution into the remaining equations; this process will reduce the number of equations by one and reduce the number of unknowns by one. This process is repeated until only a single equation remains. For instance, if we solve (1) for z, we get �� = 0.2�� + 0.3�� − 0.4 Substituting this result for z into (2) and (3) gives (2) − 5�� + 2�� + 8(0.2�� + 0.3�� − 0.4) = 1 (3) 3�� − 4�� + 2(0.2�� + 0.3�� − 0.4) = 0 or (2) − 3.4�� + 4.4�� = 4.2 (3) 3.4�� − 3.4�� = 0.8 We now have two equations and two unknowns x and y. Since we want the solution for x, we can solve (2) for y and substitute the result into (3). This gives �� = 5.23539 It is not uncommon for one or more of the equations to contain a symbolic parameter in place of a numerical value. This is always the case when solving engineering problems with the ‘D’ operator or with the Laplace transform. In such cases, the solution will contain this symbolic variable. For instance, suppose our original three equations contain the ‘D’ operator, i.e. (1) 2���� + 3�� − 10�� = 4 (2) − 5�� + 2��2�� + 8���� = 1 (3) 3�� − 4�� + 2�� = 0 If we use the same solution process while treating the symbol D as some numerical value, we get the following: �� =8��2 + 64�� + 17 4��3 + 62��2 + 36�� − 85 Although the step by step solution process is the same, close attention to detail and checking each step is required to avoid costly errors. Being able to solve simultaneous algebraic equations is often critical for analyzing and simulating the dynamics of an engineering system. It is important to note that before starting the solution process, one should always confirm that the number of unknowns matches the number of equations. If there are more unknowns than equations, the solution for the output of interest will be obtained in terms of the extra unknowns; in general, a shortage in the number of equations means the modeling process for the system is incomplete. Example: Consider the mass/spring/damper system shown below; the displacement of the mass is denoted by y(t) and the input force that moves the mass is denoted by f(t). y(t) Spring k Damper c Mass m Force Input f(t) An equation for each component can be written using modeling techniques. Note, the use of modeling techniques to obtain equations for engineering components and systems will be discussed more in depth in a later section of this notebook. The algebraic equation for the force ��1 on the spring in terms of the spring constant k=100 N/m is ��1 = 100�� The differential equation for the force ��2 on the damper in terms of the viscous damping coefficient c = 2 Ns/m is as follows: ��2 = 2��̇10 Using the ‘D’ operator, this differential equation can be converted to an algebraic equation, i.e. ��2(��) = 2����(��) And finally, the differential equation for the acceleration of the mass m = 4 kg is ��̈= (�� − ��1 − ��2)/4 Converting this differential equation to an algebraic equation gives ��2��(��) = [��(��) − ��1(��) − ��2(��)]/4 We started with one algebraic equation and two differential equations and used the ‘D’ operator to simplify the model to three algebraic equations with three unknowns: ��1(��), ��2(��), & ��(��). Don’t forget, f(D) is the input which will be a known function of time. Using the step by step reduction process to obtain a single equation, the equation for the output of interest y(D) is obtained, i.e. ��(��) = [1 Or 4��2 + 2�� + 100] ��(��) ��(��) = [1 ����2 + ���� + ��] ��(��) This single algebraic equation defines how the input force ��(��) affects the displacement y(t) of the mass. Just how to use this equation for analysis purposes is not important at this point. What is important, however, is that we were able to convert the differential equations to algebraic equations and then use simple algebraic reduction of the equations to get a single equation for the output of interest. For more complex systems consisting of many complex simultaneous equations, we can use Symbolic Math in MATLAB to obtain the single equation solution. Problems requiring this concept: 4.3, 4.20, 5.3, 5.21, Q3, Q11, H9, H23, H24, H25, H34, H43, E1a, E1c, E4a 1.4.2 Using Symbolic Math in MATLAB to Solve Simultaneous Equations Symbolic math allows a computer to obtain a solution in terms of symbols that may not have been assigned numerical values. Here we demonstrate how to use one of the symbolic commands ‘solve’ to get the solution for the three algebraic equations in section 1.4.1. The symbolic statement below provides the solution for all three unknowns x, y, and z. >> H=solve('2*x+3*y-10*z=4','-5*x+2*y+8*z=1','3*x-4*y+2*z=0','x,y,z') Note that each equation has been entered between single quotation marks in terms of the symbolic variables x, y, and z which are the three unknowns. Also, the three unknowns are listed between single quotation marks at the end of the statement. 11 To see the solution for x, type ‘H.x’, i.e. >> H.x Ans = 89/17 Since ‘solve’ is a symbolic math command, answers will be displayed as a ratio of two whole numbers. To see the answer in floating point arithmetic, you must first specify the number of digits to be displayed and then use the command ’vpa’, i.e. >> digits(5) >> vpa(H.x) Ans = 5.2353 Example: It is not uncommon for equations to contain parameters that are symbolic. The first step is to declare which parameters are to be symbolic. For the mass/spring/damper system in the section 1.5.1, note that the answer was obtained in terms of m, c, k, and D. Realizing this in advance alerts us to begin the solution process by declaring these to be symbolic. This is done with the following MATLAB command where ‘d’ has been used in place of ‘D’ to avoid a conflict since ‘D’ is used in the symbolic algorithms. The ‘syms’ command declares that the listed symbols may not have numerical values >> syms d m c k The next step is to use the MATLAB ‘solve’ command to get the solution. For the mass/spring/damper example problem, this is done as follows: >> H=solve(‘F1=k*y’,’F2=c*d*y’,’m*d^2*y=f-F1-F2’,’y,F1,F2’); Note as before, that each of the equations has been entered between single quotes. Following the equations, the three unknowns are listed also between single quotes. Executing these two commands gives the solution H which contains the solution for all three unknowns. The solution for only y has been stored in the computer memory as H.y; using the ‘pretty’ command prints the solution in a more recognizable format, i.e. >> pretty(H.y) which gives �� = [1 ����2 + ���� + ��] �� Since m, c, and k are symbolic, it is possible to assign numerical values to each. This is done with the ‘subs’ command. Note, if the numerical values are not simple integers, the ‘digits’ and ‘vpa’ commands may be needed to make the results presentable. For instance, if m = 4, c=2, and k=100:12 >> Y=subs(H.y,[m c k],[4 2 100]) Y=f/(4*d^2+2*d+100) >> pretty(Y) �� =1 4��2 + 2�� + 100 �� Additional commands that are useful when working with symbolic variables are (1) ‘collect’ for grouping like powers of a symbol such as ‘d’; (2) ‘numden’ for getting the numerator and denominator polynomials of a symbolic transfer function; (3) ‘poly2sym’ for converting a vector of polynomial coefficients to a symbolic polynomial; (4) ‘sym2poly’ for converting a symbolic polynomial to a vector of coefficients; and (5) ‘int8’ and ‘int16’ for converting and rounding off numbers with a decimal point to integers. For other useful symbolic commands, type ‘help symbolic’. Pile Driving Example Consider the following three lumped mass approximation model for a pile being driven into the ground with input force F(t); the output of interest is the force f(t) exerted on the ground resulting from the input force F(t). F(t) F(t) m=10 kg M z b=5 Ns/m k k k f M M f k=1000 N/m b v b w b It can be shown that the equations representing this model for displacements z, v, and w defined to be zero at equilibrium and assuming zero displacement at the bottom end of the pile are as 13 shown below. The ground force �� is the force that forces the bottom end of the pile into the ground. The output of interest to us is the ground force defined in the last equation. mz bz k z bv k v F + + = + + mv bv k v bz k z bw k w + + = + + + 2 2 mw bw k w bv k v + + = + 2 2 f mg bw k w = + + 3 The ‘3mg’ term represents a second input to the problem and will be neglected since it is relatively small compared to the impact force F(t) and since we are only wanting the transfer function relating f to the input F. Use the ‘D‘ operator to convert the differential equations to algebraic equations and then use symbolic math to get the transfer function relating the ground force f(t) to the input force F(t). Note, ‘d’ is used instead of ‘D’ to avoid conflict with symbolic notation for differentiation. % Symbolic math demo: Pile Driving Model clear all % clears memory format short e % numbers to be shown in scientific format syms m b k d F % variables and parameters to be assigned numerical values or to %be used in other symbolic commands G=solve('(m*d^2+b*d+k)*z=(b*d+k)*v+F','(m*d^2+2*b*d+2*k)*v=(b*d+k)*z+(b*d+k)*w',... '(m*d^2+2*b*d+2*k)*w=(b*d+k)*v','f=(b*d+k)*w','z,v,w,f'); Gf=collect(G.f,d); %group the powers of d Gf=collect(Gf,F); %group all terms multiplied times F Gf=subs(Gf,[m b k],[10 5 1000]); %substitute values for m, b, and k digits(4) % specifies the number of digits shown in answers Gf=vpa(Gf); % converts the numbers in Gf from integers to numbers with a decimal pretty(Gf) % shows Gf in standard transfer function format [num,den]=numden(Gf/F); % separates the num. and denom. and removes F from the numerator Num=sym2poly(num); % converts the symbolic num. to the standard vector for a poly Den=sym2poly(den); % converts the symbolic denom. to the standard vector for a poly H=tf(Num,Den) % generates the standard transfer function Executing this m-mile gives the following non-symbolic transfer function: H = 6.25 s^3 + 3750 s^2 + 750000 s + 5e07 ----------------------------------------------------------------------------------------------- 50 s^6 + 125 s^5 + 25088 s^4 + 30011 s^3 + 3.004e06 s^2 + 750000 s + 5e07 Note that the answer is in terms of ‘s’ instead of ‘D’. This is because transfer functions in MATLAB are expressed in terms of the Laplace operator ‘s’ instead of the ‘D’ operator. In this case, it does not matter if ‘s’ or ‘D’ is used since there are no initial conditions. Use of the Laplace transform is addressed in the next section of this notebook.14 Chapter 2. Solving Differential Equations 2.1 Analytical Solutions of Differential Equations Using the Laplace Transform Obtaining an ‘analytical’ solution to a differential equation means getting an equation for the output of interest as a function of time. Contrast this with getting a ‘numerical’ solution; in this case, only a plot is obtained of the output of interest as a function of time. In reality, it is almost always the case that a plot of the solution is desired in the end; consequently, even if an equation for the solution is obtained, a computer will be used to plot the equation. So, if only a plot is desired in the end, it might seem burdensome to first get an equation and then turn around and plot the equation when the plot could have been obtained numerically in a single step. When using a computer to obtain solutions to differential equations it is very important, when possible, to have an estimate of what the plot of the solution should look like in advance. It is a huge risk to simply believe that the computer solution is correct. Consequently, knowing how to obtain an analytical solution, even if the solution is not carried out to the very end, provides insights to what the solution should look like. For this reason, an in-depth coverage on obtaining analytical solutions is included in this document. In general, for 2nd or higher order differential equations, obtaining an analytical solution requires that the differential equation be linear. If not linear, either numerical integration is required resulting in a numerical solution or, the differential equation must be linearized if an analytical solution is desired. The following sections of this document pertain to the use of the Laplace transform to get analytical solutions to linear differential equations. 2.1.1 Linearizing Nonlinear Differential Equations This section of the notebook pertains to the use of the Laplace transform to obtain analytical solutions for ordinary linear differential equations. Just as was the case for the ‘D’ operator, the Laplace transform enables us to convert a differential equation to an algebraic equation. However, the use of the Laplace transform does not require that initial conditions be zero. To Laplace transform an ordinary differential equation, the differential equation must be linear. Quite often differential equations for engineering systems contain nonlinear functions complicating or making it impossible to get analytical solutions. If an approximate solution is acceptable, it is often possible to replace a nonlinear function within a differential equation with a straight line approximation enabling us to use the Laplace transform. The difficulty is not in getting the straight line approximation but rather determining the best straight line to use for the approximation. For example, consider the following nonlinear differential equation for x(t) with initial condition ��(0−) = 1. The minus sign on the zero, 0-, refers to time just before t = 0. Initial conditions are always found at t = 0-, not at t = 0 or at t = 0+. This is very important to remember. 3 x+ x = (0 ) =1 8 x − This equation is nonlinear because of the ����term. The input to this equation is a constant equal to 8. It is straight forward to demonstrate that the solution for x will start at the initial condition 15 of 1 and end at a constant final value since as time goes to infinity, ��̇ will become zero which allows for us to solve for ��(��)��→∞ = ��(∞) = ��. Thus, a justifiable linearization approach is obtain a straight line approximation for ���� using the two points �� = �� and �� = ��; then replace ����in the differential equation with the equation for the straight line. The result is ��3 ≈ 7�� − 6 Substituting for ��3in the differential equation gives a linear approximation for the original differential equation, i.e. − x+ (7x −6) = 8 (0 ) =1 x ��̇ + 7�� = 14 It is advised that the linearized differential equation should always be checked to confirm that the linearization points are satisfied. Setting ��̇ = 0 gives 7�� = 14 ���� �� = 2 which is correct. The initial value theorem (see the section on Laplace transforms) or actual solution ��(��) is needed to check that ��(��)��=0 = 1. So, for differential equations with solutions starting at an initial condition and ending at a different value, the best straight line to use when linearizing the nonlinear function is a line between the starting and ending values. On the other hand, it is possible for solutions of differential equations to always remain close (in the neighborhood) to certain values. For example, consider the mass spring damper system in Section 1.4.1. Suppose the input force is generated by hitting the mass with a hammer. Intuitively, the impact will cause the mass to accelerate to the left but eventually the motion will damp out with the position of the mass in the original location before the impact. Thus, the starting and ending values of the position of the mass are the same. In such cases, it may be beneficial to use the point and slope method for getting the straight line equation approximation for the nonlinear function. For example, suppose the initial condition to the previous nonlinear differential equation fluctuates but is always very close (in the neighborhood of) to the final value of 2, e.g. 0.18 ≤ ��(0−) ≤ 2.02. For such cases, a straight line approximation using the point (x = 2) and slope of the function at x = 2 is preferred, i.e. 3 = 8 ������ [����3 ���� ]��=2= 3(2)2 = 12 The result will be ��3 ≈ 12�� − 16 ����=2 It is important to confirm at this approximation for x3gives x = 2 at x = 2. Substituting this approximation into the original differential equation gives ��̇ + 12�� = 2416 Suggestions for getting straight line approximations for nonlinear equations For obtaining straight line approximations, consider the following thought process. Suppose we want a straight line approximation for a nonlinear function of x, f(x), i.e. f (x) ≈ cx + d Regardless of the approach, it is beneficial to first sketch the function to be approximated with your best guess for the location of the line. (1) The first approach using two points on the function f(x) is shown below. f (x) cx + d f x ( )2 f x ( )1 x 1 x 2 x In this case, solve two equations for the two unknown’s c and d, i.e. ��(��1) = ����1 + �� ��(��2) = ����2 + �� Example: Obtain a linear approximation for ��3 using the two point method for points y = 0 and y = 0.1. Be sure to start by sketching ��3and noting the two points for your straight line. The two equations are 0 = 0�� + �� ������ 0.13 = 0.1�� + �� Solving gives c = 0.01 and d = 0. Thus, for 0 ≤ �� ≤ 0.1, ��3 ≈ 0.01�� + 0 y 3 0.001 0 3 y ≈ y +.01 0 0.1 y 17 (2) The second approach is to use one point on the function f(x) and the slope at that point. f (x) cx + d f x ( )1 x x 1 In this case, use the truncated Taylor series approach to getting c and d, i.e. f x ≈ cx + d = f x + f x x − x ( ) ( ) '( )[ ] 1 1 1 where '( )1 c = f xis the derivative (slope) of f(x) with respect to x at 1 x = xand ( ) '( ) 1 1 1 d = f x − x f x . Example: Obtain a straight line approximation for ����for values of �� in the neighborhood of �� = ��. ��. Writing the first two terms of the Taylor series for ���� gives 2(�� − 0.1) = 0.001 + 0.03(�� − 0.1) = 0.03�� − 0.002 ��3 ≈ 0.13 + 3����=0.1 y 3 3 yy ≈ cy + d = 0.03y − 0.002 0.0010 0.1 y 18 Example: Water flows into and out of the water tank shown below. The differential equation for the height H of the water in the tank is ��̇ + 4√�� = 12 H 4 3 H H 3 9 16 9 H H Obtain a linear approximation for this differential equation by obtaining straight line approximations for √�� for the two different sets of initial conditions given below. Note, the equilibrium height ���� of the water in the tank is found by setting ��̇ = 0, ��.��. ���� = (124)2 = 9. This is the height where the flow rate out of the tank equals the flow rate into the tank. (a) Assume the initial height of the water in the tank is 16. The height will start at 16 and end at 9. So, we want to approximate √�� between the points 16 and 9 using the two point method. 4 = 16�� + �� 3 = 9�� + �� Solving for c and d gives √�� ≈ 0.14286�� + 1.71429. It is wise to check the result. For H=9, 0.14286*9+1.71429 = 3.00003. For H=16, 0.14286*16+1.71429 = 4.00005. (b) Assume the initial height of the water in the tank is between 8.99 and 9.01. Since the equilibrium height is 9, the height of the water will always be in the neighborhood of 9. Thus, the point and slope method at H = 9 will be used to approximate √��. √�� ≈ √9 +1 2√9(�� − 9) = 0.16666�� + 1.519 Checking the result: for H=9, √�� ≈ 0.16666 ∗ 9 + 1.5 = 2.99994 Using the c and d from (a) or (b), in either case, the linear approximation for the differential equation is ��̇ + 4(���� + ��) = 12 or ��̇ + 4���� = 12 − 4�� Problems utilizing this concept: 4.2, 4.25, 4.26, 5.2, 5.25, 5.26, Q2, Q8, H8, E1c, E3b, E3d 2.1.2 The Definition of the Laplace Transform Consider some general function of time ��(��). The integral definition of the Laplace transform of ��(��) is as follows: L f t F s f t e dt st ∫∞ − [ ( )] ( ) ( ) = =0 where t denotes time and s is the Laplace transform variable. Note that the integration starts at time t = 0; thus, the function f(t) prior to t = 0 does not change the resulting Laplace transform. 2.1.3 Properties of the Laplace Transform Laplace transform of derivatives considering initial conditions − L f t sF s f [ ( )] ( ) (0 ) = − 2 − − L f t s F s sf f ( ( )] ( ) (0 ) (0 ) = − − 3 2 − − L f t s F s s f sf f [ ( )] ( ) (0) (0 ) (0 ) = − − − Note that the Laplace transform of derivatives gives the same result as using the ‘D’ operator if all initial conditions are zero; the only difference is ‘s’ is used instead of ‘D’. Note also that the Laplace transforms of the derivatives contain initial condition terms at time t = 0-. f f (t) (0 )− f (0 )+ − 0+ 0time, t 0 It is important to note that a function of time may have different values at times 0+ and 0-. However, for most engineering problems, there will be no difference in the variables at t = 0-and 0+. Differences are caused by impulse inputs. For example, think of a bat hitting a baseball. At time t = 0-the speed of the ball may be 100 mph towards the batter; at time t = 0+the speed of the ball may be 200 mph away from the batter. The bat hitting the ball is effectively an impulse force input to the ball instantaneously changing the velocity of the ball. The subject of where to 20 find the initial conditions is controversial; however, initial conditions for real physical engineering systems are always found at time t = 0-regardless what you may read in another book or what someone tells you. 2.1.4 Final Value Theorem If you know the Laplace transform of a variable it can be used to find the final value of that variable using the final value theorem. That is, the final value theorem is used to compute ��(∞) which is the limit of ��(��) as t goes to infinity using the Laplace transform of ��(��), F(s). The equation for the final value theorem is f f t sF s ( ) lim ( ) ( ) ∞ =t ⎯→∞=s= 0 It is important to note that the final value theorem can only be used if none of the poles of F(s) have positive real parts; otherwise, the final value theorem will give the wrong answer. Consider the following examples: ��(��) =�� + �� ��(�� + ����)����(��)��=�� = [�� + �� �� + ����]��=��= ��. �� ���������� ���� �������������� ���������� ���� ��(��)��=∞ ��(��) =�� + �� ��(�� − ����)����(��)��=�� = [�� + �� �� − ����]��=��= −��. �� ������������������ ���������� ���� ��(��)��=∞ ��(��) =�� + ������ ��(�� + �� + ������)(�� + �� − ������)����(��)��=�� = [�� + ������ =������ (�� + �� + ������)(�� + �� − ������)]��=�� ���� + ������ = �� ���������� ���� �������������� ���������� ���� ��(��)��=∞ ��(��) =�� + ������ ��(�� − �� + ������)(�� − �� − ������)����(��)��=�� = [�� + ������ =������ (�� − �� + ������)(�� − �� − ������)]��=�� (−��)�� + ������ = �� ���������� ���� ������������������ ���������� ���� ��(��)��=∞ 2.1.5 Initial value theorem The initial value theorem is used using the Laplace transform F(s) to compute ��(0+) which is the limit of ��(��) as t approaches zero from infinity. It is important to note that (0 )+ fmay not equal ��(0−) if the input to the system contains an impulse input. +=+=s t f (0 ) lim f (t) sF(s) 0 Example: ��(��) =����+�� ��(��+��)←⎯⎯ =∞ 21 f f t sF s 0 4 + ( ) lim ( ) ( ) 0= ∞ = = == 2 ⎯→∞ s t 0 2 + 4 3 + s + f f t sF s (0 ) lim ( ) ( ) = = = 0= 3 + s ←⎯⎯ =∞ s t 1 + 2 s =∞ The final value and initial value theorems are useful for checking the accuracy of Laplace transforms of differential equations and for checking the accuracy of inverse Laplace transform solutions. Problems associated with the initial value: 4.4, 4.8, 4.21, 4.23, 5.4, 5.8, 5.20, 5.21, 5.23, 5.26, Q4, Q10, H4, H9, H18, H39, E1a, E1b, E2b, E2c, E3d, E5a, E5b, E5c, E6a, E6b, E6c 2.1.6 Converting Differential Equations to Algebraic Equations Using Laplace Transform Consider the following 2nd order linear differential equation for y with initial conditions ��(0−) ������ ��̇(0−); the input to the differential equation is ��(��). my+ cy+ ky = f Laplace transforming both sides of the equation, considering that m, c and k are constants that have not yet been assigned numerical values, gives: m s Y s − sy − y +c sY s − y + k Y s = F s 2 − − − [ ( ) (0 ) (0 )] [ ( ) (0 )] ( ) ( ) Solving for Y(s) gives: ⎢⎣⎡+ + ⎥⎦⎤ ( ) 1 ⎢⎣⎡+ + + ⎥⎦⎤ − − y ms c my (0 )( ) (0 ) + + Y s2 2 = ms cs k F s ( ) ms cs k This algebraic equation is of the form Y(s) = H(s)F(s) + I(s) where H(s)is called the transfer function relating the Laplace transform of the input F(s)to the Laplace transform of output Y(s), i.e. H s21 ( )⎢⎣⎡+ + =ms cs k ⎥⎦⎤ 22 Note that the transfer function does not contain any of the initial condition terms. Or, more importantly, if only the transfer function is needed, it doesn’t matter what the initial conditions are and the transfer function can be written simply using the coefficients of the differential equation skipping the formal process of Laplace transforming the differential equation. I(s)exists in the equation above because of the nonzero initial condition terms; if the initial conditions are zero, then I(s)is zero. F(s)remains a general expression until the input is specified. The Laplace transforms for various common inputs are presented in the next section. 2.1.7 Common Inputs and Their Laplace transforms 2.1.7.1 Step Functions and Constants One of the most common or frequently experienced inputs is a step function. f (t) A 0 time, t Using the integral definition of the Laplace transform, the Laplace transform for a step function of magnitude A is found as follows: − A [ ] [ ]sA ∞ − −∞ ∫ − A F s Ae dt st st− = ( ) 00 1 = = 0 s e = s It is important to note that a constant of A for the input has the same Laplace transform. This is because the integral for F(s) starts at t = 0; it doesn’t matter if the input was 0 or A prior to t=0. When performing simulations, however, there may be a difference in the simulation results for a step input and a constant input if the model contains input derivatives. For example, the Laplace transform of ��̇is ����(��) − ��(0−); thus, if a differential equation contains the input derivative ��̇, then the Laplace transform of the differential equation will contain ��(0−) which may or may not be zero depending on the input being either a step (��(0−) =0) or a constant (��(0−) ≠ 0). f (t) A F(s) = A s 0 time, t Problems with step and constant inputs: 4.9, 4.10, 4.13, 4.16, 4.22, 4.24, 4.25, 4.26, 4.30, 4.36, 4.37, 4.43, 5.9, 5.10, 5.13, 5.16, 5.21, 5.24, 5.25, 5.26, 5.27, 5.30, 5.36, 5.37, 5.40, 5.41, 5.43, Q11, Q12, Q15, Q18, Q21, H5, H6, H7, H9, H15, H16, H20, H21, H21, H26, H27, H32, H35, H36, H40, E2a, E2c, E3a, E3b, E3c, E3d, E5b, E5c, E6b, E6c23 2.1.7.2 Pulse and Impulse Functions An impulse function with area A is shown below. f (t) Aδ (t) 0 time, t A simple derivation for the Laplace of an impulse function is to consider a pulse with area A such as shown below f (t) A b 0 time, t b The Laplace transform for this pulse is obtained using the integral definition, i.e. A − A − A s t bs te ∞ − − ( ) 0 [ ] [ 1] 0 ∫ ∫b s b s t F s = + =− − 0− b e dt e dt b bs e = bs The limit of f (t)as b→0 is Aδ (t), an infinitely high function with a zero width; the area is still equal to A regardless of the size of b. The limit of F(s) as b→0 is found using L’Hopital’s rule to be A, i.e. ��(��)��=0 = (−��[��−���� − 1] ���� )��=0= �� Thus, the Laplace of Aδ (t)is A; if A=1, f (t)is called a unit impulse, δ (t). Problems with impulse inputs: 4.4, 4.6, 4.8, 4.24, 4.26, 4.41, 5.6, 5.8, 5.24, 5.25, 5.26, Q7, Q10, H4, H5, H9, H10, H11, H12, H15, H16, H19, H20, H21, H27, H36, H38, E3a, E3b, E3d, E4c, E5c, E6c.24 2.1.7.3 Exponential Functions The Laplace transform for an exponential input, denoted by bt ae−,and shown below is a ( ) = . F s+ s b f (t) a ae− bt 0 time, t 2.1.7.4 Periodic Functions Periodic functions repeat themselves in a finite period of time. Equations can be obtained for periodic functions using the Fourier series which is an infinite sum of sine and cosine terms. See sections 2.2.4 and 2.4.1. If ��(��) = Asin(����) ��ℎ���� ��(��) =���� ��2+��2. If ��(��) = Acos(����) ��ℎ���� ��(��) =���� ��2+��2 2.1.8 Inverse Laplace Transform 2.1.8.1 Inverse Laplace Transform Using Partial Fractions The use of partial fractions to get an inverse Laplace transform is preferred if a computer is available to factor the denominator and then compute the partial fractions. Note, if the denominator is not already factored, then a computer will definitely be required to factor any polynomial 3rd or higher order. For cases where the denominator is already factored, a much more efficient method for getting the inverse Laplace is by means of the residue theorem which will be explained in section 2.1.8.2. Any Laplace transform can be expressed as the sum of its partial fractions. If the numerator order is not greater than the denominator order and there are no repeated poles, each partial fraction will be one of three possible forms. In general, the partial fractions for some general Laplace transform H(s) consists of a constant K, plus 2nd order partial fractions for complex poles, and 1st order partial fractions for real poles, i.e. ⎢⎣⎡+ ⎥⎦⎤ a s b + a s b + C C H s Kω ω ( ) = + ... ... 1 1 2 2 + 2 2 2 2 + + 1 s p + + 2 s p ( ) ( ) s r + + + 1 1 s r + + 2 2 1 2 The constant K only occurs if the order of the numerator polynomial equals the order of the denominator polynomial. 25 For example, expressing the following F(s) in partial fractions gives the three possible partial fraction formats: 4 3 2 2 37 214 550 375 ( )4 3 2 2 s s s s + + + + 3 4 8 7 F s s + =s s = + 2 + + s s s s s s + + + + 11 55 125 0 + 0 + 5 + + 6 25 Based on the discussion pertaining to inputs in section 2.1.2.6, the first partial fraction, 2, corresponds to an impulse with area (or magnitude) of 2. The second partial fraction corresponds to a constant (or step) with magnitude of 3. The third partial fraction corresponds to an exponential term, t 4−. The fourth term can be shown to correspond to a sine wave with an e5 exponentially decaying amplitude. The general form for a 2nd order partial fraction with complex poles is denoted by as b + 2 2 ( + ) +ω s d and its inverse can be shown to be − dt, ce t sin(ω +φ) where �� = √(��−���� ��)2 + ��2 and ∅ = ����������[(�� − ����) + ������]. A sine wave with an exponentially decaying amplitude is shown below f (t) − dt ce t sin(ω +φ) 0 time, t 8 7 s is 9.0588 sin(4 2.0591). 3+ + − t Thus, using So, the inverse Laplace of 6 25 2+ + s s e t superposition, the total inverse is the sum of the individual partial fraction inverses, i.e. − − 5 3 t t f t t e e t ( ) 2 ( ) 3 4 9.0588 sin(4 2.0591) = + + + + δ26 So, if a Laplace transform is expressed in partial fraction format, then these four formulas can be used to get the inverse Laplace transform directly without a computer or tables. Consider a second example. 2 3 ⎢⎣⎡+ = + + ⎥⎦⎤ − − − 1 1 4 t L F s L [ ( )] 5 = + + δ 5 ( ) 2 3 s s 4 t e Consider the mass/spring/damper system shown below. This system was also introduced in section 1.4.2; however, in 1.4.2, b was equal to 5 and the output of interest was the ground force f(t). Changing the output of interest to z(t), it can be shown that the Laplace transform Z(s) for the displacement z(t) for a unit impulse input is given by ��(��) = [100��4 + 4��3��3 + 7��4��2 + 6��5�� + 3��6 1��3��6 + 5��4��5 + 1.1��6��4 + 1.3��7��3 + 9��7��4��2 + 3��8�� + 1��9] F(t) M=10 kg k k k M M M b=100 Ns/m z k=1000 N/m b v b w b This Laplace transform in partial fraction format can be shown to be ⎢⎣⎡+ + ⎥⎦⎤ 0.05431 ( )s s s Z s 0.03493 0.01076 =2 2 2 2 2 2 ( 0.991) 4.34 + + + ( 7.77) 9.75 + + + ( 16.2) 7.81 Using the formula for complex partial fractions above, the inverse Laplace of Z(s) can be shown to be 0.991 7.7 7 1 6.2 − t − t − t z t e t e t e t ( ) 0.0125 sin(4.34 ) 0.0036 sin(9.75 ) 0.0014 sin(7.81 ) = + + Some of these examples demonstrate the need to be able to use a computer to compute the partial fractions if partial fractions are desired. The next section introduces a homemade MATLAB m file ‘pfract’ for computing partial fractions; ‘pfract’ is listed in Chapter 13. Problems involving partial fractions: H12, H36, E1a27 Obtaining Partial Fractions Using MATLAB In general, any Laplace transform can be expressed in partial fraction format. Although it is possible to obtain the partial fractions manually, a computer will be required to factor the denominator polynomial if the order is greater than 2. Since a computer is to be used to factor the denominator, it might as well be used at the same time to obtain the partial fractions. An m file named ‘pfract’ for computing partial fractions has been written and is included with this document. This program factors the denominator and then computes the partial fraction numerator coefficients. The program is limited to first order (non-repeated) real or complex poles. For example, consider the following Laplace transform H(s): ��(��) =10��3 + 6��2 + 28�� + 150 2��3 + 16��2 + 74�� + 100 The MATLAB commands for getting the partial fractions and the results are >> num=[10 6 28 150]; >> den=[2 16 74 100]; >>pfract(num,den) % ‘pfract’ is not a MATLAB M-file but must be copied from Chapter 13 and saved. ��(��) =10��3 + 6��2 + 28�� + 150 2��3 + 16��2 + 74�� + 100 = 5 +1.118 �� + 2+−38.12�� − 101.5 Or >> [R]=pfract(num,den) ��(1) = 5 ��(2) =1.118 �� + 2 ��(3) =−38.12�� − 101.5 ��2 + 6�� + 25 ��2 + 6�� + 25 Note that one of each of the three possible partial fraction formats occur in this example. The poles are -2 and −3 ± ��4. The partial fraction ‘5’ occurs because both the numerator and denominator are 3rd order. 2.1.8.2 Inverse Laplace Transform Using the Residue Theorem The residue theorem is the simplest way to get the inverse Laplace transform if you do not already have the function in partial fraction format and if the denominator is already factored; the numerator order must be at least one less than the denominator order. For cases where the numerator order equals the order of the denominator, long division must first be used to get the appropriate format for the residue theorem. The formula for the inverse Laplace transform of a general function F(s) is as follows:28 ��(��) = ��−1{��(��)} = ∑���������������� ���� ��(��)���������� ��ℎ�� ���������� ���� ��(��) To get the residue for a first order real pole ‘–p’ where p ≥ 0: �������������� ������ ��ℎ�� �������� ′ − ��′ ���� (�� + ��)��������(��)��=−�� For example: ⎢⎣⎡+ + = + ⎥⎦⎤ s t s t − L 2 1 ( 3)( 5) s s ⎢⎣⎡+ + 2 ( 3) s e ( 3)( 5) s s + + ⎥⎦⎤ ⎢⎣⎡+ + ⎥⎦⎤ 2 ( 5) s e ( 3)( 5) s − − 3 5 t t = − e e s s =− =− 3 5 s Consider the following Laplace transform with a 1st order pole, a 2nd order pole and a 3rd order pole. ��(��) =5 (�� + 2)(�� + 4)2(�� + 10)3 The pole -2 is first order; the pole -4 is second order; and the pole -10 is third order. To get the residue for an nth order pole –p in the denominator of F(s), the formula is (�� − 1)![����−1[(�� + ��)����(��)������] �������������� ������ ���� ����ℎ ���������� �������� − �� =1 ������−1]��=−�� An example demonstrating the procedure for a 2nd order pole is as follows: ⎢⎣⎡+ + 2 s t ⎤ ⎢⎣⎡+ + ⎤ s t ⎢⎣⎡+ + ⎥⎦⎤ − L 1 ⎥ = + ( 3) s e 2 ⎥ + 1 d 2 2 ( 5) ( 3)( 5) s s 2 2 ( 3)( 5) s s ⎦ ⎦ s (2 1)! − ds s e + ( 3)( 5) s s 2 =− =− 3 s 5 − − − 3 5 5 t t t = − − 0.5 0.5 e te e Note, it is advised to write the denominator in factored format before beginning the residue process to avoid errors. For example 4��2 + 28�� + 40 =8�� + 16 8�� + 16 4(�� + 2)(�� + 5)=2�� + 4 (�� + 2)(�� + 5) Shortcut for Complex Poles, − r ± jω Complex Poles With Nonzero Real Parts, r ≠ 0 For Laplace transforms F(s) with complex poles, the two residues for a complex pair will combine and simplify to an exponentially decaying sine wave. Realizing this end result, enables one to use a shortcut and obtain the resulting combined solution in a single step instead of the two steps that would be required for each of the two complex conjugates. Consider the following Laplace transform with a pair of complex poles −�� ± ����:29 N s ( ) F s+ + ω + − ω ( )s r j s r j D s = ( )( ) ( ) where N(s) is the numerator of F(s) and D(s) is what remains of the denominator of F(s) after removing the complex pair that we are working with. The sum of the residues for the complex pair of poles can be shown to simplify to the following expression which we will use for the short cut: ert − [residue for − r + jω] + [residue for − r − jω]=sin(ω φ) ω+ M t where = + + =− +=( )( ) N s 2 2 M s r F s [( ) ] ( ) ω ω and s r jD s s r j =− + ω ⎪⎨⎧⎟⎟⎠⎞ ⎪⎬⎫ ⎜⎜⎝⎛ angle s r F s angle( )( ) 2 2 φ ω ω = + + = {[( ) ] ( ) } N s s r j D s =− + ⎪⎩ ω s r j =− + ⎪⎭ For an example, consider the following Laplace transform s + 5 N s ( ) Y s+ + + − ( )D s s j s j = s s j s j ( 2 3)( 2 3) + + + − = ( )( 2 3)( 2 3) Thus, for this problem, N(s) = s+5 and D(s) = s. The inverse Laplace of Y(s) consists of the sum of three residues, i.e. y(t) =[residue_ for _ s = 0]+[residue_ for _ s = −2− j3]+[residue_ for _ s = −2+ j3] But, we will use the short cut for the sum of the residues for the complex pair, i.e. et − 2 y(t) =[residue _ for _ s = 0]+ sin(3 ) e y t 5 − 2 3 t M t +φ ( ) = + +φ M t 13 3 �� = |�� + �� sin(3 ) ��|��=−��+����= ��. �������� ∅ = ������−�� ([�� + �� ��]��=−��+����) = −��. �������� ������.30 Complex Poles with Zero Real Parts, r = 0 If any of the poles are complex with a zero real part, then it is important to note that the exponential term ��−���� = ��−���� = ��; thus, the multiplier on the sine function does not go to zero as �� → ∞. This will always happen if the input to a transfer function is a sine or cosine wave since the denominator of the Laplace transform of sine and cosine functions have a zero real part. Consider the previous example but with the real part of the complex poles equal to zero. =s ss Y s s + 5 + 5 ( )2 2 s s j s j ( 0 3)( 0 3) + + + − = ( 3 ) + et − 0 y(t) =[residue _ for _ s = 0]+ sin(3 ) e y t 5 − 0 3 t M t +φ ( ) = + + = + − M t t sin(3 ) 0.556 0.6479sin(3 1.03) 2 3 3 �� = |�� + �� φ ��|��=−��+����= ��. �������� ∅ = ������−�� ([�� + �� ��]��=−��+����) = ������−������−����= −��. ���� ������. For cases of a sine function input to a transfer function, it is not uncommon to be interested in the inverse Laplace transform of the output after the transient portion of the solution has gone to zero or only for large values of time (t→ ∞). If the poles of the transfer function all have negative real parts, then the portion of the inverse Laplace transform associated with the poles of the transfer function will all go to zero as t → ∞. Thus, the only thing left will be the portion of the inverse associated with the input poles which will have zero real part which means this portion does not go to zero but remains a sine or cosine function forever. For example, suppose we have ��(��) = ��(��)��(��) = [��������(�� + ��) (�� + ��)(�� + ����)[(�� + ��)�� + ����]]��(��) where the input ��(��) = ��������(������) ������ ��(��) =������ ����+������. The transfer function G(s) is the expression in the brackets; note that all of the poles of G(s) have negative nonzero real parts. Thus, the portions of the inverse Laplace associated with these poles will go to zero as t → ∞. So, as t → ∞, the only portion left will be the portion of the inverse associated with the poles of the input R(s). Thus, ��(��)�� → ∞ = �� + �� + �� +������|��(��)|��=��+������������ ������(������ + ��) = ��|��(������)| ������(������ + ��) where �� = ���������� ���� ��(������). Thus, if a transfer function has a sine input with frequency 20 rad/sec and amplitude 30, that is r(t)=30sin(20t), then as t → ∞, the output will be a sine function with the same frequency 20 rad/sec and with an amplitude ����|��(������)| and with a phase shift of �� = ���������� ���� ��(������). That is, 31 |��(������)| = |��������(������ + ��) (������ + ��)(������ + ����)[(������ + ��)�� + ����]| = |−��. �������� + ����. ��������| = ��. ������ �� = ����������(−��. �������� + ����. ��������) = ��. ���� ������ Thus, if r(t)=30sin(20t), then as t → ∞, ��(��) = ���� ∗ ��. ������ ∗ ������(������ + ��. ����) = ��. ����������(������ + ��. ����) The main points to be learned are that for a system (transfer function) with a sinusoidal input, the output will converge to a sinusoidal function at the same frequency as the input with the amplitude of the output equal to the input amplitude multiplied by the magnitude of the transfer function using the input frequency; there will be a phase difference in the input and output calculated using the input frequency. This also applies to cosine inputs. Example The following system has r(t) = 4sin(0.3t) for the input. What will be the output y(t) at steady state conditions (after all transients have died out)? How much greater or smaller is the output amplitude compared to the input amplitude? ��(��) = [�������� + ��. ������ + �� �������� + �������� + ���� + ��]��(��) Ans: y(t)=2*4*sin(0.3t-1.35)=8sin(0.3t-1.35) Note the output amplitude is 2 times greater than the input amplitude. Problems containing inverse Laplace transforms: 4.6, 4.23, 4.27, 4.30, 4.35, 4.36, 4.37, 4.43, 5.5, 5.6, 5.9, 5.22, 5.23, 5.25, 5.26, 5.27, 5.30, 5.35, 5.36, 5.40, 5.43, Q5, Q6, Q8, Q12, H2, H3, H4, H6, H10, H12, H13, H19, H26, H27, H37, H38, H39, E2a, E2b, E2c, E3a, E3b, E3c, E3d, E4B, E4c,E5c,E6C 2.1.8.3 Using MATLAB to Plot the Inverse Laplace Transform Without Getting an Equation for the Inverse Laplace Transform The previous sections of this document have focused on using the Laplace transform to get an equation for the solution to a linear ordinary differential equation using the inverse Laplace transform. This section introduces how to use the Laplace transform and MATLAB to get a plot of the solution without having to first get an equation for the solution. For cases where only a plot of the solution is desired, this approach to getting a plot of the solution represents the ultimate in simplicity since there is minimum chance for error and minimum effort required to get a plot. Also, this method can be used for cases with nonzero initial conditions and/or unusual inputs; the Laplace of the input must be known. Suppose you have a Laplace transform Y(s). If an equation for the inverse Laplace, y(t), is desired, then the only options are to use the partial fraction or residue methods mentioned in 32 previous sections and use inverse Laplace to get the equation. If however, an equation is not needed but only a plot of y(t), then it is very convenient to use the ‘impulse’ command in MATLAB to get the plot directly from the Laplace transform. If you pretend that Y(s) is the transfer function for a system with a unit impulse input, then Y(s)*1 = Y(s) is the Laplace of the output. Thus, entering Y(s) into MATLAB as a transfer function (pretending that it is a transfer function) and then using the command ‘impulse(Y)’ gives a plot of y(t). For example, to get a plot of the inverse of Y s ( )3 2 2 =s s s + + + 13 55 75 >> Y=tf(2,[1 13 55 75]) %This MATLAB command defines the transfer function above. >> impulse(Y) %This MATLAB command generates a plot of y(t) A more versatile use of the impulse command allowing labeling and line width and color control is as follows: >> Y=tf(2,[1 13 55 75]); >> [y,t]=impulse(Y,3); % Instead of generating a plot, this generates values for y and t. >> plot(t,y,'r','LineWidth',3) >> ylabel('y') >> xlabel('time, sec.') >> title('Using the impulse command to plot the inverse Laplace of Y(s)') Using the impulse command to plot the inverse Laplace of Y(s)0.03 0.025 0.02 0.015 y 0.01 0.005 0 0 0.5 1 1.5 2 2.5 3 time, sec. 33 For a second example, consider the problem in section 2.1.6, i.e. ⎢⎣⎡+ + ⎥⎦⎤ ( ) 1 ⎢⎣⎡+ + + ⎥⎦⎤ − − y ms c my (0 )( ) (0 ) + + Y s2 2 = ms cs k F s ( ) ms cs k If the input is a unit step, then F(s)=1/s; Y(s) becomes ⎢⎣⎡+ + ⎤ − − − − 1 1 (0 )( ) (0 ) [ (0 )( ) (0 )] 1 y ms c my ⎢⎣⎡+ + + + + ⎥⎦⎤ y ms c my s + + + Y s+ + ( )2 2 2 = ms cs k s ms cs k ⎥ = ⎦ s ms cs k ( ) If numerical values are assigned to the initial conditions and the parameters m, c, and k and then the final expression for Y(s) is entered as transfer function (pretend), then the command ‘impulse(Y)’ gives a plot of the inverse of Y(s) which is y(t). Thus, for m=10, c=0.05, k=50, y(0-)=6, and ��̇(0−) = 2: ��(10��2 + 0.05�� + 50)=60��2 + 20.3�� + 1 ��(��) =[6(10�� + 0.05) + 10(2)]�� + 1 10��3 + 0.05��2 + 50�� + 0 >> Y=tf([60 20.3 1],[10 0.05 50 0]); % Pretending that Y(s) is a transfer function. >> impulse(Y) This gives us a plot of the solution to the original differential equation even though the input was not an impulse but rather a step input and for a case of nonzero initial conditions. Problems utilizing this concept: H4, H5, H10, H11, H12, H15, H16, H19, H20, H27, H36, H38 2.1.9 Modes of a System The initial condition or impulse response of a dynamic system is madeup of the sum of the responses of the individual modes of the system. It is not uncommon for some of the modes to dominate the response of the system; consequently, relatively high order system models are sometimes approximated by lower order models consisting of only the dominate modes. The way to determine the individual modes is to express the Laplace transform in partial fractions. Each partial fraction represents a mode. If all initial conditions are zero and the input is an impulse, the modes are simply the partial fractions of the transfer function of the system. Modes are either be first order or second order. The dominant mode usually has the largest eigenvalue time constant. A complex mode with a small damping ratio may dominate the other modes depending on the time constants. The procedure for computing the modes is to use the MATLAB m-file ‘pfract’ on the numerator and denominator polynomials of the Laplace transform or of the system transfer function. Consider the mass/spring/damper system introduced in section 2.1.8.1. 34 Suppose that the output of interest is the displacement z(t) and the input F(t) is a unit impulse. It can be shown that the transfer function H(s) relating Z(s) to F(s) is the following: ��(��) = ��(��)��(��) = [100��4 + 4��3��3 + 7��4��2 + 6��5�� + 3��6 1��3��6 + 5��4��5 + 1.1��6��4 + 1.3��7��3 + 9��7��4��2 + 3��8�� + 1��9] 1 The following MATLAB commands generate the individual modes of this system: >> num=[100 4000 70000 600000 3e6]; >> den=[1000 5e4 1.1e6 1.3e7 9e7 3e8 1e9]; >> pfract(num,den) 0.05431 ( )2 2 2 0.03493 0.01076 Z s ⎢⎣⎡+ + ⎥⎦⎤ =32.47 324.7 + + s s s s s s + + 1.981 19.81 + + 15.55 155.5 To determine if any of these modes dominate the response, the ‘damp’ command is used to generate the eigenvalues, i.e. >> damp(H)
Eigenvalues
Damping Ratio
Time Constant, s
Undamped Natural Frequency, rad/s
−0.99 ± ��4.34
0.223
1.01
4.45
-7.77±��9.75
0.623
0.128
12.5
-16.2±7.82
0.901
0.062
18.0

## What will be the output y(t) at steady state conditions (after all transients have died out)?

## (1) What are the natural variables of this system?

## How much greater or smaller is the output amplitude compared to the input amplitude?

Don't forget about the age old question of corey johnson unc

Don't forget about the age old question of college psychology notes

Don't forget about the age old question of 3e8/1e9

Don't forget about the age old question of research on social loafing demonstrated that european participants worked harder on a task when working alone, whereas chinese participants worked harder on a task when they were part of a group. these results illustrate the importance of the _____ perspe

We also discuss several other topics like cc 302

Don't forget about the age old question of draw the shear diagram for the shaft.

References: