MATLAB 입문 CHAPTER 8 Numerical Calculus and Differential Equations

1. MATLAB 입문CHAPTER 8 Numerical Calculus and Differential Equations ACSL, POSTECH

2. Numerical Methods for Differential Equations • The Euler Method • Consider the equation (8.5-1) r(t) : a known function • From the definition of derivative, ( is small enough) (8.5-2) • Use (8.5-2) to replace (8.5-1) by the following approximation: (8.5-3)

3. Numerical Methods for Differential Equations • Assume that the right side of (8.5-1) remains constant over the time interval . Then equation (8.5-3) can be written in more convenient form as follows: (8.5-4) , where : step size • The Euler method for the general first-order equation is (8.5-5) • The accuracy of the Euler method can be improved by using a smaller step size. However, very small step sizes require longer runtimes and can result in a large accumulated error because of round-off effects.

4. Euler method for the free response of dy/dt = –10y, y(0) = 2 r= -10; delta = 0.01 ; % tau = 0.1, step size 20% of tau y(1)=2; k=0; for time=[delta:delta:0.5] k=k+1; y(k+1) = y(k) + r*y(k)*delta; end t=[0:delta:0.5]; y_true = 2*exp(-10*t); plot(t,y,'o',t,y_true), xlabel('t'), ylabel('y');

5. Euler method solution for the free response of dy/dt= –10y, y(0) = 2. Figure 8.5–1 8-19

6. Euler method solution of dy/dt= sin t, y(0) = 0. Figure 8.5–2 8-20 More? See pages 490-492.

7. Numerical Methods for Differential Equations • The Predictor-Corrector Method • The Euler method can have a serious deficiency in problems where the variables are rapidly changing because the method assumes the variables are constant over time interval . (8.5-7) (8.5-8) • Suppose instead we use the average of the right side of (8.5-7) on the interval . (8.5-9) , where (8.5-10)

8. Numerical Methods for Differential Equations • Let Euler predictor: (8.5-11) Trapezoidal corrector: (8.5-12) • This algorithm is sometimes called the modified Euler method. However, note that any algorithm can be tried as a predictor or a corrector. • For purposes of comparison with the Runge-Kutta methods, we can express the modifided Euler method as (8.5-13) (8.5-14) (8.5-15)

9. Modified Euler solution of dy/dt = –10y, y(0) = 2 r= -10; delta = 0.02 ; % same tau as before y(1)=2; k=0; for time=[delta:delta:0.5] k=k+1; x(k+1) = y(k) + delta*r*y(k); y(k+1) = y(k) + (delta/2)*(r*y(k) + r*x(k+1)); end t=[0:delta:0.5]; y_true = 2*exp(-10*t); plot(t,y,'o',t,y_true), xlabel('t'), ylabel('y');

10. Modified Euler solution of dy/dt= –10y, y(0) = 2. Figure 8.5–3 8-21

11. Modified Euler solution of dy/dt= sin t, y(0) = 0. Figure 8.5–4 8-22 More? See pages 493-496.

12. Numerical Methods for Differential Equations • Runge-Kutta Methods • The second-order Runge-Kutta methods: (8.5-17) , where : constant weighting factors (8.5-18) (8.5-19) • To duplicate the Taylor series through the h2 term, these coefficients must satisfy the following: (8.5-19) (8.5-19) (8.5-19)

13. Numerical Methods for Differential Equations • The fourth-order Runge-Kutta methods: (8.5-23) (8.5-24) Apply Simpson’s rule for integration

14. Numerical Methods for Differential Equations • MATLAB ODE Solvers ode23 and ode45 • MATLAB provides functions, called solvers, that implement Runge-Kutta methods with variable step size. • The ode23 function uses a combination of second- and third-order Runge-Kutta methods, whereas ode45 uses a combination of fourth- and fifth-order methods. <Table 8.5-1> ODE solvers

15. Numerical Methods for Differential Equations • Solver Syntax <Table 8.5-2> Basic syntax of ODE solvers

16. Lets try these on the problem we've been solving

