Transcript MA375 - Rice U - Computational and Applied Mathematics

MA557/MA578/CS557
Lecture 23
Spring 2003
Prof. Tim Warburton
[email protected]
1
Basis Ordering
There has been some confusion about the numbering of the
PKDO basis functions. There are now a number of basis
functions. We will choose an ordering:
 1  0,0
 p 2  0,1
 2  1,0 
 p 3  1,1
 3   2,0 
 p 1  
 M  0,p 
 2 p 2   p,1
p ,0 
where M=(p+1)(p+2)/2 and we have normalized each of the { i }
2
Accuracy of Polynomial Interpolation
• We introduce the operator supremum-norm and the
Linfinity norm:
I p  sup
f 0
Ip f
f


where f

 max f  x, y 
 x , y T
for the p’th order interpolation operator.
• The interpolation error is bounded below by:
f  x, y   f *  x, y 

 f  x, y   I p f  x, y 

where f * is the optimal polynomial approximation
in the .  norm
3
cont
• This tells us that the interpolating polynomial which
approximates a function can only be less or equal in
accuracy to the best polynomial approximation.
• Finding the optimal polynomial which approximates a
function is an unsolved problem.
• However, there is a theorem by Lebesgue which
allows us to bound the interpolation error in terms of
the error one would attain by choosing the best
possible polynomial approximation.
4
Theorem (Lebesgue)
• Assume that f  C T  and we consider a set of M
interpolation points then:
f Ipf

 1   p  f  f *
where  p  I p  max
 x , y T

iM
 h  x, y  is
i 1
i
termed the Lebesgue constant
• i.e. the Lebesgue constant gives us an idea of the
p
optimality of the M points. The smaller  is the
better.
5
Proof
f Ipf

 f  f*
 ff
 ff

 f*Ipf

triangle inequality
nM
*
 f   f  xn , yn  hn  x, y 
*

*

n 1

nM
definition of interpolant

nM
 f  x , y  h  x, y    f  x , y  h  x , y 
*
n
n
n
n
n 1
n
n
n 1

ok since f *  P p  T 
 ff
 ff
 ff

 1 

*

*

*
nM

n 1


nM
  f  x , y   f  x , y  h  x, y 
*
n
n
n
n
n
n 1

nM


f *  xn , yn   f  xn , yn  hn  x, y 
n 1
 ff

nM
*

 h  x, y 
n
n 1

hn  x, y   f  f *



6
Comments
• Note – the Lebesgue number only depends on the
distribution of nodes. The only condition on the function
being approximated is that it is continuous.
p

• Here’s the tricky part. It is difficult to compute
exactly
as it requires us to search the continuum of points in the
triangle.
• However, we can approximate this number by randomly
sampling a large number of points in the triangle. i.e. for
some, large, N compute the following:
iM
 p  max  hi  x j , y j 
j 1,.., N
i 1
for some finite subset
 x , y 
j
j
j 1,.., N
T
7
Open Problem
• It is still an open problem to find the set of nodes for the
triangle which minimize the Lebesgue constant.
• Various attempts have been made:
• Hesthaven:
• http://epubs.siam.org/sam-bin/getfile/SINUM/articles/30587.pdf
(You will need access permission – so download on on campus)
• Wingate and Taylor:
• http://citeseer.nj.nec.com/taylor98fekete.html
• And others.. See comparison tables in Hesthaven paper.
8
Fekete Nodes For the Triangle
1) Read the Wingate and Taylor paper over the break.
2) The Fekete nodes are chosen according to:
the condition that they are the nodes which maximize the
determinant of the VDM matrix.
3) Wingate and Taylor suggest a method for calculating the
Fekete nodes based on the method of steepest ascent.
4) Nodes are moved in the direction which most increase the
log of the determinant of the generalized Vandermonde
matrix.
9
Steepest Descent Approach
Suppose we start with a list of M initial node locations we will
solve the following ODE to steady state. We denote the
determinant of the Vandermonde matrix by: V
   
 rn , sn    ,   log V  r1, s1, , rM , sM  
