Slide 1Slide 2Slide 3Slide 4Slide 5Slide 6Slide 7Slide 8Slide 9Slide 10Slide 11Slide 12Slide 13Slide 14Slide 15Slide 16Slide 17Slide 18Slide 19Slide 20Slide 21Slide 22Slide 23Slide 24Slide 256&7-2Dynamics, two examples of the use of ode45.Numeric Integration using trapz and quad/quadl functions.Readings: Matlab by Pratap Chapter 5.4,5.56&7-31. Problem Definition Use Matlab to plot the velocity of a free-falling object. Assume that the object is near the earth’s surface, so that the force due to gravity is given by mass * g where 9.8 m/s2. Further assume that air friction is present and that the force due to air friction satisfies Fair friction = b * v2 , where b is constant and v is the velocity of the falling object (downward is negative velocity).2. Refine, Generalize, Decompose the problem definition(i.e., identify sub-problems, I/O, etc.) Mass of object is 70Kg , acceleration due to gravity is -9.8 m/s2 and the coefficient of air friction b = .25 Kg/s. At time = 0 assume v=0, that is,v(0) = 0 (initial condition)6&7-4mb v2-mgFrom Newton’s 2nd law, the sum of the forces on the falling object equal it’s mass times acceleration:2. (continued))(**)(*2tvbgmdttdvm )(*)/()(2tvmbgdttdv6&7-52. (continued)Using calculus we can rewrite the differential equation in the form:mbv2-mgand integration of both sides gives,dttvmbgdv )(*2bmggtbmgtv tanh)(6&7-62. (continued)Since we can write tanh(x) as:mbv2-mgwhen x -> infinity tanh(x) -> 1 and so the falling body has a terminal velocity of:xxxxxxeeeeeex2211)tanh(bmgv )(6&7-73. Develop Algorithm (processing steps to solve problem)[t, v] = ode45(@aux_function, time, initial_condition);Rather than using calculus to solve this problem we can use a built-in Matlab function… ode45. The format of this function is:aux_function: user created function that computes the derivate time: a vector [start, stop] for the range of time of the solutioninitial_condition: initial value v(start)t: solution column vector of time values on range [start, stop]v: solution column vector velocity values at times t6&7-84. Write the “Function" (Code)(instruction sequence to be carried out by the computer)Use the Matlab editor to create a file vprime.m .function vp = vprime(t,v)% function vp = vprime(t,v)% compute dv/dtm = 70;g = 9.8;b = .25;vp = -g +(b/m)*v.^2;)(*)/()(2tvmbgdttdv6&7-95. Test and Debug the Code To test our program in Matlab, we will plot velocity versus time for t=0 seconds to t=10 seconds, plot(t,v)6. Run CodeTo run our program we will use ‘ode45’ as follows:>> [t, v] = ode45(@vprime, [0,30], 0);>> plot(t,v) (See plot on next slide.) [t, v] = ode45(@aux_function, time, initial_condition);As time increases the velocity levels out around -52.3818 m/s. This velocity is called the terminal velocity of the object. This agrees with the solution since v(infinity) = -sqrt(m*g/b) = -52.38326&7-111. Problem Definition Use Matlab to plot both the velocity and position(height) of a free-falling object. 2. Refine, Generalize, Decompose the problem definition(i.e., identify sub-problems, I/O, etc.) Mass of object is 70Kg , acceleration due to gravity is 9.8 m/s2 and the coefficient of air friction b = .25 Kg/s. Initial conditions: at time = 0 assume v=0 and y = 2000m. Plot the height y and the velocity v versus time for t = 0 to t = 30 seconds.6&7-123. Develop Algorithm (processing steps to solve problem)The ODE describing the motion of a falling object is:222dtdybmgdtydmmbv2-mg2vmbgdtdvvdtdywhich is equivalent to the following system of ODEs:6&7-134. Write the “Function" (Code)(instruction sequence to be carried out by the computer)Use the Matlab editor to create a file yvprime.m .Function yvprime has two inputs:t: a scalar value for timeyv: a vector containing [y , v] values at time = tFunction yvprime has one output: yvp: a vector containing [ dy/dt ; dv/dt]function yvp = yvprime(t, yv)% function yvp = yvprime(t, yv)m = 70;g = 9.8;b = .25;y = yv(1);v = yv(2);yvp = [v ; (-g +(b/m)*v.^2)];2vmbgdtdvvdtdy6&7-146. Run CodeTo run our program we will use ‘ode45’ in a slightly different way as follows:[t, yv] = ode45(@yvprime, [0,30], [2000;0]); [t, yv] = ode45(@aux_function, time, initial_conditions);aux_function: user created function that computes the derivatestime: a vector [start stop] for the range of time of the solution.initial_condition: initial value(s) [ y(start) ; v(start)]t: solution column vector of time values from [start, stop]yv: solution matrix with two columns [ y , v] representing values y in the first column and v in the second at time t6&7-155. Test and Debug the Code To test our program in Matlab, we will plot velocity versus time >> plot(t,yv( : , 2 )) and we will plot y versus time,>> plot(t,yv( : , 1))6&7-165. Test and Debug the Code Matlab produces a discretized solution >> yv>> t6&7-17Matlab offers a family of ODE solvers that use different methods to solve ODEs:ode45, ode15i, ode15s, ode23, ode23s, ode23t, ode23tb, ode113Matlab does NOT guarantee that the solution any of the solvers produces is accurate. Knowledge that a particular solver will produce an accurate solution for a given ODE is beyond the scope of this course. CS 357 / CS 457 are courses in numerical analysis that cover methods of solving ODEs.6&7-18Most integrals arising from solutions of problems in engineering and the physical sciences cannot be represented in "closed form"- they must be evaluated numerically.• For a function of a single variable, we seek an approximation to the area "under" the curve:abSum of the “signed” areasf ( x ) dxabxf(x)+-Definite IntegralThe Matlab functions: quad (quad8) and trapz only apply to continuous functions f(x) of one variable - x ( a scalar).6&7-19f(x)x0 = a x1xn = b x2x3Xn-1. . .In this case, we approximate the integral by adding the areas of a set of approximating rectangles. yxApproximation of f(x) by using constant functions on the intervals [xk-1 , xk ] (not a good approximation).y1y2y3ynf ( x ) f ( xk) xk1 x xk, k 1, 2 ,, n6&7-20In this case, we approximate the integral by adding the areas of a set of approximating trapezoids.f(x)x =a0x1x =b n2xn-1x. . .y0y1y2Approximation of f(x) by using a linear function on the intervals
View Full Document