Transcript MA375 - Rice U - Computational and Applied Mathematics

MA557/MA578/CS557
Lecture 5a
Spring 2003
Prof. Tim Warburton
[email protected]
1
Homework 2 correction
Q4) Implement the numerical approximation of:


3
0
t
x
 ( x, t  0)  e
 x2
Geometry:
By:
 N i 
 i 1 
x1  4, xN  4, xi  
x

 1 
 xN
 N 1 
 N 1 
Scheme:
in1  1    in     in1
dt
 3
dx
Initial Condition :
  xi  xi 1 

   
,t  0 
 2 

Boundary Condition :
0
i
 0
0
0
2
Homework 2 cont
Q4 cont)
For N=10,40,160,320,640,1280 run to t=10, with:
dx = (8/(N-1))
dt = dx/6
Use

n 1
0
dt  n  dt  n

 1  u   0   u   N 1
dx 

 dx 
On the same graph, plot t on the horizontal axis and error on the
vertical axis. The graph should consist of a sequence of 6
curves – one for each choice of dx.
Comment on the curves.
NOTE: For the purposes of this test we define error as:

x x

error n  max  in    i i 1 , ndt  
1i  N
 2


3
Lecture 5
• We will define stability for a numerical
scheme and investigate stability for the
upwind scheme.
• We will compare this scheme with a finite
difference scheme.
• We will consider alternative ways to
approximate the flux functions.
4
Recall
Basic Upwind Finite Volume Method
d
 qi dx   uq  xi1 , t   uq  xi , t 
dt
Approximate fluxes with upwind flux
d
 qi dx   uqi  uqi1
dt

n1
i
 in  dx
dt
Approximate time derivative and look for
solution in  qin
 u in  u in1
 u
dt
dx
in1  1    in     in1
Note we must supply a value for the left most average at each time step:

n
0 5
Convergence
• We have constructed a physically reasonable numerical scheme to
approximate the advection equation.
• However, we need to do some extra analysis to determine how
good at approximating the true PDE the discrete scheme is.
• Let us suppose that the i’th subinterval cell average of the actual
solution to the PDE at time T=n*dt is denoted by
1
qin 
xi 1  xi
xi 1

xi
x
1 i1
q  x, ndt  
q  x, ndt 

dx xi
where q satisfies:
d
dt
xi 1

xi
q  x, t dx 
d
 dxqi   uq  xi1 , t   uq  xi , t 
dt
6
Error Equation
• The goal is to estimate the difference of the exact
solution and the numerically obtained solution at
some time T=n*dt.
• So we are interested in the error: E  q  
n
i
n
i
n
i
T
,n
dt
• For the given finite volume scheme dt and dx will
be related in a fixed manner (i.e. dt = Cdx for
some C, independent of dx).
• Suppose we let dt  0 and
E n  O  dt s 
the scheme is said to be of order s.
then
7
Norms and Definitions
• We define the discrete p-norms:
E
p
1/ p

p
  dx  Ei 
 i

i 
• We say that the scheme is convergent at time
T in the norm ||.|| if: lim E n  0
dtndt0T 
• It is said to be accurate of order s if:
E n  O  dt s  as dt  0
8
Local Truncation Error
• Suppose at the beginning of a time step we actually have the exact
solution -- one question we can ask is how large is the error
committed in the evaluation of the approximate solution at the end
of the time step.
• i.e. choose
• Then
in  qin
 dt  n  dt  n
  1  u  qi   u  qi 1
 dx 
 dx 
1 n1
n
Ri :  i  qin1 
dt
n 1
i
• We expand qin1 , qin1about xi,tn with Taylor series:
 ____
2 

q
dx

q
 
n
n
3

qi 1  qi  dx  

O
(
dx
)
2



x
2

x
 




____
2
 ____
2 
 q  dt   q 
qin 1  qin  dt  
 O (dt 3 )
2
 t  2  t 
____
2
9
Estimating Truncation Error
• Inserting the formulas for the expanded q’s:
Rin :
1   dt  n  dt  n
n 1 
1

u
q

u
q

q

 i 
 i 1 i 
dt   dx 
 dx 





 1  dt u  q n

i


  dx 



____
____



2
2


1
 dt  
 q  dx  q
3 
    u  qin  dx  

O
(
dx
)

2 
dt   dx  
 x  2  x 




 

____
____

2
2
  n

 q  dt   q 
3 
 O (dt ) 
  qi  dt  


2 
 t  2  t 
 



 

10
Estimating Truncation Error
• Removing canceling terms:




 1  dt u  q n

i


  dx 



____
____



2
2


1
 dt  
 q  dx  q
3 
Rin :    u  qin  dx  

O
(
dx
)

2 
dt   dx  
 x  2  x 




 

____
____

2
2
  n

 q  dt   q 
3 
 O (dt ) 
  qi  dt  


2 
 t  2  t 
 



 

11
Estimating Truncation Error
• Simplifying:
____
____



