Transcript The Method of Fundamental Solutions and Compactly

Numerical Solution of Partial Differential Equations
Using Compactly Supported Radial Basis Functions
C.S. Chen
Department of Mathematics
University of Southern Mississippi
7/7/2015
1
PURPOSE OF THE LECTURE
TO SHOW HOW COMPACTLY SUPPORTED RADIAL
BASIS FUNCTIONS (CS-RBFs) CAN BE USED TO
EXTEND THE APPLICABILITY OF BOUNDARY
METHODS TO PROVIDE 'MESH FREE' METHODS
FOR THE NUMERICAL SOLUTION OF PARTIAL
DIFFERENTIAL EQUATIONS.
7/7/2015
2
7/7/2015
3
MESHLESS METHOD
MESH METHOD
7/7/2015
4
The Method of Particular Solutions
Consider the Poisson’s equation
Lu  f (x), x  ,
u  g (x), x  ,
(1)
(2)
Where   Rd , d  2,3, is a bounded open nonempty domain
with sufficiently regular boundary .
Let v  u  u p where u p satisfying Lu p  f (x) (3)
but does not necessary satisfy the boundary condition in (2).
v satisfies
7/7/2015
Lv  0, x  ,
u  g (x)  u p (x),
x  .
(4)
(5)
5
• Domain integral Domain integral
u p ( P)   G ( P, Q;  ) f (Q)dv

G( P, Q;  ) is a fundamental solutionfor
L
• Atkinson’s method (C.S. Chen, M.A. Golberg & Y.C. Hon, The
MFS and quasi-Monte Carlo method for diffusion equations, Int. J. Num.
Meth. Eng. 43,1421-1435, 1998)
ˆ 
u p ( P)  ˆ G( P, Q;  ) f ( P)dv, 

(Requires no meshing if ˆ = circle or sphere)
f
• Radial Basis Function Approximation of
• 7/7/2015
Others
6
The Method of Particular Solutions
Assume that
f (x)  fˆ (x)
and that we can obtain an analytical solution u
 p to
uˆ p  fˆ (x).
(6)
Then
u p  uˆ p .
To approximate f by f we usually require fitting the given
data set
x 
N
i 1
of pairwise distinct centres with the imposed
conditions
f (x)  fˆ (x),
7/7/2015
1  i  N.
7
The linear system

N

f (xi )   ai xi  x j ,
i 1
1 i  N,
(7)
is well-posed if the interpolation matrix is non-singular

A   xi  x j

Once f in (6) has been established,
1i  N
N
u p   ai i
(8)
i  i
(9)
i 1
where
and
7/7/2015


i   x  xi ,


i   x  xi .
8
Compactly Supported RBFs
Let : R  R be a continous function with  (0)  0. If xi , let


i xi    x  xi ,
where  is the Euclidean norm. Then i is called the RBF.
References
• Z. Wu, Multivariate compactly supported positive definite radial
functions, Adv. Comput. Math., 4, pp. 283-292, 1995.
• R. Schaback, Creating surfaces from scattered data using radial basis
functions, in Mathematical Methods for Curves and Surface,
eds. M. Dahlen, T. Lyche and L. Schumaker, Vanderbilt Univ. Press,
Nashville, pp. 477-496, 1995
• W. Wendland, Piecewise polynomial, positive definite and compactly
supported RBFs of minumal degree, Adv. Comput. Math., 4, pp 389396, 1995.
7/7/2015
9
References
• C.S. Chen, C.A. Brebbia and H. Power, Dual receiprocity method
using compactly supported radial basis functions, Communications in
Numerical Methods in Engineering, 15, 1999, 137-150.
• C.S. Chen, M. Marcozzi and S. Choi, The method of fundamental
solutions and compactly supported radial basis functions - a meshless
approach to 3D problems, Boundary Element Methods XXI, eds. C.A.
Brebbia, H. Power, WIT Press, Boston, Southampton, pp. 561-570,
1999.
• C.S. Chen, M.A. Golberg, R.S. Schaback, Recent developments of the
dual reciprocity method using compactly supported radial basis
functions, in: Transformation of Domain Effects to the Boundary, ed.
Y.F. Rashed, WIT PRESS, pp183-225, 2003.
• M.A. Golberg, C.S. Chen, M. Ganesh, Particular solutions of the 3D
modified Helmholtz equation using compactly supported radial basis
functions, Engineering Analysis with Boundary Elements, 24, pp. 539547, 2000.
7/7/2015
10
Wendland’s CS-RBFs
Define
n

