Winter 2017 Lab 2 1 Mechanical Vibrations The exercises in this lab explore the use of MATLAB to approximate second order nonlinear equations and analyze mechanical vibrations. To complete this lab, 1. Open a new script in MATLAB titled Lab_1_Firstname_Lastname.m , using your name. (For example, if your name is Stuart Dent, please title the file Lab_1_Stuart_Dent.m .) 2. Write the lab number and your name as a comment at the top of the script: % Lab 1, Firstname Lastname. 3. For each exercise, write the comment % Exercise N , where 'N' is replaced by the number of that exercise. If your answer includes a graph, precede the plot command by figure('Name','Exercise N','NumberTitle','off'); , where again N is replaced by the number of the exercise, but the command is otherwise entered exactly as shown. This command tells MATLAB that the name of the graph should be 'Exercise N' and that the window does not need to be numbered. 4. Include working commands and any explanatory comments to answer the questions. Your answers should appear as the output of the script. 5. Be sure that your commands work as you intended! Before you submit your assignment, clear your workspace and run your script to test it. Introduction We cannot draw a slope field for differential equations of order 2 or higher (unless we impose some conditions on the solutions), so a graphical analysis of such an equation will be more limited than for first order equations. MATLAB can find explicit solutions to some higher order equations using dsolve, but unless the equation has constant coefficients or a very particular form, solutions cannot be written explicitly using elementary functions. We can also use MATLAB to find numerical approximations to higher order equations. Let's start with one we know how to solve explicitly: y + 3y + 2y = 0 >> syms y(t); >> dsolve(diff(y, 2) + 3*diff(y, 1) + 2*y == 0) We can give partial initial conditions: >> dsolve(diff(y, 2) + 3*diff(y, 1) + 2*y == 0, y(0) == 2) Winter 2017 Lab 2 2 or complete initial conditions: >> Dy = diff(y, t); >> dsolve(diff(y, 2) + 3*diff(y, 1) + 2*y == 0, y(0) == 2, Dy(0) == 1) Linear equations with variable coefficients are typically unsolvable, or have solutions expressed using special functions defined for this purpose: >> dsolve(diff(y, 2) + t*y == 0) and nonlinear equations are worse: >> dsolve(diff(y, 2) + diff(y,1)/y == 0) To approximate a higher order equation, MATLAB requires that we rewrite it as a linear system of first order equations. First, we replace each derivative of y with a new variable. y1 = y, y2 = y, y3 = y, ... yn = y(n1) Then each yi is the derivative of yi1, so we have n 1 equations: y2 = y, y3 = y, ... yn = y(n1) and the original differential equation can be written in terms of y1, ..., yn and the first derivative of yn (which is the n-th derivative of y): y + 3y + 2y = 0 y2 + 3y2 + 2y1 = 0, y2 = y1 MATLAB can read this as a vector: >> [V] = odeToVectorField(diff(y,2) + 3*diff(y, 1) + 2*y == 0) Notice that the entries of V are the expressions for y1 and y2. However, MATLAB's numerical solvers do not like symbolic functions like y(t). Before proceeding, we have to convert the system to a MATLAB function: >> M = matlabFunction(V, 'vars', {'t','Y'}) Now we can ask for a solution, where the second input to ode45 is the range of values for t, and the third input gives the initial conditions for y(t0) and y(t0), in this case y(0) = 2 and y(0) = 1: Winter 2017 Lab 2 3 >> sol = ode45(M, [0 20], [2 1]) This looks highly unsatisfying. MATLAB has combined the results of the numerical solver with some data about how the approximations were obtained. We can extract some information piece by piece. Here are the t-coordinates of the points: >> sol.x And here are the values of y1 = y and y2 = y computed for each t: >> sol.y We can plot this as usual, extracting the coordinates for t and y: >> plot(sol.x,sol.y(1,1:29)) However, to plot this way, we have to check first how many entries sol.x has: >> size(sol.x) Instead, we can pass the entire solution object with its data to the fplot command. >> hold on >> fplot(@(t) deval(sol,t,1), [0, 20]) Plotted this way, MATLAB recomputes the solution at 100 points instead of 29, and we see a slightly smoother graph. We can compare this to the exact solution: >> fplot(@(t) 5*exp(-t) - 3*exp(-2*t), [0, 20]) These will be indistinguishable. Winter 2017 Lab 2 4 Exercises 1. a. Use the command dsolve to solve the differential equations below. i. y + 3y = 0 ii. y + 3y = 0 iii. y + 3y = 0 iv. y + 3y = 0 v. y + 3y = 0 b. Graph each of the solutions in a. in the same window with 0 t 10. c. Which of these equations does not represent a mechanical vibration? Why not? In your comments, explain how to recognize that the equation cannot describe a mechanical vibration i. from the graph of the solution, and also ii. from the coefficients of the original equation. (Recall how we have interpreted the coefficients m, b, and k where my + by + ky = 0.) d. In your comments, classify the other four solutions as undamped, underdamped, critically damped, or overdamped. 2. Consider the differential equations below, where y(0) = 1 and y(0) = 1. i. y + 2y + 0.99y = 0 ii. y + 2y + y = 0 iii. y + 2y + 1.01y = 0 a. How do you predict the solutions will behave? b. Find the formula and graphs for each solution on the interval 0 t 10. c. Is it easier to compare their behaviors by studying the symbolic expressions of the differential equations or their solutions? 3. a. In one window, graph four different solutions to y + 10 y + y = sin t by using different initial conditions. (Be sure that all four graphs are clearly visible in the window.) b. In your comments, describe the apparent behavior of the solutions as t . 4. a. Graph solutions to y + ay = sin 3t, y(0) = 1, y(0) = 1 for each of the values a = 0.09, 0.9, 9, and 90. b. Which equation (if any) shows resonance? (Hint: You may need to change the range for x to see the full behavior of each solution.) 1 syms ytt) : 2 dy= diff{y, t); 3 dsolve{diff{y, 2) + 3*diff{y, 1) == 0, y{0) == 2, dy{0) == 1) 4 5 6 fplot{@{t) 7/3 exp{3*t)/3, [0, 10]),title{'solution Plot') 7 xlabelI't') 8 ylabel{'gjt1 Command Window }} Untitled Error: File: Untitled.m Line: 8 Column: 8 Character vector is not terminated properly. }) Untitled Error: File: Untitled.m Line: 8 Column: 8 Character vector is not terminated properly