2
2
 q  dx   q 
  dt  
3 
   dx u  dx  x  2  x 2  O(dx )  





1


Rin : 

____
____
dt  


2
2



q
dt

q
3 
   dt  


O
(
dt
)

2 
 

 t  2  t 


 

____
  ____


   q  u  q  



  t   x 



 


n
Ri : 

____
____


  dt   2 q 

dx   2 q 
2
2 
    2  O(dt )   u   2  O(dx )  
2  x 
  2  t 


 
12
Final Form
• Using the definition of q
____
  ____


   q  u  q  







  t   x 



 


n
Ri : 

____
____


2
  dt   2 q 


dx

q
2
2 
    2  O(dt )   u   2  O(dx )  
2  x 
  2  t 


 
____
____


2
2




dt  q
dx  q
Rin :    2  O(dt 2 )   u   2  O(dx 2 ) 
 2  t 

2  x 



 ____
____
2
2

q
Using that:  q
2

u
t 2
x 2
 ____
2 
udx  udt    q 
2
Rin :
1


O
(
dx
)

 2 
2 
dx   x 


13
Interpretation of Consistency
 ____
2 
udx
udt


 q 
n
2
Ri :
1


O
(
dx
)

 2 
2 
dx   x 


So the truncation error is O(dx) under the assumption that
dt/dx is a constant..
This essentially implies that the numerical solution diverges
from the actual solution by an error of O(dx) every time step.
If we assume that the solution q is smooth enough then the
truncation error converges to zero with decreasing dx. This
property is known as consistency.
14
Error Equation
• We define the error variable:
in  qin  Ein
• We next define the numerical iterator N:
 n1  N  n
• Then:
E
n 1
 N q  E
q
 N q  E   N q   N q   q
  N  q  E   N  q   dtR
n
n
n
n
n
n 1
n
n
n
n
n 1
n
• So the new error consists of the action of the numerical
scheme on the previous error and the error commited in the
approximation of the derivatives.
15
Abstract Scheme
• Without considering the specific construction of the
scheme suppose that the numerical N operator satisfies:
N  P   N Q   P  Q
• i.e. N is a contraction operator in some norm then…
16
Estimating Error in Terms of Initial Error and
Cumulative Truncation Error


E n1  N  q n  E n   N  q n   dtR n

E n1  N  q n  E n   N  q n   dt R n
triangle inequality

E n1  q n  E n  q n  dt R n
contraction property of N

E n1  E n  dt R n


 E n1  dt R n1  dt R n
 .....

E
n 1
mn
 E  dt  R m
0
m 1
by induction
17
Error at a Time T (independent of dx,dt)
E
n 1
mn
 E  dt  R m
0
m 1


E n1  E 0  T max R m
m 1,..n

• If the method is consistent (and the actual
solution is smooth enough) then:
E n1  E 0  T O  dx 
• As dx -> 0 the initial error ->0 and
consequently the numerical error at time T
tends to zero with decreasing dx (and dt).
18
Specific Case: Stability and Consistency for
the Upwind Finite Volume Scheme
• We already proved that the upwind FV
scheme is consistent.
• We still need to prove stability of.
in1  1    in     in1
19
Stability in the Discrete 1-norm
in1  1    in     in1


N 1
n 1
1
 dx in1
i 1
N 1
 dx 1    in     in1
triangle inequality
i 1
N 1
N 1
 dx 1      dx  in1 assuming 0    1
n
i
i 1
i 1
N 1
 dx   dx in
n
0
i 1
 dx 0n   n
1
So here’s the interesting story. In the case of a zero boundary
condition then we automatically observe that the operator is a
20
contraction operator.
Boundary Condition
• Suppose  , are two numerical solutions
with  n   n
0
• Then:
0
 n1   n1 1  dx 0n   0n   n   n
  n  n
1
1
• i.e. if we are spot on with the left boundary
condition the N iterator is indeed a
contraction.
21
Relaxation on Stability Condition
• Previous contraction condition on the
numerical iterator N
N  P   N Q   P  Q
• A less stringent condition is:
N  P   N  Q   1   dt  P  Q
• Where alpha is a constant independent of dt
as dt->0
22
Relaxation on Stability Condition
• In this case the stability analysis yields:


E n 1  N  q n  E n   N  q n   dtR n

E n 1  q n  E n  q n  dt R n

E n 1  1   dt  E n  dt R n


 1   dt  1   dt  E n 1  dt R n 1  dt R n
 .....

E
n 1
 1   dt 
 eT

n 1
mn
E  dt  1   dt 
0
m 1

E 0  T max R m
m 1,..n

nm
Rm
Ndt=T
23
Interpretation


E n1  eT E 0  T max R m
m1,..n
 , Ndt=T
• Relaxing the stability yields a possible
exponential growth – but this growth is
independent of T so if we reduce dt (and dx)
then the error will decay to zero for fixed T.
24
Next Lecture (6)
• Alternative flux formulations
• Alternative time stepping schemes
25