(
1

r
)
, 0  r  1,
n
(1  r )   
r  1.
 0,
For d=1,
(10)
  (1  r )   C 0
  (1  r ) 3 (3r  1)  C 2
  (1  r ) 5 (8r 2  5r  1)  C 4
For d=2, 3,
F
o
r
7/7/2015
  (1  r ) 2 C 0
  (1  r ) 4 (4r  1) C 2
  (1  r ) 6 (35r 2  18r  3) C 4
  (1  r )8 (32r 3  25r 2  8r  1) C 6
11
Globally Supported RBFs
φ=1+r
7/7/2015
  r 2  c2
12
Compactly Supported RBFs
  (1  r )2
7/7/2015
  (1  r ) 4 (4r  1)
13
Compact support cut-off parameter (scaling factor)

7/7/2015
14
Analytic Particular Solutions L=  in 2D
r
 

4

r  r 
1    4  1 , r   ,
        

0,
r  .

1  d  d 
    r  ,
r  dr  dr 
in 2D,
 d  r 
 r  r

r      =  r 1    4  1 dr
    
 dr    
4r 7 5r 6 4r 5 5r 4 r 2
 5  4  3  2   A,
r  .
7
2

2
2
4
7/7/2015
15
 4r 6 5r 5 4r 4 5r 3 r

r
   =   5  4  3  2   C1  dr
2

2
2
 
 7

4r 7
5r 6
4r 5 5r 4 r 2 A


 3  2    B,
r  .
5
4
49 12
5 8
4 r
Choose A= B = 0, we have
 4r 7
5r 6
4r 5 5r 4 r 2

 3  2  , r  ,
r 
5
4
   =  49 12
5 8
4
  
C ln(r )  D,
r  ,

(*)
where C and D are to be determined by matching  and ' at r   .
7/7/2015
16
Note that
d  r   4r 6 5r 5 4r 4 5r 3 r
  =  5  4  3  2  , r  ,
dr     7
2

2
2

C

,
r 

r
4 6 5 5 4 4 5 3  C
For r   ,
 4 3  2 
5
7
2

2
2 
 C
2


 C
14 
14
From (*), we have
7/7/2015
529 2  2
D

ln 
5880 14
17
 4r 7
5r 6
4r 5 5r 4 r 2

 3  2  , r  ,

5
4
5 8
4
 r   49 12
 =
529 2  2  r 
  

ln   ,
r  ,

5880 14   

7/7/2015
18
Example in 2D
7/7/2015
19
7/7/2015
20
7/7/2015
21
Example 2D
5 2
 y 
u ( x, y ) 
sin( x) cos 
,
4
 2 
 y 
u ( x, y )  sin( x) cos 
,
 2 
( x, y)  ,
( x, y )  .
1
1
0.5
0.5
0
0
-0.5
-0.5
-1
-1
-1
-0.5
0
0.5
1
10,000 uniform grid points
7/7/2015
-1
-0.5
0
0.5
1
22
Sparse matrix for uniform grid points.
7/7/2015
23
 r  r 
 1    4    1
       
4
CS-RBF used:
MFS: Fictitious circle with r = 4.
Interpolation points: 100x100 uniform grid points
7/7/2015
α
RMSE
0.10
3.76E-4
561,013
1.65
0.20
1.55E-5
2,171,057
3.36
0.25
8.43E-6
3,376,025
5.32
NZ
C PU (sec)
24
Analytic Particular Solutions L=  in 3D
Recall
r
 

Since
4

r  r 
1    4  1 , r   ,
        

0,
r  .

1
 2
r
 d  2 d 
 dr  r dr  ,


(11)
in 3D,
we have
r2
r4
2r 5
5r 5
r7
 r   6  2 2  3 3  16 4  14 5 , r  
   
  
2 3

,
r  .

14 42r
7/7/2015
25
Numerical Example in 3D
Consider the following Poisson’s problem
u( x, y, z)  3 cos( x) cos( y) cos( z), ( x, y, z) ,
u( x, y, z)  cos( x) cos( y) cos( z), ( x, y, z) .
Physical Domain
  
0.5
0.25
0
-0.25
-1
-0.5
0.5
0.25
0
7/7/2015
1
0
-0.25
-0.5
26
4
We choose   1  r  (4r  1) to approximate the forcing term.
To evaluate particular solutions, we choose 300 quasi-random points
in a box [-1.5,1.5]x[-0.5,0.5]x[-0.5,0.5]. The numerical results are
compute along the x-axis with y=z=0.
8
Relative Errors (%)
7
6
 = 0.7
