Introduction to Numerical Solutions of Ordinary Differential Equations

1
Introduction to Numerical Solutions of Ordinary
Differential Equations

• Larry Caretto
• Mechanical Engineering 501AB
• Seminar in Engineering Analysis
• October 29, 2003

2
Outline

• Review last class and homework
• Midterm Exam November 4 covers material up to and
  including last weeks lecture
• Overivew of numerical solutions
• Initial value problems in first-order equations
• Systems of first order equations and initial
  value problems in higher order equations
• Boundary value problems
• Stiff systems and eigenvalues

3
Review Last Class

• Systems of equations
• Combine to higher order equation
• Matrix approaches
• Reduction of order
• Laplace transforms for ODEs
• Transform ODE from y(t) to Y(s)
• Solve for Y(s)
• Rearrange solution
• Transform back to y(t)

4
Simultaneous Solution

• Differentiate first (y1) equation to obtain
  derivative(s) of y2 that are present in second
  equation
• Solve first equation for y2
• Substitute equations for y2 and its deriv-atives
  into the result of the first step
• Result is a higher order equation that can be
  solved (if possible) by usual methods

5
Matrix Differential Equations II

• Matrix components in

6
Solving

• Assume that A(n x n) has n linearly independent
  eigenvectors
• Eigenvectors are columns of matrix, X
• Define a new vector s X-1y (y Xs)
• Substitute y Xs into the matrix differential
  equation
• Obtain equation for s using L X-1AX that
  gives independent equations

7
Terms in Solution of

• With these definitions, si Cielit pi/li
  becomes s EC L-1p (E(0) I)

8
Apply Initial Conditions on y(0)

• Get y Xs XEC XL-1p
• At t 0, E I, and y y0 initial y
  components, giving y0 XIC XL-1p
• Premultiply by X-1 to obtain X-1y0 X-1XC
  X-1XL-1p C L-1p
• Constant vector, C X-1y0 - L-1p
• Result y XE X-1y0 - L-1p XL-1p
• Homogenous(p 0) y XEX-1y0

9
Reduction of Order

• An nth order equation can be written as a system
  of n first order equations
• We can write a general nonlinear nth order
  equation as shown below

• Redefine y as z0 and define the following
  derivatives

10
Reduction of Order II

• Continue this definition up to zn-1

• z0 to zn-1 are n variables that satisfy
  simultaneous, first-order ODEs

11
Reduction of Order III

• The main application of reduction of order is to
  numerical methods
• Numerical solutions of ODEs develop methods to
  solve a system of first order equations
• Higher order equations are solved by converting
  them to a system of first order equations

12
Simple Laplace Transforms
f(t) F(s) f(t) F(s)
ctn cn!/sn1 ceatsin wt
ctx cG(x1)/sx1 ceatsin wt
ceat c/(s a) ceatcos wt
csin wt cw/(s2 w2) ceatcos wt
ccos wt cs/(s2 w2) Additional transforms in Table 5.9, pp 297-299 of Kreyszig Additional transforms in Table 5.9, pp 297-299 of Kreyszig
csinh wt cw/(s2 - w2) Additional transforms in Table 5.9, pp 297-299 of Kreyszig Additional transforms in Table 5.9, pp 297-299 of Kreyszig
ccosh wt cs/(s2 - w2) Additional transforms in Table 5.9, pp 297-299 of Kreyszig Additional transforms in Table 5.9, pp 297-299 of Kreyszig
13
Transforms of Derivatives

• The main formulas for differential equations are
  the Laplace transforms of derivatives
• These have transform of desired function, L
  f(t) F(s) and the initial conditions, f(0),
  f(0), etc.

14
Solving Differential Equations

• Transform all terms in the differential equation
  to get an algebraic equation
• For a differential equation in y(t) we get the
  transforms Y(s) L y(t)
• Similar notation for other transformed functions
  in the equation R(s) L r(t)
• Solve the algebraic equation for Y(s)
• Obtain the inverse transform for Y(s) from tables
  to get y(t)

15
Partial Fractions

• Method to convert fraction with several factors
  in denominator into sum of individual factors (in
  denominator)
• Necessary to modify complicated solution for Y(s)
  into simpler functions whose transforms are known
• Special rules for repeated factors and for
  complex conjugates

16
Shifting Theorems

• Sometimes required for transform table
• First theorem is defined for function f(t) with
  transform F(s)
• F(s - a) has inverse of e-atf(t)
• Example Y(s) 2(s1)4/(s 1)2 1
• For F(s) s/(s2 w2) f(t) cos wt and f(t)
  cos wt for F(s) w/(s2 w2)
