MATLAB CHAPTER 8 Numerical Calculus and Differential Equations

1
MATLAB ??CHAPTER 8 Numerical Calculus
and Differential Equations

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
• k0
• for timedeltadelta0.5
• kk1
• y(k1) y(k) ry(k)delta
• end
• t0delta0.5
• y_true 2exp(-10t)
• plot(t,y,'o',t,y_true), xlabel('t'), ylabel('y')

5
8-19
Euler method solution for the free response of
dy/dt 10y, y(0) 2. Figure 8.51
6
8-20
Euler method solution of dy/dt sin t, y(0) 0.
Figure 8.52
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
• k0
• for timedeltadelta0.5
• kk1
• x(k1) y(k) deltary(k)
• y(k1) y(k) (delta/2)(ry(k) rx(k1))
• end
• t0delta0.5
• y_true 2exp(-10t)
• plot(t,y,'o',t,y_true), xlabel('t'), ylabel('y')

10
8-21
Modified Euler solution of dy/dt 10y, y(0)
2. Figure 8.53
11
8-22
Modified Euler solution of dy/dt sin t, y(0)
0. Figure 8.54
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 Simpsons 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.
• ltTable 8.5-1gt ODE solvers

15
Numerical Methods for Differential Equations

• Solver Syntax
• ltTable 8.5-2gt Basic syntax of ODE solvers

16
Lets try these on the problem we've been solving
17
Numerical Methods for Differential Equations

• ltexamplegt Response of an RC Circuit
• (RC0.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.

19
(No Transcript)
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

• ltexamplegt Solve (8.6-1) for with
  the initial conditions .
• ( Suppose that
  and use ode45.)
• Define the function
• function xdotexample1(t,x)
• xdot(1)x(2)
• xdot(2)(1/5) (sin(t) - 4x(1) - 7x(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
• lt a mass and spring with viscous surface
  friction gt
• 
  

(8.6-6)
( Letting , )
(Matrix form)
(8.6-7)
(Compact form)
23
(No Transcript)
24
(No Transcript)
25
(No Transcript)
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)
• ltexamplegt 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
• lt An armature-controlled dc motorgt

? Apply Kirchhoffs voltage low and Newtons low
(8.6-10)
(8.6-11)
motors current , rotational velocity
( matrix form )
Letting
L inductance, R resistance, I inertia,
torque constant, back emf constant, c
viscous damping constant, applied voltage
29
(No Transcript)
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

• lt Table 8.7-1gt LTI object functions

33
ODE Solvers in the Control System Toolbox

• ODE Solvers
• ltTable 8.7-2gt 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
• lt Table 8.7-3 gt Predefined input functions

35
Advanced Solver Syntax

• The odeset Function
• lt Table 8.8-1 gt 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).

37
(No Transcript)
38
(No Transcript)