MATLAB CHAPTER 8 Numerical calculus and differential equations ACSL, POSTECH 1 Acsl, Postech MATLAB : Chapter 8. Numerical calculus and differential equations Review of integration and differentiation f(x) Engineering applications A Acceleration and velocity: b v (b) a (t )dt v (0) 0 Velocity and distance: b x(b) v (t ) x(a ) a a b A f ( x)dx a b

x Capacitor voltage and current: 1 b v(b) i (t )dt Q(a ) C a Work expended: d W f ( x)dx 0 d W kxdx 0 2 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech Integrals Properties of integrals b b b a a a

[cf ( x) dg ( x)]dx c f ( x)dx d g ( x)dx b c b a a c f ( x)dx f ( x)dx f ( x)dx Definite integrals x n 1 x b b n 1 a n 1 a x dx n 1 x a n 1 n 1 b b n 1 n 1 x b x dx ln x x a ln b ln a a 2 sin xdx cos x x 2 [cos 2 cos ] 2 x

This week's content handles definite integrals only the answers are always numeric Indefinite integrals cos xdx sin x see chapter 9 3 Acsl, Postech MATLAB : Chapter 8. Numerical calculus and differential equations Derivatives Integration g ( x) f ( x)dx differentiation f ( x) dg dx h( x ) f ( x ) / g ( x ) dh dg df f ( x) g ( x) dx dx dx df dg g ( x) f ( x) dh dx

dx 2 dx g z f ( y ) y g (x) dz dz dy dx dy dx h( x ) f ( x ) g ( x ) : product rule : quotient rule : chain rule example dx n nx n 1 1) dx d ln x 1 dx x d sin x cos x dx d cos x sin x dx 2)

d ( x 2 sin x) x 2 cos x 2 x sin x dx d (sin 2 x) dy 2 y 2 sin x cos x 3) dx dx 4 Acsl, Postech MATLAB : Chapter 8. Numerical calculus and differential equations Numerical integration Rectangular integration trapezoidal integration y y sin dx 0 y=f(x) 0 y=f(x) 0

cos 0 cos 2 (exact solution) a sin dx cos x x b a x (use of the trapz function) b Numerical integration functions Command Description quad (function,a,b,tol) Uses an adaptive Simpsons rule to compute the integral of the function function with a as the lower integration limit and b as the upper limit. The parameter tol is optional. Tol indicates the specified error tolerance. quadl (function, a,b,tol) Uses Lobatto quadrature to compute the integral of the function function. The rest of the syntax is identical to quad. trapz (x,y)

Uses trapezoidal integration to compute the integral of y with respect to x, where the array y contains the function values at the points contained in the array x. 5 Acsl, Postech MATLAB : Chapter 8. Numerical calculus and differential equations Rectangular, Trapezoidal, Simpson Rules Rectangular rule N1 I h f i f(x2) f(x1) f(x) Simplest, intuitive i 0 wi = 1 Error=O(h) Mostly not good enough! x0 =a x1 x2 Trapezoidal rule Two point formula, h=b-a Linear interpolating polynomial I h( f 0 f1 ) / 2 O(h 3 ) f0 f(x)

xN=b f1 Simpsons Rule Three point formula, h=(b-a)/2 Quadratic interpolating polynomial Lucky cancellation due to right-left symmetry x0 =a x1=b 1 4 1 I h( f 0 f1 f 2 ) O(h 5 ) 3 3 3 Significantly superior to Trapezoidal rule Problem for large intervals -> Use the extended formula 6 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech Matlab's Numerical Integration Functions trapz(x,y) trapezoidal integration method not as accurate as Simpson's method

very simple to use quad('fun',a,b,tol) When you have a function file When you have a vector of data adaptive Simpsons rule integral of fun from a to b tol is the desired error tolerance and is optional quad1('fun',a,b,tol) Alternative to quad adaptive Lobatto integration integral of fun from a to b tol is the desired error tolerance and is optional 7 Acsl, Postech MATLAB : Chapter 8. Numerical calculus and differential equations Comparing the 3 integration methods test case: sin dx cos x

cos 0 cos 2 0 0 (exact solution) Trapezoid Method: x=linspace(0,pi, 10); y=sin(x); trapz(x,y); Simpson's Method: quad('sin',0,pi); ( Lobatto's Method quadl('sin',0,pi); The answer is A1=A2=0.6667, A3=0.6665 8 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech To obtain a vector of integration results: For example, sin(x) is the integral of the cosine(z) from z=0 to z=x In matlab we can prove this by calculating the quad integral in a loop: for k=1:101 x(k)= (k-1)* pi/100; sine(k)=quad('cos',0,x(k)); end

plot(x, sine) % x vector will go from 0 to pi % this calculates the integral from 0 to x % this shows the first half of a sine wave % calculated by integrating the cosine 9 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech Using quad( ) on homemade functions 1) Create a function file: function c2 = cossq(x) % cosine squared function. c2 = cos(x.^2); Note that we must use array exponentiation. 2) The quad function is called as follows: quad('cossq',0,sqrt(2*pi)) 3) The result is 0.6119. 10 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech 11 Acsl, Postech MATLAB : Chapter 8. Numerical calculus and differential equations Numerical differentiation dy y lim dx x 0 x True slope y3

y=f(x) y2 B A y1 x x1 x x2 y 2 y1 y 2 y1 x 2 x1 x y3 y 2 y3 y 2 mB x3 x 2 x mA C x3 mC : backward difference : forward difference mA mB 1 y 2 y1 y 3 y 2 y 3 y1 : central difference 2 2 x x 2x

