Transcript Stability of ODEs Numerical Methods for PDEs Spring 2007 Jim E. Jones

Stability of ODEs
Numerical Methods for PDEs
Spring 2007
Jim E. Jones
References:
•Numerical Analysis, Burden & Faires
•Scientific Computing: An Introductory Survey, Heath
Stability of the ODE
The Continuous Problem
Ordinary Differential Equation: Initial Value
Problem (IVP)
y(t )  f (t , y ), a  t  b,
y (a)  
Last lecture we saw that f satisfying a Lipschitz
condition was enough to guarantee that the IVP is
well-posed.
Lipschitz condition definition
A function f(t,y) satisfies a Lipschitz condition in the
variable y on a set D in R2 if a constant L > 0 exists with
| f (t , y1 )  f (t , y0 ) | L | y1  y0 |
whenever (t,y0) and (t,y1) are in D. The constant L is called
a Lipschitz constant for f.
Well posed IVP definition
The IVP
y(t )  f (t , y ), a  t  b,
y (a)  
is a well-posed problem if:
• A unique solution y(t) exists, and
• There exists constants e0 >0 and k > 0 such that for any e in (0,e0),
whenever d(t) is continuous with |d(t)| < e for all t in [a,b], and when
|d0| < e, the IVP
z (t )  f (t , z )  d (t ), a  t  b,
z (a)    d 0
has a unique solution z(t) satisfying
| z (t )  y (t ) | ke , t  [a, b]
Difference between solutions may still be large
The IVP
y(t )  f (t , y ), a  t  b,
y (a)  
The perturbed IVP
z (t )  f (t , z )  d (t ), a  t  b,
z (a)    d 0
One can show that the difference in solutions is bounded
L (t a )
e
1
L (t a )
| z (t )  y (t ) | e
| d0 | 
max t[ a ,b ] | d (t ) |
L
Difference between solutions may still be large
The IVP
y(t )  f (t , y ), a  t  b,
y (a)  
The perturbed IVP
z (t )  f (t , z )  d (t ), a  t  b,
z (a)    d 0
One can show that the difference in solutions is bounded
L (t a )
e
1
L (t a )
| z (t )  y (t ) | e
| d0 | 
max t[ a ,b ] | d (t ) |
L
The k, and thus the difference between
solutions, may be large if L is and/or b >> a
Example
The IVP
The perturbed IVP
y(t )  100 y,0  t  10,
y (0)  1
z (t )  100 z ,0  t  10,
z (0)  1.000001
•Is the IVP well-posed?
•What’s the difference between solutions, z(t)-y(t)?
Example
The IVP
The perturbed IVP
y(t )  100 y,0  t  10,
y (0)  1
z (t )  100 z ,0  t  10,
z (0)  1.000001
•Is the IVP well-posed?
•What’s the difference between solutions, z(t)-y(t)?
To characterize the time growth (or decay) of initial
perturbations, we need the concept of stability.
Stability definition
The IVP
y(t )  f (t , y ), a  t  ,
y (a)  
is stable if:
• A unique solution y(t) exists, and
• For every e >0 there exists a d > 0 such that whenever 0 < d0 < d, the
IVP
z (t )  f (t , z ), a  t  ,
z (a)    d 0
has a unique solution z(t) satisfying
| z (t )  y (t ) | e , t  [a, )
Stability definition
The IVP
y(t )  f (t , y ), a  t  ,
y (a)  
is stable if:
• A unique solution y(t) exists, and
• For every e >0 there exists a d > 0 such that whenever 0 < d0 < d, the
IVP
z (t )  f (t , z ), a  t  ,
z (a)    d 0
has a unique solution z(t) satisfying
| z (t )  y (t ) | e , t  [a, )
For a stable ODE the difference between the
solutions is bounded for all time.
Absolute Stability definition
The IVP
y(t )  f (t , y ), a  t  ,
y (a)  
is absolutely stable if:
• A unique solution y(t) exists, and
• For every d0 the IVP
z (t )  f (t , z ), a  t  ,
z (a)    d 0
has a unique solution z(t) satisfying
lim t  | z (t )  y(t ) | 0
For an absolutely stable ODE the difference
between the solutions goes to zero as t increases.
Example
The IVP
y(t )  y, a  t  ,
y (a)  
•If  is real, what can we say about the stability,
absolute stability?
•If  is complex, what can we say about the stability,
absolute stability?
Example
•If  is real, what can we say about the stability,
absolute stability?
>0 unstable
<0 absolutely stable
•If  is complex, what can we say about the stability,
absolute stability?
Re()>0 unstable
Re()<0 absolutely stable
Re()=0 oscillating solution, stable
Systems of Ordinary Differential Equation: Initial
Value Problem (IVP)
u(t )  3u (t )  w(t )
, a  t  b,
w(t )  u (t )  3w(t )
u (a)  5
w(a)  7
Let y=(u,w)t
 

y(t )  f (t , y ), a  t  b,


y (a)  
and in this case the rhs can be described by a matrix,
3 1 


y(t )  Ay (t )  
y (t ), a  t  b,

1 3


y (a)  
Stability of Systems of Ordinary Differential
Equation
If we assume A is diagonalizable, so that its
eigenvectors are independent, then the IVP



y (t )  Ay (t ), a  t  b,


y (a)  
has solution components corresponding to each
eigenvalue i
Stability of Systems of Ordinary Differential
Equation
If we assume A is diagonalizable, so that its
eigenvectors are independent, then the IVP



y (t )  Ay (t ), a  t  b,