5
4
3
 = 1.0
 = 1.4
2
1
0
-1.0
-0.5
0.0
0.5
1.0
1.5
X Axis
7/7/2015
The effect of various scaling factor α
27
Modified Helmholtz Equation in 3D
n

r
1  d 2 d 
1   p (r ), 0  r   ,
2
r




  

2 
r  dr
dr 

0,
r  .

j
w(r )
d
 w
j
Let (r ) 
, r  0,  (0)  lim j   , j  0, 1, 2.
r 0 dr  r 
r
1
r2
7/7/2015
2
 d 2 d  1 d w
 r

 dr
dr  r dr 2
(15)
(16)
(17)
28
Hence, (15) is equivalent to
n
 
r r
d w
r 1   p  ,
2
 w  
 
dr 2

0,

2
0  r  ,
r  .
(18)
The general solution of (18) is of the form
 Ae  r  Be r  q (r ), 0  r   ,
w(r )  
 r
r
Ce

De
,
r  .

(19)
q(r) can be obtained by the method of undetermined coefficients,
or by symbolic ODE solver (MATHEMATICA or MAPLE).
A,B,C and D in (19) are to be chosen so that  in (16) is twice
differentiable at r = 0 and 
7/7/2015
29
Theorem 1. Let w be a solution of (18) with w(0)=0. Then 
defined by (16) is twice continuously differentiable at 0 with
 (0)=w'(0),  '(0)=0, ''(0)=[s2 w'(0)+p(0)]/3.
Furthermore,  (r) satisfies (15) as limr→0+

A  B  q (0)  0

Ae    Be   q ( )  Ce    De 

 Ae    Be   q ' ( )  Ce    De 

(20)
Choose D = 0. Consequently, we get
7/7/2015
A  [ B  q ( 0)]


e   ( q ' ( )  q ( )]

 B
2

2 

C

B
(
e

1
)

q
(

)
e
 q ( 0)


(21)
30
Hence, the particular solution Φ is given by

 s( 2 B  q ( 0)  q ' ( 0),
r  0,

 Ae  sr  Be sr  q ( r )
( r )  
,
0  r  ,
r

 sr
Ce

,
r  .

r

(22)
Notice that for r > α and large wave number λ
Ce  r e ( r ) q '( )  q ( ) q  e  ( r )  q (0)e  r 


r
2r
r
 0 for r   .
7/7/2015
31
Example
The general solution of
2
d 2w
r

2
 ,
2   w  r 1 

dr

0 r 
is given by
w(r ) 
4  r2 2  6r  2r 2 2  r 32
4 2
 Ae  r  Be r
Hence,
q (r ) 
4  r2 2  6r  2r 2 2  r 32
4 2
A,B,C can be obtain from (21).
7/7/2015
32
4
r  r


For  (r )  1    4  1 ,

    
1800 1 
 60
q (r )   6 3  8 5  r  4 2  6 4  2 
 



 
480
2880
300 
 240 1400  3  10
 r  4 3  6 5  r  2 2  4 4 
 
 
 
 
2
120 
15 5
4 6
 20
r  2 3  4 5  2 4 r  2 5 r
 
  

4
7/7/2015
33
Numerical Results
Consider
   400I  u  f ,
u  g,
in ,
on .
(23)
(24)
  ( x, y, z)  R3: H( x, y, z)  1 with
2
2

3 
3  2 2
H ( x , y , z)  min x   ,  x     y  z
4 
4 

7/7/2015
34
We choose f and g such that the u(x,y,z)=ex+y+z /400 is the
exact solution of (23)-(24).
To evaluate particular solutions, we choose N=400 quasirandom points . We choose the CS-RBFs
2
r

 (r )   1  
 
7/7/2015
35
Error estimates for various compact support cut-off parameter

0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
7/7/2015
||u-uN ||
2.2864E-02
1.6194E-02
9.7795E-03
6.4514E-03
4.5061E-03
3.8208E-03
3.3006E-03
2.9089E-03
2.6375E-03
CPU-time (sec.)
01.40
04.43
07.48
10.63
13.71
16.98
20.53
24.08
27.60
36
Multilevel Interpolation
Given a set of X scattered points, we decompose X into a nested
sequence
X1  X 2    X M  X
M subsets


X k  x1( k ) ,, x N( kk)  X , 1  k  M .
Algorithm:
S1 | X 1  f | X 1 ,
Nk

Sk ( x)   c(j k )k x  x (j k )
S2 | X 2  ( f  S1 )| X 2 ,
j 1


M 1
S M |X M  ( f 
7/7/2015
S
k 1
k
)| X M .
S
1
 S2  S M | X  f | X
37
The basic idea is to set α1 relatively large with few interpolation
points, and to let the αk decrease as k increases with more points.
Reference:
M.S. Floater, A. Iske, Multistep scattered data interpolation using compactly
supported radial basis functions, J. Comp. and Appl. Math., 73, 65-78, 1996.
30 interpolation points
100 interpolation points




7/7/2015


38
300 interpolation points

7/7/2015

39
Numerical Results
u  2e x  y ,
( x, y)  D,
u  e x  y  e x cos y,
( x, y) D,
where D is the Oval of Cassinii.
0.4
0.2
-1
-0.5
00
-0.2
0.5
1
-0.4
7/7/2015
40
α/1.5, 1.0, 0.6, 0.2/ n/30, 60, 200, 400/
7/7/2015
x
y
Numer.
Exact
Abs. Error
0.00
0.00
1.99994
2.00000
0.000061
-0.63
0.00
3.74097
3.74098
0.000012
1.34
0.00
7.65247
7.65261
0.000028
0.17
0.17
2.17206
2.35207
0.000082
0.00
0.20
1.79712
1.79707
0.000052
-0.37
0.37
1.11578
1.11594
0.000156
-1.34
0.00
0.5228
0.52270
0.000130
41
For time-dependent problems, we consider two
approaches to convert problems to Helmholtz
equation
• LAPLACE TRANSFORM
• FINITE DIFFERENCES IN TIME
ALSO POSSIBLE
• OPERATOR SPLITTING (RAMACHANDRAN AND
BALANKRISHMAN)
7/7/2015
42
(HEAT EQUATION)
Consider the BVP
u( P, t )  ut (P, t )  f ( P, t ), P  D  Rd , d  2,3
u ( P, t )  g ( P, t ), P  D, t  0
u ( P,0)  h( P,0), P  D
(25)
(26)
(27)
where D is a bounded domain in 2D and 3D.
Let

uˆ ( P, s)   est u( P, t )dt
0
(28)
Then uˆ ( P, s) satisfies
uˆ ( P, s)  suˆ ( P, s)  fˆ ( P, s)  h( P)  M ( P, s)
uˆ ( P, s)  gˆ ( P, s), P  D
7/7/2015
(29)
(30)
43
(WAVE EQUATION)
Consider BVP
u ( P, t )  utt ( P, t ), P  D, t  0
u ( P, t )  g ( P, t ), P  D, t  0
u ( P,0)  u0 ( P), ut ( P,0)  v0 ( P)
(31)
(32)
(33)
Then uˆ ( P, s) satisfies
uˆ ( P, s)  s 2uˆ ( P, s)  su0 ( P)  v0 ( P)
uˆ ( P, s)  gˆ ( P, s), P  D
(34)
(35)
Similar approach can be applied to hyperbolic-heat
equation and heat equations with memory
7/7/2015
44
(HEAT EQUATION)
Let tn  n , n  0. For tn  t  tn1
u( P, t )  u( P, tn1 )  (1   )u( P, tn )
u( P, t )  u( P, tn1 )  (1   )u( P, tn )
and
ut (P, t )  [u(P, tn1 )  u( P, tn )]/ 
Then vn ( P)  u( P, tn ) satisfies
vn1  vn1 /   vn /   (1  )vn /   f n /   mn (P)
For   1 (Rothe's Method)
vn1  vn1 /   vn /   f n
For  1/2 (Crank- Nicholson)
vn1  2vn1 /   2vn /   vn  2 f n  mn (P)
7/7/2015
(36)
(37)
(38)
(39)
(40)
(41)
45
(J. Su & B. Tabarrok, A time-marching integral equation method for unsteady
state problems, Comp. Meths. Appl. Mech. Eng. 142, 203-214, 1997)
(WAVE EQUATION)
utt  (un1  2un  un1 ) /  2
(41)
Then vn  un satisfies
vn1  vn1 /  2  (vn1  2vn ) /  2  (1  )vn1
(42)
CONVECTION-DIFFUSION EQUATION


u  V  u  ut (V  velocityfield)
(43)
(44)
ut  (un 1  un 1 ) / 2

vn1  (1  )vn1  (vn1  vn1 ) / 2  V  vn
(45)
Similar approach works for non-linear equation
u( P, t )  ut ( P, t )  f ( P, u, u)
7/7/2015
(46)
46