• Here Y(s1) 2(s1)4/(s 1)2 1 4/(s
  1)2 1 so y(t) e-t times Y(s) inverse
• y(t) e-t2 cos t 4 sin t

17
Shifting Theorems II

• Second shifting theorem
• Applies to e-asF(s), where F(s) is known
  transform of a function f(t)
• Inverse transform is f(t a) u(t a) where u(t
  a) is the unit step function
• If we have e-ass/(s2 w2), we recognize s/(s2
  w2) as F(s) for f(t) cos wt
• Thus e-ass/(s2 w2) is the Laplace transform for
  coswt(t a) u(t a)

18
Numerical Analysis Problems

• Numerical solution of algebraic equations and
  eigenvalue problems
• Solution of one or more nonlinear algebraic
  equations f(x) 0
• Linear and nonlinear optimization
• Constructing interpolating polynomials
• Numerical quadrature
• Numerical differentiation
• Numerical differential equations

19
Interpolation

• Start with N data pairs xi, yi
• Find a function (polynomial) that can be used for
  interpolation
• Basic rule the interpolation polynomial must fit
  all points exactly
• Denote the polynomial as p(x)
• The basic rule is that p(xi) yi
• Many different forms

20
Newton Polynomials

• p(x) a0 a1(x x0) a2(x x0)(x x1)
  a3(x x0)(x x1)(x x2) an-1(x x0)(x
  x1)(x x2) (x xn-2)
• Terms with factors of x xi are zero when x xi
• Use this and rule that p(xi) yi to find ai
• a0 y0, a1 (y1 y0) / (x1 x0)
• y2 a0 a1(x2 x0) a2(x2 x0)(x2 x1)
• Solve for a2 using results for a0 and a1

21
Newton Polynomials II

• y2 a0 a1(x2 x0) a2(x2 x0)(x2 x1)

• y2 a0 a1(x2 x0) a2(x2 x0)(x2 x1)
• Data determine coefficients
• Develop scheme known as divided difference table
  to compute ak

22
Divided Difference Table
x0 y0 ?a0
?a1
x1 y1 ?a2

x2 y2 ? a3

x3 y3
23
Divided Difference Example
0 0 ?a0
?a1 ?a1
10 10 ?a2

20 40 ? a3

30 100
24
Divided Difference Example II

• Divided difference table gives a0 0, a1 1, a2
  .1, and a3 1/600
• Polynomial p(x) a0 a1(x x0) a2(x x0)(x
  x1) a3(x x0)(x x1)(x x2) 0 1(x 0)
  0.1(x 0)(x 10) (1/600)(x 0)(x 10)(x
  20) x 0.1x(x 10) (1/600)x(x 10)(x 20)
• Check p(30) 30 .1(30)(20) (1/600)
  (30)(20)(10) 30 60 10 100 (correct)

25
Constant Step Size

• Divided differences work for equal or unequal
  step size in x
• If Dx h is a constant we have simpler results
• Fk Dyk/h (yk1 yk)/h
• Sk D2yk/h2 (yk2 2yk-1 yk)/h2
• Tk D3yk/h3 (yk3 3yk2 3yk1 yk)/h3
• Dnyk is called the nth forward difference
• Can also define backwards and central differences

26
Interpolation Approaches

• When we have N data points how do we interpolate
  among them?
• Order N-1 polynomial not good choice
• Use piecewise polynomials of lower order (linear
  or quadratic)
• Can match first and or higher derivatives where
  piecewise polynomials join
• Cubic splines are piecewise cubic polynomials
  that match first and second derivatives (as well
  as values)

27
(No Transcript)
28
(No Transcript)
29
Polynomial Applications

• Data interpolation
• Approximation functions in numerical quadrature
  and solution of ODEs
• Basis functions for finite element methods
• Can obtain equations for numerical
  differentiation
• Statistical curve fitting (not discussed here)
  usually used in practice

30
Derivative Expressions

• Obtain from differentiating interpolation
  polynomials or from Taylor series
• Series expansion for f(x) about x a

• Note d0f/dx0 f and 0! 1

• What is error from truncating series?

31
Truncation Error

• If we truncate series after m terms

Terms used Truncation error, em

• Can write truncation error as single term at
  unknown location (derivation based on the theorem
  of the mean)

32
Derivative Expressions

• Look at finite-difference grid with equal
  spacing h Dx so xi x0 ih
• Taylor series about x xi gives f(xi kh)
  fx0 (ik)h fik in terms of f(xi) fi

