6&7-2 Dynamics, 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-3 1. 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-4 m b v2 -mg From Newton’s 2nd law, the sum of the forces on the falling object equal it’s mass times acceleration: 2. (continued) )(**)(*2tvbgmdttdvm )(*)/()(2tvmbgdttdv6&7-5 2. (continued) Using calculus we can rewrite the differential equation in the form: m bv2 -mg and integration of both sides gives, dttvmbgdv )(*2bmggtbmgtv tanh)(6&7-6 2. (continued) Since we can write tanh(x) as: m bv2 -mg when x -> infinity tanh(x) -> 1 and so the falling body has a terminal velocity of: xxxxxxeeeeeex2211)tanh(bmgv )(6&7-7 3. 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 solution initial_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-8 4. 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/dt m = 70; g = 9.8; b = .25; vp = -g +(b/m)*v.^2; )(*)/()(2tvmbgdttdv6&7-9 5. 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 Code To 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-11 1. 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-12 3. Develop Algorithm (processing steps to solve problem) The ODE describing the motion of a falling object is: 222dtdybmgdtydmm bv2 -mg 2vmbgdtdvvdtdywhich is equivalent to the following system of ODEs:6&7-13 4. 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 time yv: a vector containing [y , v] values at time = t Function 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)]; 2vmbgdtdvvdtdy6&7-14 6. Run Code To 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 derivates time: 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-15 5. 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-16 5. Test and Debug the Code Matlab produces a discretized solution >> yv >> t6&7-17 Matlab offers a family of ODE solvers that use different methods to solve ODEs: ode45, ode15i, ode15s, ode23, ode23s, ode23t, ode23tb, ode113 Matlab 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-18 Most 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: a b Sum of the “signed” areas  f ( x ) dx a b  x f(x) + - Definite Integral The Matlab functions: quad (quad8) and trapz only apply to continuous functions f(x) of one variable - x ( a scalar).6&7-19 f(x) x0 = a x1 xn = b x2 x3 Xn-1 . . . In this case, we approximate the integral by adding the areas of a set of approximating rectangles. y x Approximation of f(x) by using constant functions on the intervals [xk-1 , xk ] (not a good approximation). y1 y2 y3 yn f ( x )  f ( x k ) x k  1  x  x k , k  1 , 2 , L , n6&7-20 In this case, we approximate the integral by adding the areas of a set of approximating trapezoids. f(x) x =a 0 x 1 x =b n 2 x n-1 x . . . y0 y1 y2 Approximation of f(x) by using a linear