12 Acsl, Postech MATLAB : Chapter 8. Numerical calculus and differential equations The diff function The diff Function d= diff (x) example x=[5, 7, 12, -20]; d= diff(x) d=2 5 -32 Backward difference & central difference method example 13 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech Try that: Compare Backward vs Central Construct function with noise x=[0:pi/50:pi]; n=length(x); % y=sin(x) +/- .025 random error y=sin(x) + .05*(rand(1,51)-0.5); td=cos(x); % true derivative Backward difference using diff: d1= diff(y)./diff(x); subplot(2,1,1)

plot(x(2:n),td(2:n),x(2:n),d1,'o'); xlabel('x'); ylabel('Derivative'); axis([0 pi -2 2]) title('Backward Difference Estimate') Central Difference d2=(y(3:n)-y(1:n-2))./(x(3:n)-x(1:n-2)); subplot(2,1,2) plot(x(2:n-1),td(2:n-1),x(2:n-1),d2,'o'); xlabel('x'); ylabel('Derivative'); axis([0 pi -2 2]) title('Central Difference Estimate'); 14 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech Analytical solutions to differential equations(1/6) Solution by Direct Integration Ordinary differential equation (ODE) example dy 2 t dt t dy t3 t t3 2 0 dt dt 0 t dt 3 0 3 t y (t ) y (0) t3

3 dy dt d2y y (t ) 2 dt y (t ) Partial differential equation (PDE) -- not covered 15 Acsl, Postech MATLAB : Chapter 8. Numerical calculus and differential equations Analytical solutions to differential equations(2/6) Oscillatory forcing function dy f (t ) dt Forcing function example f (t ) sin wt t dy t dt dt sin wtdt 0 A second-order equation

0 cos wt t 1 cos wt w 0 w t 1 cos wt y (t ) y (t ) y (0) 0 w 1 cos wt y (t ) y (0) w d2y 3 t dt 2 t d2y dy t4 3 0 dt 2 dt dt y (0) 0 t dt 4 t dy t 4 y (0) dt 4 4 tt dy t5 dt y ( t ) y

( 0 ) y ( 0 ) dt ty (0) 0 dt 0 4 20 t y (t ) t 5 / 20 ty (0) y (0) 16 Acsl, Postech MATLAB : Chapter 8. Numerical calculus and differential equations Analytical solutions to differential equations(3/6) Substitution method for first-order equations dy y f (t ) ( y (t ) Ae st , dy / dt sAe st dt ) f (t ) 0 suppose f (t ) 0

t 0 for M at t=0 Such a function is the step function y (t ) Ae st B dy y sAe st Ae st 0 dt B y (0) A s 1 0 Characteristic equation y (t ) Ae st y (0) A s 1 / Characteristic root ( s 1 0 , y (0) Ae 0 A y (t ) y (0)e t / to find A A y (0) M ) y (t ) y (0)e t / M (1 e t / ) Forced response Solution, Free response

The solution y(t) decays with time if >0 is called the time constant. Natural (by Internal Energy) v.s. Forced Response (by External Energy) Transient (Dynamics-dependent) v.s. Steady-state Response (External Energydependent) 17 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech solve on paper then plot (next week we solve with Matlab) Use the method in previous slide 18 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech solve on paper then plot (next week we solve with Matlab) Re-arrange terms: mv' + cv = f , OR (m/c) v' + v = f/c now it's in the same form as 2 slides back 19 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech Analytical solutions to differential equations(4/6) Nonlinear equations

example yy 5 y y 0 y sin y 0 y y 0 Substitution method for second order equations(1/3) 1) y c 2 y 0 ( y (t ) Ae st ( s 2 c 2 ) Ae st 0 s c y (t ) A1e ct A2 e ct suppose that ) general solution c 2, y (0) 6, y (0) 4 y (0) 6 A1 A2 y (0) 2 A1 2 A2 4 solution A1 4, A2 2 20 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech Analytical solutions to differential equations(5/6) Substitution method for second order equations(2/3) 2)

y 2 y 0 ( y (t ) Ae st ) y (t ) A1e it A2 e it Eulers identities general solution e it cos t i sin t y (t ) B1 sin t B2 cos t B1 y (0) / , B2 y (0) y (t ) period 3) y (0) sin t y (0) cos t 2 frequenc f 1 / P P y my cy ky f (t ) ( f (t ) 0, y (t ) Ae st ) (ms 2 cs k ) Ae st 0 1. Real and distinct: s1 and s2. y (t ) A1e s1t A2 e s2t 2. Real and equal: s1. y (t ) ( A1 tA2 )e s1t

3. Complex conjugates: s a i at y (t ) B1e sin t B2 e at cos t 21 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech Analytical solutions to differential equations(6/6) Substitution method for second order equations(3/3) 1. real, distinct roots: m=1,c=8, k=15. characteristic roots s=-3,-5. y (t ) A1e 3t A2 e 5t 2. complex roots: m=1,c=10,k=601 characteristic rootss 5 24i y (t ) B1e 5t sin 24t B2 e 5t cos 24t P 2 / 24 / 12 3. unstable case, complex roots: m=1,c=4,k=20 characteristic roots s 2 4i y (t ) B1e 2t sin 4t B2 e 2t cos 4t P 2 / 4 / 2 4. unstable case, real roots: m=1,c=3,k=-10 characteristic roots s=2,-5. y (t ) A1e 2t A2 e 5t 22 MATLAB : Chapter 8. Numerical calculus and differential equations

Acsl, Postech Problem 22 Just find the roots to the characteristic equation and determine which form the free response will take (either 1, 2, 3 or 4 in previous slide) This alone can be a helpful check on matlab's solutionsif you don't get the right form, you can suspect there is an error somewhere in your solution 23 MATLAB : Chapter 8. Numerical calculus and differential equations Acsl, Postech Next Week: we learn how to solve ODE's with Matlab 24