• Compact derivative notation

33
Derivative Expressions II

• Combine all definitions for compact series
  notation

• Use this formula to get expansions for various
  grid locations about x xi and use results to
  get derivative expressions

34
Derivative Expressions III

• Apply general equation for k 1 and k 1

Forward Backward
35
Derivative Expressions IV

• Subtract fi1 and fi-1 expressions

• Result called central difference expression

36
Order of the Error

• Forward and backward derivative have error term
  that is proportional to h
• Central difference error is proportional to h2
• Error proportional to hn called nth order
• Reducing step size by a factor of a reduces nth
  order error by an

37
Order of the Error Notation

• Write the error term for nth error term as O(hn)
• Big oh notation, O, denotes order
• Recognizes that factor multiplying hn may change
  slightly with h

First order forward
First order backward
Second order central
38
Higher Order Derivatives

• Add fi1 and fi-1 expressions

• f is second-order, central difference

39
Higher Order Directional

• We can get higher order derivative expressions at
  the expense of more computations
• Get second order forward and backward derivative
  expressions from fi2 and fi-2, respectively
• Combine with previous expressions for fi1 and
  fi-1 to eliminate first order error term

40
Specific Taylor Series

• General equation

• k 2

• k -2

• k 3

• k -3

41
Second Order Forward

• Subtract 4fi1 from fi2 to eliminate h2 term

Second order error
42
Second Order Backwards

• Add 4fi-1 to fi-2 to eliminate h2 term

Second order error
43
Other Derivative Expressions

• Can continue in this fashion
• Write Taylor series for fi1, fi-1, fi2, fi-2,
  fi3, fi-3, etc.
• Create linear combinations with factors that
  eliminate desired terms
• Eliminate fi term to obtain central difference
• Keep only terms in fk with k ? i for forward
  difference expressions
• Keep only terms in fk with k ? i for forward
  difference expressions
• See results on page 271 of Hoffman

44
Order of Error Examples

• Table 2-1 in notes shows error in first
  derivative for ex around x 1
• Using first- and second-order forward and
  second-order central differences
• Step h 0.4, 0.2, and 0.1
• Error ratio for doubling step size
• 4.01 to 4.02 for central differences
• 2.07 to 2.15 for first-order forward differences
• 4.32 to 4.69 for second-order forward

45
Roundoff Error

• Possible in derivative expressions from
  subtracting close differences
• Example f(x) ex f(x) ? (exh ex-h)/(2h) and
  error at x 1 is (e1h e1-h)/(2h) - e

Second order error
46
(No Transcript)
47
Numerical ODE Solutions

• The initial value problem
• Euler method as a prototype for the general
  algorithm
• Local and global errors
• More accurate methods
• Step-size control for error control
• Applications to systems of equations
• Higher-order equations

48
The Initial Value Problem

• dy/dx f(x,y) (known) with y(x0) y0
• Basic numerical approach
• Use a finite difference grid xi1 xi h
• Replace derivative by finite-difference
  approximation dy/dx ? (yi1 yi) / (xi1 xi)
  (yi1 yi) / h
• Derive a formula to compute favg the average
  value of f(x,y) between xi and xi1
• Replace dy/dx f(x,y) by (yi1 yi) / h favg
• Repeatedly compute yi1 yi h favg

49
Notation

• xi is the value of the independent variable at
  point i on the grid
• Determined from the user-selected value of step
  size (or the series of hi values)
• Can always specify exactly the independent
  variables value, xi
• yi is the value of the numerical solution at the
  point where x xi
• fi is derivative value found from xi and the
  numerical value, yi. I.e., fi f(xi, yi).

50
More Notation

• y(xi) is the exact value of y at x xi
• Usually not known but notation is used in error
  analysis of algorithms
• f(xi,y(xi)) is the exact value of the derivative
  at x xi
• e1 y(x1) y1 is the local truncation error
• This is error for one step of algorithm starting
  from known initial condition

51
Local versus Global Error

• At the initial condition we know the solution, y,
  exactly
• First step introduces some error
• Remaining steps have single step error plus
  previous accumulated error
• Ej y(xj) yj is global truncation error
• Difference between numerical and exact solution
  after several steps
• This is the error we want to control

52
Eulers Method

• Simplest algorithm, example used for error
  analysis, not for practical use
• Define favg f(xi, yi)
• Eulers method algorithm is yi1 yi hifi yi
  hi f(xi, yi)
• Example dy/dx x y, y 0 at x 0
• Choose h 0.1
• We have x0 0, y0 0, f0 x0 y0 0 x1 x0
  h 0.1, y1 y0 hf0 0