t
 ri si 

We will now show that this is equivalent to solving:
hn
 hn

 rn , sn   
 rn , sn  ,  rn , sn  
t
s
 r


10
We need to compute the gradient of log(V)
   
 rn , sn    ,   log V  r1, s1, , rM , sM  
t
 ri si 

Step 1:
expand the determinant of the Vandermonde matrix
with respect to the co-determinants:
M
V    1  n  r1, s1  V1n
n
n 1
11
Step 2: as an example of how to compute the gradient of |V|
with respect to the node coordinates we first
differentiate with respect to r1
M

n
 r1  V    r1    1  n  r1, s1  V 
 n 1

M
   1
n 1
n
     r , s   V
r
n
   r 1  r1, s1 

 1  r2 , s2 



   r ,s 
 1 M M
1
1
1n
  r M  r1, s1  

 M  r2 , s2  
 M  rM , sM 



12
Step 3: Expand derivative terms with respect to
the Lagrange interpolating basis for the M nodes
 n  r ,s  
m M
 h  r ,s   r
m 1
m M
m
n
m
, sm 
  r n  r , s      r hm  r , s   n  rm ,sm 
m 1
13
Step 4: replace derivative terms in differentiation
of V:


   r ,s
 r 1 1 1
  1 r2,s2
r  V   
1

  r ,s
 1 M M





 r M  r1,s1 
 M  r2,s2  

 M  rM ,sM 
 m M
    r hm  r1, s1   1  rm , sm 
 m1
 1  r2 , s2 
 


 1  rM , sM 





s
,
r

s
,
r
h

 r m  1 1  M  m m  

m 1

 M  r2 , s2 



 M  rM , sM 

m M
14
Step 5: factorize matrix
  r h1  r1, s1   r h2  r1, s1 
 r hM  r1, s1  


0
1
0
r  V   

1




1
0
0


  1  r1, s1 
 M  r1, s1  


s
,
r

s
,
r



M  2 2 
 1 2 2



  r , s 
s
,
r



M 
M
M
 1 M M

r h1  r1,s1  V
15
Conclusion:
 V 
r1
 h1 
 V
 r1, s1 

 r 
• This is not exactly as reported in the Wingate &
Taylor paper.
• However, since |V| simply scales the node
velocities uniformly it will not make very much
difference as we are integrating in a “pseudo-time”
manner.
16
Interpretation of Algorithm
At each node, the algorithm pushes the node in the
direction of the slope of the interpolating polynomial
at that node:
  rn 
 s   hn  rn , sn 
t  n 
At steady state, the interpolating polynomial
associated with a node will have (approximately)
zero slope at that node..
This is a property shared by the Gauss-LobattoLegendre nodes in 1D !.
17
Download The Code
• Download the Fekete node routine computing routine from
the class web page.
• Over the break use this code to compute the Fekete nodes
for p=1 to p=10.
• Save the node coordinates in a file for each polynomial
order..
• These will be the nodes we use for future computations.
18
To Run:
• Unzip the files and set the current directory to the
umFEKETE directory.
• In Matlab:
p = 7;
[R,S] = am282fekete2d(p);
scatter(R,S);
• This might take quite a while. The code will output the
approximate Lebesgue constant as the nodes move.
• Run now with p=4
19
P=7 nodes
• Note there are 8 nodes on each edge and a total of 36 nodes.
• The approximate Lebesgue number is: 4.92707087069006
20
am282fekete2d
am282derivmatrix2d
 
ˆ V 1
Dr  D
r
am282ab
a 
2 1  r 
b s
1 s
1
am282vdm2dab
V  a, b 
am282jacobideriv2dab
am282jacobi2dab
 n  a, b 
am282jacobi1d
Pn , ( x )
   a, b  ,   a, b  


s
 r

