MAE 384. Advanced Mathematical Methods for Engineers. Programming Assignment 3. 1. (25 pts.) Consider the differential equation: 𝑑𝑦 1 = − 𝑦 sin2 𝑡 𝑑𝑡 2 with initial condition given by 𝑦(0) = 1 Solve this equation from t = 0 to t = 8π using the following methods: (a) (b) (c) Solve analytically by separating variables and integrating. Solve using the 4th-order Runge-Kutta method (write your own code for this, do not use the MATLAB provided ODE solvers) for the following two step sizes: I. Maximum step size for stability (don’t try and do this analytically – try out your code for different step sizes to find the stability limit). II. Maximum step size for a time-accurate solution. “Good” accuracy can be defined in several ways, but use the definition that the numerical solution remains within 2% of the true solution a t = nπ. Solve using the MATLAB function ode45. Present your results for the above problem in an appropriate fashion. For problem 1, be sure to include a comparison of the numerical methods with each other and with the true solution. Be sure to discuss your findings with respect to the notions of stability and accuracy of the numerical methods. 2. (25 pts.) A car and its suspension system traveling over a bumpy road can be modeled as a mass/spring/damper system. In this model, 𝑦1 is the vertical motion of the wheel center of mass, 𝑦2 is the vertical motion of the car chassis, and 𝑦0 represents the displacement of the bottom of the tire due to the variation in the road surface. Applying Newton’s law to the two masses yields a system of second-order equations: 𝑚1 𝑦̈ 1 + 𝑐2 (𝑦̇ 1 − 𝑦̇ 2 ) + 𝑘2 (𝑦1 − 𝑦2 ) + 𝑘1 𝑦1 = 𝑘1 𝑦0 𝑚2 𝑦̈ 2 − 𝑐2 (𝑦̇ 1 − 𝑦̇ 2 ) − 𝑘2 (𝑦1 − 𝑦2 ) + 𝑘1 𝑦1 = 0 (a) Convert the two second-order ODE’s into a system of 4 first-order ODE’s. Write them in standard “state-space” form. (b) Assume the car hits a large pothole at t = 0 so that −0.2 m 0 ≤ 𝑡 < 0.2 s 𝑦0 (𝑡) = { 0 𝑡 > 0.2 s Create a MATLAB function that returns the right hand sides of the state-space equations for an input t and an input state vector. (c) Solve the system on the time interval [0 60] seconds using the MATLAB function ode45. Find the displacement and velocity of the chassis and the wheel as a function of time. Use the following data: 𝑚1 = 90 kg, 𝑚2 = 1700 kg, 𝑘1 = 140 N/m, 𝑘2 = 30 N/m, 𝑐1 = 145 N-s/m 3. (25 pts.) Write a MATLAB program to simulate the dynamics of a helicopter lifting a survivor. When lifting the survivor into the helicopter with a constant speed winch, the resulting dynamics are non-linear, and stability is dependent upon the winch speed. Using polar coordinates, we can find the equations of motion to be: −𝑚𝑔 sin 𝜃 = 𝑚(𝑟𝜃̈ + 2𝑟̇ 𝜃̇) 𝑟̇ = constant (negative) Notice that the mass of the survivor factors out and thus the solution is independent of the mass of the person being lifted. In these equations, r is the instantaneous length of the winch cable, g, is the gravitational constant, and θ is the angle of the swing. You may choose to use either your Runge-Kutta solver from problem 1 or ode45 to integrate the equations of motion. This problem is of particular interest to the survivor since an unstable condition can cause the angle of the swing to exceed 90⁰, essentially placing him/her in danger of being beheaded by the rotor blades of the rescue helicopter. Also, it is desirable to retrieve the survivor as fast as possible to get away from the danger. Use your program to determine the maximum winch speed for which the survivor will not swing above the helicopter attach point for a lift from the initial conditions: 𝜃0 = 0.15 𝑟𝑎𝑑 𝜃0̇ = 0 𝑟0 = 40 𝑚 And ending when 𝑟 = 0.5 𝑚. The maximum lifting speed of the winch is 5 m/s. Create a plot of the stable winch speeds by plotting a 1 for each stable condition and 0 for an unstable condition (so the stability condition on the y-axis and the winch speed on the x-axis). 4. (50 pts.) Heat transfer through a long slender rod is governed by the equation 𝑑2 𝑇 + ℎ ′ (𝑇𝑎 − 𝑇) = 0 𝑑𝑥 2 Where 𝑇𝑎 is the ambient temperature and ℎ ′ is the convective heat transfer coefficient. The convective heat transfer coefficient depends on a variety of different things including velocity changes in the air passing over the bar, the materials of the bar, even the orientation of the bar. Say that a slender bar is placed between two large plates spaced 0.25 m apart. One plate is held at a constant 250⁰C and the other is held at 30⁰C. Air is blown between the gap (across the bar) between the plates. Because of the profile of the air velocity between the plates, the convective heat transfer coefficient varies across the rod according to the equation ℎ ′ (𝑥) = (−0.64𝑥 2 + 0.16𝑥) × 104 Where x is the distance along the bar starting from the hot side and going to the cold side. If the ambient temperature is 22⁰C, solve the governing equation using (a) Shooting method (you may use ode45 to propagate solutions) (b) Finite difference method (you should write a MATLAB code that will setup the coefficient matrix for you and use the backslash operator to solve this linear set of equations.). You should do numerical experiments by changing the step size in order to determine an appropriate step size. You should plot the relative error against step size on a log-log scale (you may calculate the error relative to the ode45 solution using the interp1 function using the spline option to get the ode45 solution at the same grid points as your finite difference method). Discuss the error plot and what would be an appropriate step size for this computation in your results section. (c) Plot the final solutions for parts (a) and (b) and compare the results. If you receive at least 70% of the points on this assignment, you will receive credit for outcomes #5 and #6.