53
Euler Example Continued

• Next step is from x1 0.1 to x2 0.2
• f1 x1 y1 0.1 0 0.1
• y2 y1 hf1 0 (.1)(.1) .01
• For dy/dx x y, we know the exact solution is
  (x0 y0 1)ex-x0 x 1
• For x0 y0 0, y ex x 1
• Look at application of Euler algorithm for a few
  steps and compute the error

54
Euler Example
xi yi fi y(xi) E(xi)
0 0 0 0 0
.1 0 .1 .0052 .0052
.2 .01 .21 .0214 .0114
.3 .031 .331 .04986 .01886
.4 .0641 .4641 .091825 .027725
55
Euler Example Plot
56
Error Propagation

• Behavior of Euler algorithm is typical of all
  algorithms for numerical solutions
• Error grows at each step
• We usually do not know this global error, but we
  would like to control it
• Look at local error for Euler algorithm
• Then discuss general relationship between local
  and global error

57
Taylor Series to Get Error

• Expand y(x) in Taylor series about x a

• Look at one step from known initial condition, a
  x0, to x0 h so x a h

• In ODE notation, dy/dx0 f(x0, y(x0))

58
Local Euler Error

• Result of Taylor series on last chart

Euler Algorithm Truncation Error

• This is only the Euler algorithm for the first
  step when we know f(x0,y(x0))
• This gives the local truncation error
• Local truncation error for Euler algorithm is
  second order

59
Global Error

• We will show that a local error of order n, has a
  global error of order n-1
• To show this consider the global error at x x0
  kh after k algorithm steps
• Is approximately k times the local error
• If local error is O(hn), approximate global error
  after k steps is k O(hn) ? kAhn
• A new step size, h/r, takes kr steps to get to
  the same x value

60
Global Error Concluded

• Compare error for same x kh with step sizes h
  and h/r
• Exkh(h) ? kAhn
• Exkh(h/r) krA(h/r)n

• When we reduce the step size by a factor of 1/r
  we reduce the error by a factor of 1/rn-1 this
  is the behavior of an algorithm whose error is
  order n-1

61
Euler Local and Global Error

• Previously showed Euler algorithm to have second
  order local error
• Should have first order global error
• Results for previous Euler example at x 1 with
  different step sizes

62
Better Algorithms

• Seek high accuracy with low computational work
• Could improve Euler accuracy by cutting step
  size, but this is not efficient
• Use other algorithms that have higher order
  errors
• Runge-Kutta methods typically used
• This is a class of methods that use several
  function evaluation methods per step

63
Second-order Runge Kutta

• Huens method

• Modified Euler method

64
Fourth-order Runge Kutta

• Uses four derivative evaluations per step

65
Comparison of Methods

• Look at Euler, Heun, Modified Euler and
  fourth-order Runge-Kutta
• Solve dy/dx e-y-x with y(0) 1
• Compare numerical values to exact solution y
  ln( ey0 e-x0 e-x)
• Look at errors in the methods at x 1 as a
  function of step size
• Compare error propagation (increase in error as x
  increases)

66
(No Transcript)
67
(No Transcript)
68
Some Basic Concepts

• A numerical method is convergent with the
  solution of the ODE if the numerical solution
  approaches the actual solution as h ? 0 (with
  increase in numerical precision at smaller h)
• Mainly a theoretical concept algorithms in use
  are convergent
• Stability refers to the ability of a numerical
  algorithm to damp any errors introduced during
  the solution

69
More on Stability

• Finite difference equations in numerical
  algorithms, when iterated, may be numerically
  increase without bound
• Stability usually is obtained by keeping step
  size h small, sometimes smaller than the h
  required for accuracy
• For most ODEs stability is not a problem, but it
  is for stiff systems of ODEs and for partial
  differential equations

70
(No Transcript)
71
(No Transcript)
72
Error Control

• How do we choose h?
• Want to obtain a result with some desired small
  global error
• Can just repeat calculations with smaller h until
  two results are sufficiently close
• Can use algorithms that estimate error and adjust
  step size during the calculation based on the
  error

73
Runge-Kutta-Fehlberg

• Uses two equations to compute yn1, one has
  O(h5), the other O(h6) error
• Requires six derivative evaluations per step
  (same evaluations used for both equations)
• The error estimate can be used for step size
  control based on an overall 5th order error
• Cask-Karp version different coefficients

74
Runge-Kutta-Fehlberg II