n
n
am282jacobderiv1d
dPn ,
(x)
dx
21
Surface Mass Matrices
• Recall the DG scheme:
Find Cm  t  , m  1,..., M  such that:
m M
 h ,h 
m 1
n
m T

 mM   hm 
dCm mM   hm 
  u  hn ,
Cm .....
 Cm    v  hn ,

dt
x T  m1  
y T 
m 1  
     
  hn , 
  2
 


 C    0



T
where    u  n  ,
C   C   C 
• We now have defined the set of nodes which we can
build the hn Lagrange interpolating polynomials.
• The only remaining task is to discretize the boundary
terms.
22
Surface Term
• We now break the boundary integral into the three
edge components:
     
 hn , 
  2
 

m  M f 3 

  f   f
 C       hn , 



2
m 1 f 1 


T


m  M f 3
 h , h 
m 1 f 1
n
m T
f
 
 hm  Cm 
 
 T f
  f   f


2


  Cm 


• Last step due to the face normals being constant
• The last step is to compute:
 hn , hm T
f
23
 hn , hm T
f
• Starting with edge 1 (-1<=a<=1,b=-1), we will work with
the PKDO basis:
  ,  
ij
kl
Tˆ1
1
   ij   a, 1 kl   a, 1 da
1
i
k
 0,0


1

(

1)
1

(

1)
 0,0

 2i 1,0

 2 k 1,0
   Pi  a  
 1 Pk  a  
 1da
 Pj
 Pl
 2 
 2 


1 


1
1
 Pj2i 1,0  1 Pl 2 k 1,0  1 Pi 0,0  a  Pk0,0  a da
1
2 i 1,0
j
P
 1 Pl
2 k 1,0
2
 1  ik
2i  1
24
  ,  
ij
kl
Tˆ1
2 i 1,0
j
P
 1 Pl
2 k 1,0
2
 1  ik
2i  1
i  p j  p i
nM
ij   r , s    ij   rn , sn  hn  r , s   hn  r , s   
n 1
i 0
 V    
1
n ij
j 0
ij

 hn , hm Tˆ
1
k  p l  p k
 i  p j  p i 1

1
     V   ij  ,    V   kl  
n ij 
n kl 
k 0 l 0
 i 0 j 0
Tˆ1
i  p j  p i k  p l  p  k

i 0

j 0 k 0 l 0
1
V
 
i  p j  p i k  p l  p  k

i 0
i 0
n kl 
   V    V  
1
j 0 k 0 l 0
i  p j  p  i l  p i

n ij 
1
V
 

j 0
  V 1 
l 0
n ij 
1
  ,  
ij
2 i 1,0
j
P
kl
Tˆ1
 1 Pl
2 k 1,0
2
 1  ik
2i  1
n ij
n kl 
 V 1 
Pj2i 1,0  1 Pl 2i 1,0  1
n il 
2
2i  1
25
Edge 2
• If we try the same thing with the second edge it gets
more complicated:
  ,  
ij
kl
Tˆ2
1
   ij  1, b   kl  1, b  db
1
 0,0  1  b i 2i 1,0
  0,0  1  b  k 2 k 1,0

   Pi 1 
 b   Pk 1 
 b db
 Pj
 Pl
 2 
 2 
 

1 

1
 1  b i k 2i 1,0

2 k 1,0
  
 b  Pl  b db  ...?
 Pj

1 
 2 
1
26
Edges 2 & 3
• We can use the surface mass matrix from edge 1 by noting
that the Fekete nodes have the same distribution on edge
1 as they do on edge 2.
• Identify which nodes correspond and permute the edge 1
matrix.
27
Computing Jump Terms
• Note we need [Cm] which only makes sense for the
nodes which lie on the f’th edge.
28
Next Lecture
• Next class we will start examining the consistency of
the DG scheme
• Also – finally the consistency analysis.
• Review the Hesthaven-Warburton paper for details.
29