MATLAB 입문 CHAPTER 8 Numerical calculus and differential equations

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. h 1- h 1- h 1- h 1- Improper integrals and singularities • Improper integrals • singularity

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

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

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

8. 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 • 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

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

10. Integration near a Singularity Trapezoid: x=[0:0.01:1]; y=sqrt(x); A1=trapz(x,y); Simpson: A2=quad('sqrt',0,1); Lobatto: A3=quadl('sqrt',0,1); (The slope has a singularity at x=0)

11. 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. Problem 1 and 5 -- Integrate Hint: position is the integral of velocity*dt Hint: Work is the integral of force*dx

13. Problem 11 • First find V(h=4) which is the Volume of the full cup • Hint: Flow = dV/dt, so the integral of flow (dV/dt) from 0 to t = V(t) • Set up the integral using trapezoid (for part a) and Simpson (for b) • Hint: see next slide for how to make a vector of integrals with changing b values • So you can calculate vector V(t), i.e. volume over time for a given flow rate • V = quad(....., 0, t) where t is a vector of times • and then plot and determine graphically when V crosses full cup

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

15. Problem 14 Hint: The quad function can take a vector of values for the endpoint Set up a function file for current Set up a vector for the time range from 0 to .3 seconds Then find v = quad('current',0,t) will give a vector result and plot v

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

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

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

19. Problem 7: Derivative problem • Hint: velocity is the derivative of height

20. Problem 17

21. Polynomial derivatives -- trivia Example p1= 5x +2 p2=10x2+4x-3 Result: der2 = [10, 4, -3] prod = [150, 80, -7] num = [50, 40, 23] den = [25, 20, 4] • b = polyder (p) p = [a1,a2,…,an] b = [b1,b2,…,bn-1] • b = polyder (p1,p2) • [num, den] = polyder(p2,p1) Numerical differentiation functions

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

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

24. 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)

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

26. Problem 21 – 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

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

28. 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:

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

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

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