y (a)  
has solution components corresponding to each
eigenvalue i
Re(i) > 0 components grow exponentially
Re(i) < 0 components decay exponentially
Re(i)=0 oscillatory components
Stability of Systems of Ordinary Differential
Equation
If we assume A is diagonalizable, so that its
eigenvectors are independent, then the IVP


y(t )  Ay (t ), a  t  b,


y (a)  
has solution components corresponding to each
eigenvalue i
Unstable if Re(i) > 0 for any eigenvalue
Asymptotically stable if Re(i) < 0 for all eigenvalues
Stable if Re(i)  0 for all eigenvalues
Stability of the numerical method
The Discrete Problem
Ordinary Differential Equation: Initial Value
Problem (IVP)
y(t )  f (t , y ), a  t  b,
y (a)  
• Numerical Solution: rather than finding an analytic
solution for y(t), we look for an approximate discrete
solution wi (i=0,1,…,n) with wi approximating y(a+ih)
t=a
h
t=b
t
Local truncation error definition
The difference method
w0  
wi 1  wi  h(c f (ti 1 , wi 1 )  c f (ti , wi ))
has local truncation error
 i 1 
y (ti 1 )  y (ti )
 [c f (ti 1 , y (ti 1 ))  c f (ti , y (ti ))]
h
for each i=0,1,…,n-1
The local truncation error is a measure of the degree
to which the true IVP solution fails to satisfy the
difference equation.
Local truncation error definition
The difference method
w0  
wi 1  wi  h(c f (ti 1 , wi 1 )  c f (ti , wi ))
has local truncation error
 i 1 
y (ti 1 )  y (ti )
 [c f (ti 1 , y (ti 1 ))  c f (ti , y (ti ))]
h
for each i=0,1,…,n-1
The local truncation error is error in single step,
assuming the previous step is exact, scaled by the
mesh size h.
Definition of consistency and convergence
A difference method is consistent with the differential equation
if
lim h0{max 0i n |  i (h) |}  0.
A difference method is convergent (or accurate) with respect to
the differential equation if
lim h0{max 0i n | wi  y(ti ) |}  0.
Definition of consistency and convergence
A difference method is consistent with the differential equation
if
lim h0{max 0i n |  i (h) |}  0.
This is the error in a single
step: the local error
A difference method is convergent (or accurate) with respect to
the differential equation if
lim h0{max 0i n | wi  y(ti ) |}  0.
This is the total error : the
global error
See: http://www.cse.uiuc.edu/iem/ode/eulrmthd/
Numerical Stability definition
Apply the numerical method to
y(t )  f (t , y ), a  t  ,
y (a)  
to generate discrete solution w and apply the same method to the perturbed
Problem to generate solution u
z (t )  f (t , z ), a  t  ,
z (a)    d
The method is stable if for every e there exists a K such that
di | ui  wi | K , i
whenever d < e.
Stability of Euler’s method
Forward Euler
wi  wi 1  hf (ti 1 , wi 1 )
Apply each method to the IVP
y(t )  y, a  t  ,
y (a)  
Backward Euler
wi  wi 1  hf (ti , wi )
Solutions generated by Euler’s method
Forward Euler
wi  (1  h )i 
Backward Euler
i
 1 
wi  

 1  h 
The quantity inside the parens is called the growth
factor r. For the difference between the solution and
the perturbed solution, we have
d i | r |i d
and the requirement for stability is that |r|  1
Stability analysis of forward Euler’s method
Forward Euler
wi  (1  h )i 
For complex , stability requires that h must lie inside the
circle with radius 1 in the complex plane centered at -1.
If we consider real  for which the ODE is stable,  < 0,
the stability requirement is
h  2 / 
Stability analysis of backward Euler’s method
Backward Euler
i
 1 
wi  

 1  h 
If we consider  for which the ODE is stable, Re() < 0,
the stability of backward Euler is assured: the growth factor
is less than 1 in magnitude. Backward Euler is
unconditionally stable.
Examples from last time
• The example (2.a and 2.b) from last time illustrate the
difference in stability of forward and backward Euler.
Ordinary Differential Equation: Example 2.a
y(t )  2 y,0  t  100,
y ( 0)  1
Forward Euler with h=0.1
wi 1  wi  2hwi  (1  2h) wi
w1 
4
w0
5
Backward Euler with h=0.1
wi 1  wi  2hwi 1  wi 1  (1  2h) 1 wi
w1 
5
w0
6
Both computed solutions go to zero as t increases like
the true ODE solution
y  e 2t
Ordinary Differential Equation: Example 2.b
y(t )  2 y,0  t  100,
y ( 0)  1
Forward Euler with h=1.1
wi 1  wi  2hwi  (1  2h) wi
w1 
6
w0
5
Backward Euler with h=1.1
wi 1  wi  2hwi 1  wi 1  (1  2h) 1 wi
w1 
5
w0
16
Backward Euler go to zero as t increases. Forward
Euler blows up.
Examples from last time
• The example (2.a and 2.b) from last time illustrate the
difference in stability of forward and backward Euler.
• Another example at
http://www.cse.uiuc.edu/iem/ode/stiff/
– Here the step size for stability (h=0.02) is tighter than
one needs to control truncation error if one is not
interested in resolving the fast decaying initial transient
component of the solution.
Numerical Stability definition
A numerical method may be unstable, using the previous definition, because
the underlying ODE itself is unstable. To focus specifically on the numerical
method, we can alternatively define stability as:
A method is stable if the numerical solution at any arbitrary but fixed time t
remains bounded as h goes to zero.