• See page 377 in Hoffman for algorithm
• Typical formula components below
• yn1 yn (16k1/135 6656k2/12825
• k3 hf(xn 3h/8, yn 3k1/32 9k2/32)
• Error k1/360 - 128k3/4275
• hnew holdEDesired/Error1/5
• EDesired is set by user
• RKF45 code by Watts and Shampine

75
Other Error Controls

• Do fourth-order Runge-Kutta two times, with step
  sizes different by a factor of two and compare
  results
• Runge-Kutta-Verner is similar to
  Runge-Kutta-Fehlberg, but uses higher order
  expressions

76
Systems of Equations

• Any problem with one or more ODEs of any order
  can be reduced to a system for first order ODEs
• Last week we reduced a system of two second order
  equations to a system of four first order
  equations
• In this process the sum of the orders is constant
• Example from last week

77
Reduction of Order Example

• Define variables y3 dy1/dt and y4 dy2/dt
• Then dy3/dt d2y1/dt2 and dy4/dt d2y2/dt2
• Have four simultaneous first-order ODEs

78
Solving Simultaneous ODEs

• Apply same algorithms used for single ODEs
• Must apply each step and substep to all equations
  in system
• Key is consistent determination of fi(x,y)
• All yi values in y must be available at the same
  x point when computing the fi
• E.g., in Runge-Kutta we must evaluate k1 for all
  equations before finding k2

79
Runge-Kutta for ODE System

• y(i) is vector of dependent variables at x
  xi
• k(1), k(2), k(3), and k(4), are vectors
  containing intermediate Runge-Kutta results
• f is a vector containing the derivatives
• k(1) hf(xi, y(i))
• k(2) hf(xi h/2, y(i) k(1)/2)
• k(3) hf(xi h/2, y(i) k(2)/2)
• k(4) hf(xi h, y(i) k(3))
• y(i1) (k(1) 2k(2) 2k(3) k(4))/6

80
ODE System by RK4

• dy/dx -y z and dz/dx y z with y(0) 1
  and z(0) -1 with h .1
• k(1)y h-y z 0.1-1 (-1) -.2
• k(1)z hy - z 0.11 - (-1) .2
• k(2)y h-(y k(1)y/2) z k(1)z/2 0.1
  -(1 -0.2/2) (-1 .2/2) -.18
• k(2)z h(y k(1)y/2) (z k(1)z/2) 0.1(1
  -0.2)/2 - (-1 .2/2) .18

81
ODE System by RK4 II

• k(3)y h-(y k(2)y/2) z k(2)z/2 0.1
  -(1 -0.18/2) (-1 .18/2) -.182
• k(3)z h(y k(2)y/2) (z k(2)z/2) 0.1(1
  -0.18)/2 - (-1 .18/2) .182
• k(4)y h-(y k(3)y) z k(3)z 0.1 -(1
  -0.182) (-1 .182) -.1636
• k(4)z h(y k(3)y) (z k(3)z) 0.1 (1
  -0.182) - (-1 .182) .1636

82
ODE System by RK4 III

• yi1 yi (k(1)y 2k(2)y 2k(3)y k(4)y )/6
  1 (.2) 2(.18) 2(.182)
  (.1636)/6 .8187
• zi1 zi (k(1)z 2k(2)z 2k(3)z k(4)z )/6
  1 (.2) 2(.18) 2(.182) (.1636)/6
  .8187
• Continue in this fashion until desired final x
  value is reached
• No x dependence for f in this example

83
Numerical Software for ODEs

• Usually written to solve a system of N equations,
  but will work for N 1
• User has to code a subroutine or function to
  compute the f array
• Input variables are x and y f is output
• Some codes have one dimensional parameter array
  to pass additional information from main program
  into the function that computes derivatives

84
Derivative Subroutine Example

• Visual Basic code for system of ODEs at right is
  shown below

Sub fsub( x as Double, y() as Double, _

f() as Double
) f(1) -y(1) Sqr(y(2)) y(3)Exp(2x)
f(2) -2 y(1)2 f(3) -3 y(1) y(2) End
Sub
85
Derivative Subroutine Example

• Fortran code for system of ODEs at right is shown
  below

subroutine fsub( x, y, f ) real(KIND8) x,
y(), f() f(1) -y(1) sqrt(y(2))
y(3)exp(2x) f(2) -2 y(1)2 f(3) -3
y(1) y(2) end subroutine fsub
86
Derivative Subroutine Example

• C code for system of ODEs at right is shown
  below

void fsub(double x, double y, double f)
f1 -y1 sqrt(y2) y3exp(2x)
f2 -2 y1 y1 f3 -3 y1
y2