10/22/2019

Ali M. Sahlodin
Department of Chemical Engineering
AmirKabir University of Technology
1397 S.H

 Properties and challenges of DAEs

 Index reduction
 Constraint stabilization
 Consistent initialization

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 2

1
 10/22/2019

g
 Pendulum in ODE form    sin( )
L

 Use explicit-Euler method x(tk 1 )  x(tk )  hf (tk , x(tk ))

h=0.025 h=0.001

Solution eventually goes unstable, regardless of step size

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 3

 Implicit-Euler method x(tk 1 )  x(tk )  hf (tk 1, x(tk 1 ))

Simulation continues with no problem (solution does not diverge)!

What made the difference?

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 4

2
 10/22/2019

 Linear time-invariant ODE

x = Ax, A nn matrix
 Find eigenvalues of A

Solve characteristic equation A  I  0

Matrix of eigenvalues (diagonal)   diag  1, 2 ,, n 

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 5

 Find eigenvectors of A
(A  1I)v1  0 Matrix of eigenvectors
(A  2 I)v2  0
R   v1 | v2 |  | vn 

(A  n I)vn  0

Column vectors
 Eigenvalue decomposition of A

A  RR 1  x  A x  RR 1 x
 R 1 x  R 1 x

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 6

3
 10/22/2019

 Introduce state variables y

y  R 1 x  y   y
 This transformation resulted in a decoupled
system of differential equations

y j (t )   j y j (t ), j  1,..., n
 Solution of the decoupled ODEs

y j (t )  y j (t0 )exp( j t ), j  1,..., n

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 7

 Solution of the original ODE

n
x(t )  Ry(t ) x(t )   y j (t0 )v j exp( j t )
j 1

Generally, eigenvalue is a complex number   Re( )  i Im g (  )

et  e(Re( )i Img( ))t  eRe( )t ei Img(  )t

Euler formula eiz  cos( z )  i sin( z )

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 8

4
 10/22/2019

 Write the exponential term in the ODE solution as

bounded

et  eRe( )t  cos(I mg( )t )  i sin(I mg( )t ) 

 Stability of the ODE

if Re( )>0  lim eRe( )t   ODEs are unstable

t 

if Re( )  0  lim eRe( )t  0 ODEs are stable

t 

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 9

 Even if the ODE itself is stable, its numerical

solution can be unstable.

Given ODE: y   y
Stability of Forward Euler

yn1  yn  h yn  (1  h ) yn
 (1  h )( yn1  h yn1 )  yn1 (1  h ) 2
 y0 (1  h ) n1

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 10

5
 10/22/2019

For stability, we need lim yn1  0 |1  h | 1

n

Img
i i h
1  h  e  h  e  1
Re
-2 -1
Stability region of forward Euler

2
For real eigenvalue h

Step size is limited by stability: too big a step  unstable solution

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 11

yn1  yn  h yn1
Using a similar exercise
Img h

Stability region is outside the circle Re

The entire left half plane is covered.

No condition on step size for stable ODE.
Backward Euler can be stable even for unstable ODEs (how?)

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 12

6
 10/22/2019

 Why was forward Euler solution unstable?

 Rewrite as first-order ODE
 g
 g    sin 
   sin( )  L
L
  

 Linearize for stability analysis
 g 
f  0  cos   f g
 L   I  0   i cos 
x   x L
1 1 
Eigenvalues do not have a real part  Forward Euler is not
stable with any step size (see stability circle)
Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 13

 Solve using forward Euler: x  1000 x  100sin t , x(0)  1

2
For stability h   ,   1000  h  0.002


h=0.0025 h=0.0020 h=0.0015

 Explicit solvers need very small steps for stiff

ODEs  inefficient integration
 DAEs are similar to stiff ODEs
 Implicit solvers preferred

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 14

7
 10/22/2019

Numerical solution can be

stable, but very inaccurate numerical solution

true solution

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 15

Both stable solutions

h=0.001 h=0.1

Oops! Too inaccurate

The pendulum will continue to oscillate The pendulum is going to stop!

Which conclusion is correct?

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 16

8
 10/22/2019

 k1 0.04
C1  C2 C (0)  (1,1,1)
 k2 1e 4
 Stiff system (why?)
C2  C3  C1
 k3 3e7
C  C
 2 3 Implicit Euler

h=0.5 h=0.01

Both stable solutions, but completely different trends!!

Copyright © Ali M. Sahlodin, Dept. of Chemical Engineering, AmirKabir Univ. of Tech. 17