OR682Math685CSI700

1
OR682/Math685/CSI700

• Lecture 8
• Fall 2000

2
Ordinary Differential Equations

• Given an equation involving the derivatives of a
  function, solve for the function
• True solution is a function we compute a finite
  approximation
• Ordinary all derivatives are with respect to a
  single variable, often time

3
Boundary Conditions

• Solve y? f (t,y) where t is a real variable,
  and f and y are vectors of length n
• Solution is not unique
• Boundary conditions must be specified
• Initial value problem y (t0) y0
• For example (1-dimensional)
• y? y, y (t0) y0 y(t) y0et

Matlab functions ode15s Matlab m-files pop.m,
pop1.m
4
Higher-Order Problems

• Equations with higher derivatives can be
  transformed to first-order problems
• If u(n) f (t,u,u?,,u(n1)) we define y1 u,
  y2 u?, , yn(t) u(n 1)), so that

5
Example 2nd-order Equation

• Falling object with air resistance (k constant)
• Here g gravity, u height, v velocity
• Transform to (y1 u(t), y2 u?(t))

Matlab m-files para.m, para1.m, pararun.m
6
Stability of ODEs

• Unstable members of the solution family move
  apart with time errors are magnified
• Stable members of the solution family move
  together with time errors reduced
• Neutrally stable solution curves neither
  converge nor diverge
• All forms of behavior can occur in a single
  problem as time varies

7
Measuring Stability

• For the system y? f (t, y), stability is
  measured in terms of the Jacobian J fy
• Look at the maximum of the real parts of the
  eigenvalues of J
• Positive implies unstable
• Negative implies stable
• Zero implies neutrally stable
• The Jacobian J changes with time

Matlab m-files stab_ex.m
8
Stability Linear System

• Consider y? Ay, y(0) y0, for a matrix A
• Let A have eigenvalues ?i with eigenvectors ui,
  and write y0 in the form
• Then

9
Linear System (continued)

• The Jacobian is J A
• Eigenvalues of A with
• positive real part yield exponentially growing
  components (unstable)
• negative real part yield exponentially decaying
  components (stable)
• Zero real part yield oscillatory components
  (neutrally stable)

Matlab m-files prey.m, prey1.m, runprey.m,
runprey1.m, ploteig.m, ploteig1.m
10
Numerical Solution of ODEs

• Analytical solution formula (Matlab dsolve)
• Numeric solution table of values
• Starting at t0, the algorithm tracks the
  trajectory defined by the ODE
• The ODE defines the slope of the trajectory
• Use slope to predict y(t1)
• Errors cause the estimates to lie on different
  curves in the solution family
• If ODE unstable, errors are magnified as t
  increase

Matlab function dsolve Matlab m-file pop2.m
11
Eulers Method

• Simple, but representative algorithm
• Based on Taylor series approximation
• Eulers method drop high-order terms

Matlab m-file eul.m
12
Types of Errors

• Two distinct sources of error
• Rounding due to finite precision of arithmetic
• Truncation due to method used (formula is not
  exact, high-order terms are ignored)
• Typically, truncation error is much larger than
  rounding error (many orders of magnitude larger)
• We will ignore rounding error

13
Local Global Error

• Truncation error comes in two forms
• Local error made in one step of a numerical
  method (error in yk if yk1 were exact) Lk
• Global total error in yk (local error plus
  effect of cumulative error) Ek yk y(tk)
• Unstable equation global error will usually
  dominate
• Stable equation local error will usually
  dominate
• We can only control local error

14
Order of Numerical Method

• A numerical method is said to have accuracy of
  order p if Lk O(h p1)
• Reason a rough approximation of the global error
  is given by the sum of the local errors is
  proportional to h p1/ h h p

15
Analysis Eulers Method

• Formulas for Taylor series Eulers method
• Global error (consider 1-dimensional case)
• Eulers method is a first-order method (letyk
  y(tk))

16
Growth Factor

• Use the mean-value theorem to show that
• for some unknown ?
• Then (from our analysis above)
• The global error is multiplied by the growth
  factor (1hk J)

17
Stability of Eulers Method

• A numerical method is stable if small errors do
  not cause numerical solutions to diverge
• Can be due to ODE itself
• Eulers method is stable if the growth factor
  satisfies 1 h J lt 1
• Equivalent to 2 lt h J lt 0

18
Stability (continued)

• If equation is stable (J lt 0) then Eulers method
  is stable if 0 lt h lt 2 / J
• For systems Eulers method is stable if the
  maximum real part of the eigenvalues of (I h
  J) is less than 1.

Matlab m-files runpreye.m, ploteige.m
19
Step-size Control

• Have algorithm automatically choose h
• Make h as large as possible
• Maintain stability and accuracy
• Based on estimates of local error, and stability
  analysis of algorithm
• For Eulers method, local error is h2 y??/2, so

20
For Next Class

• Homework see web site
• Reading
• Heath chapter 9 (rest of chapter)
• Heath chapter 10