17. Numerical Methods for Differential Equations <example>Response of an RC Circuit (RC=0.1s, v(0)=0V, y(0)=2V)

18. Numerical Methods for Differential Equations • Effect of Step Size • The spacing used by ode23 is smaller than that used by ode45 because ode45 has less truncation error than ode23 and thus can use a larger step size. • ode23 is sometimes more useful for plotting the solution because it often gives a smoother curve. • Numerical Methods and Linear Equations • It is sometimes more convenient to use a numerical method to find the solution. • Examples of such situations are when the forcing function is a complicated function or when the order of the differential equation is higher than two. • Use of Global Parameters • The global x y z command allows all functions and files using that command to share the values of the variables x, y, and z.

20. Extension to Higher-Order Equations • The Cauchy form or the state-variable form • Consider the second-order equation (8.6-1) (8.6-2) • If , , then The Cauchy form or state-variable form

21. Extension to Higher-Order Equations <example>Solve (8.6-1) for with the initial conditions . ( Suppose that and use ode45.) %Define the function function xdot=example1(t,x) xdot(1)=x(2); xdot(2)=(1/5) * (sin(t) - 4*x(1) - 7*x(2)); xdot = [xdot(1) ; xdot(2)]; %then, use it: [t, x] = ode45('example1',[0,6],[3,9]); plot(t,x); The initial condition for the vector x. xdot(1) xdot(2) x(1) x(2)

22. Extension to Higher-Order Equations • Matrix Methods < a mass and spring with viscous surface friction > (8.6-6) ( Letting , ) (Matrix form) (8.6-7) (Compact form)

26. Extension to Higher-Order Equations • Characteristic Roots from the eig Function • Substituting , • Cancel the terms • Its roots are s = -6.7321 and s = -3.2679. (8.6-8) (8.6-8) A nonzero solution will exist for A1 and A2 if and only if the deteminant is zero

27. Extension to Higher-Order Equations • MATLAB provides the eig function to compute the characteristic roots. • Its syntax is eig(A) <example> The matrix A for the equations (8.6-8) and (8.6-9) is • To find the time constants, which are the negative reciprocals of the real parts of the roots, you type tau = -1./real (r). The time constants are 0.1485 and 0.3060.

28. Extension to Higher-Order Equations • Programming Detailed Forcing Functions < An armature-controlled dc motor> ▪ Apply Kirchhoff’s voltage low and Newton’s low (8.6-10) (8.6-11) : motor’s current , : rotational velocity ( matrix form ) Letting L: inductance, R: resistance, I: inertia, : torque constant, : back emf constant, c: viscous damping constant, : applied voltage

30. The following is extra credit (2 pts per problem)

31. ODE Solvers in the Control System Toolbox • Model Forms • The reduced form • The state-model of the same system Both model forms contain the same information. However, each form has its own advantages, depending on the purpose of the analysis. (8.7-1) (8.7-2) (8.7-3) The specification of the output: (8.7-4) (8.7-5)

32. ODE Solvers in the Control System Toolbox < Table 8.7-1> LTI object functions

33. ODE Solvers in the Control System Toolbox • ODE Solvers <Table 8.7-2> Basic syntax of the LTI ODE solvers • LTI Viewers: It provides an interactive user interface that allows you to switch between different types of response plots and between the analysis of different systems. The viewer is invoked by typing ltiview.

34. ODE Solvers in the Control System Toolbox • Predefined Input Functions < Table 8.7-3 > Predefined input functions

35. Advanced Solver Syntax • The odeset Function < Table 8.8-1 > Complete syntax of ODE solvers

36. Advanced Solver Syntax • Stiff Differential Equations • A stiff differential equation is one whose response changes rapidly over a time scale that is short compared to the time scale over which we are interested in the solution. • A small step size is needed to solve for the rapid changes, but many steps are needed to obtain the solution over the longer time interval, and thus a large error might accumulate. • The four solvers specifically designed to handle stiff equations: ode15s (a variable-order method), ode23s (a low-order method), ode23tb (another low-order method), ode23t (a trapezoidal method).