Title: Introduction to Numerical Solutions of Ordinary Differential Equations
1Introduction to Numerical Solutions of Ordinary
Differential Equations
- Larry Caretto
- Mechanical Engineering 501AB
- Seminar in Engineering Analysis
- October 29, 2003
2Outline
- 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
3Review 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)
4Simultaneous 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
5Matrix Differential Equations II
6Solving
- 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
7Terms in Solution of
- With these definitions, si Cielit pi/li
becomes s EC L-1p (E(0) I)
8Apply 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
9Reduction 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
10Reduction of Order II
- Continue this definition up to zn-1
- z0 to zn-1 are n variables that satisfy
simultaneous, first-order ODEs
11Reduction 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
12Simple 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
13Transforms 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.
14Solving 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)
15Partial 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
16Shifting 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
17Shifting 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)
18Numerical 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
19Interpolation
- 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
20Newton 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
21Newton 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
22Divided Difference Table
x0 y0 ?a0
?a1
x1 y1 ?a2
x2 y2 ? a3
x3 y3
23Divided Difference Example
0 0 ?a0
?a1 ?a1
10 10 ?a2
20 40 ? a3
30 100
24Divided 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)
25Constant 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
26Interpolation 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)
29Polynomial 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
30Derivative Expressions
- Obtain from differentiating interpolation
polynomials or from Taylor series - Series expansion for f(x) about x a
- What is error from truncating series?
31Truncation 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)
32Derivative 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
33Derivative 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
34Derivative Expressions III
- Apply general equation for k 1 and k 1
Forward Backward
35Derivative Expressions IV
- Subtract fi1 and fi-1 expressions
- Result called central difference expression
36Order 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
37Order 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
38Higher Order Derivatives
- Add fi1 and fi-1 expressions
- f is second-order, central difference
39Higher 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
40Specific Taylor Series
41Second Order Forward
- Subtract 4fi1 from fi2 to eliminate h2 term
Second order error
42Second Order Backwards
- Add 4fi-1 to fi-2 to eliminate h2 term
Second order error
43Other 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
44Order 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
45Roundoff 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)
47Numerical 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
48The 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
49Notation
- 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).
50More 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
51Local 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
52Eulers 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
53Euler 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
54Euler 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
55Euler Example Plot
56Error 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
57Taylor 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))
58Local 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
59Global 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
60Global 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
61Euler 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
62Better 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
63Second-order Runge Kutta
64Fourth-order Runge Kutta
- Uses four derivative evaluations per step
65Comparison 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)
68Some 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
69More 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)
72Error 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
73Runge-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
74Runge-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
75Other 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
76Systems 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
77Reduction 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
78Solving 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
79Runge-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
80ODE 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
81ODE 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
82ODE 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
83Numerical 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
84Derivative 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
85Derivative 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
86Derivative 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