Numerical Calculus Review and Integration Methods in MATLAB

1. MATLAB 입문CHAPTER 8 Numerical calculus and differential equations ACSL, POSTECH

2. f(x) A x a b Review of integration and differentiation • Engineering applications • Acceleration and velocity: • Velocity and distance: • Capacitor voltage and current: • Work expended:

3. Integrals • Properties of integrals • Definite integrals • Indefinite integrals This week's content handles definite integrals only the answers are always numeric see chapter 9

4. Derivatives • Integration differentiation example : product rule : quotient rule : chain rule 2) 1) 3)

5. Numerical integration Rectangular integration trapezoidal integration Numerical integration functions y y (exact solution) y=f(x) y=f(x) ··· ··· (use of the trapz function) x x ··· ··· a b a b

6. f(x2) f(x1) f(x) x1 x2 xN=b x0 =a Rectangular, Trapezoidal, Simpson Rules • Rectangular rule • Simplest, intuitive • wi = 1 • Error=O(h) Mostly not good enough! • Trapezoidal rule • Two point formula, h=b-a • Linear interpolating polynomial • Simpson’s Rule • Three point formula, h=(b-a)/2 • Quadratic interpolating polynomial • Lucky cancellation due to right-left symmetry • Significantly superior to Trapezoidal rule • Problem for large intervals -> Use the extended formula f0 f1 f(x) x1=b x0 =a

7. Matlab's Numerical Integration Functions • trapz(x,y) When you have a vector of data • 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 • adaptive Simpson’s 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

8. Comparing the 3 integration methods test case: 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); (exact solution) ( The answer is A1=A2=0.6667, A3=0.6665

9. 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; % x vector will go from 0 to pi sine(k)=quad('cos',0,x(k)); % this calculates the integral from 0 to x end plot(x, sine) % this shows the first half of a sine wave % calculated by integrating the cosine

10. 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.

12. True slope y3 y=f(x) y2 B C A y1 Δx Δx x1 x2 x3 Δx 0 Numerical differentiation : backward difference : forward difference : central difference

13. The diff function • The diff Function d= diff (x) • Backward difference & central difference method example example x=[5, 7, 12, -20]; d= diff(x) d = 2 5 -32

14. 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');

15. Analytical solutions to differential equations(1/6) • Solution by Direct Integration • Ordinary differential equation (ODE) example • Partial differential equation (PDE) -- not covered

16. Analytical solutions to differential equations(2/6) • Oscillatory forcing function • A second-order equation example Forcing function

17. suppose for M at t=0 , ( ) Such a function is the step function Characteristic equation Characteristic root , ) ( to find A Forced response Solution, Free response ∵The solution y(t) decays with time if τ>0 τ is called the time constant. Analytical solutions to differential equations(3/6) • Substitution method for first-order equations Natural (by Internal Energy) v.s. Forced Response (by External Energy) Transient (Dynamics-dependent) v.s. Steady-state Response (External Energy-dependent)

18. solve on paper then plot(next week we solve with Matlab) • Use the method in previous slide

19. 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

20. Analytical solutions to differential equations(4/6) • Nonlinear equations • Substitution method for second order equations(1/3) example 1) ( ) general solution suppose that solution

21. Analytical solutions to differential equations(5/6) • Substitution method for second order equations(2/3) 2) ( ) general solution Euler’s identities period frequency ( ) 3) 1. Real and distinct: s1 and s2. 2. Real and equal: s1. 3. Complex conjugates:

22. Analytical solutions to differential equations(6/6) • Substitution method for second order equations(3/3) • real, distinct roots: m=1,c=8, k=15. characteristic roots s=-3,-5. 2. complex roots: m=1,c=10,k=601 characteristic roots 3. unstable case, complex roots: m=1,c=-4,k=20 characteristic roots 4. unstable case, real roots: m=1,c=3,k=-10 characteristic roots s=2,-5.

23. 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 solutions—if you don't get the right form, you can suspect there is an error somewhere in your solution

24. Next Week: • we learn how to solve ODE's with Matlab