Numerical Methods for PDEs
Integral Equation Methods, Lecture 1
Discretization of Boundary Integral Equations
Notes by Suvranu De and J. White
April 23, 2003
1 Outline for this Module
Slide 1
Overview of Integral Equation Methods
Important for many exterior problems
(Fluids, Electromagnetics, Acoustics)
Quadrature and Cubature for computing integrals
One and Two dimensional basics
Dealing with Singularities
1
st
and 2
nd
Kind Integral Equations
Collocation, Galerkin and Nystrom theory
Alternative Integral Formulations
Ansatz approach and Green’s theorem
Fast Solvers
Fast Multipole and FFT-based methods.
2 Outline
Slide 2
Integral Equation Methods
Exterior versus interior problems
Start with using point sources
Standard Solution Methods
Collocation Method
Galerkin Method
Some issues in 3D
Singular integrals
3 Interior Vs Exterior Problems
Slide 3
Interior Exterior
Temperature
known on surface
2
0T?=
inside
2
0T?=
outside
Temperature
known on surfac
"Temperature in a tank" "Ice cube in a bath"
What is the heat ?ow?
Heat ?ow = Thermalconductivity
integraltext
surface
?T
?n
Note 1 Why use integral equation methods?
For both of the heat conduction examples in the above ?gure, the temperature,
T, is a function of the spatial coordinate,x, and satis?es ?
2
T(x)=0.Inboth
1
problemsT(x) isgivenonthesurface,de?nedby Γ,andthereforebothproblems
are Dirichlet problems. For the “temperature in a tank” problem, the problem
domain, ? is the interior of the cube, and for the “ice cube in a bath” problem,
the problem domain is the in?nitely extending region exterior to the cube. For
suchanexteriorproblem,oneneedsanadditionalboundaryconditiontospecify
whathappens su?ciently farawayfromthe cube. Typically,itis assumedthere
are no heat sources exterior to the cube and therefore
lim
bardblxbardbl→∞
T(x) → 0.
For the cube problem, we might only be interested in the net heat ?ow from
the surface. That ?ow is given by an integral over the cube surface of the
normal derivative of temperature, scaled by a thermal conductivity. It might
seemine?cienttousethe?nite-elementor?nite-di?erencemethodsdiscussedin
previous sections to solve this problem, as such methods will need to compute
the temperature everywhere in ?. Indeed, it is possible to write an integral
equation which relates the temperature on the surface directly to its surface
normal, as we shall see shortly.
In the four examples below, we try to demonstrate that it is quite common
in applications to have exterior problems where the known quantities and the
quantities of interest are all on the surface.
4Examples
4.1 Computation of Capacitance
Slide 4
v
+
-
2
0 Outsi?Ψ=
is given on SΨ
potential
What is the capacitance?
Capacitance = Dielectric Permittivity
integraltext
?Ψ
?n
Note 2 Example 1: Capacitance problem
In the example in the slide, the yellow plates form a parallel-plate capacitor
with an applied voltageV. In this 3-D electrostatics problem, the electrostatic
potential Ψ satis?es ?
2
Ψ(x)=0in the region exterior to the plates, and the
potential is known on the surface of the plates. In addition, far from the plates,
Ψ → 0. What is of interest is the capacitance,C, which satis?es
q =CV
2
where q, the net charge on one of the plates, is given by the surface normal of
the potential integrated over one plate and scaled by a dielectric permittivity.
4.2 Drag Force in a Microresonator
Slide 5
Resonator Discretized Stru
Computed Forces
Bottom View
Computed Forces
Top View
Note 3 Example 2: Drag force in a MEMS device
The example in the slide is a microresonator,it is a structure that can be made
to vibrate using electrostatic forces. The changing character of those vibra-
tions can be used to sense rotation. The particulars of how the microresonator
operates is not directly relevant to our discussion of integral equations, except
for one point. In order to determine how much energy is needed to keep the
microresonator vibrating, it is necessary to determine the ?uid drag force on
comb structures shown in the bottom part of theslide. The ?uid is the air sur-
roundingthe structure,andat the micron-scaleofthese devices, airsatis?es the
incompressible Stokes equation,
?
2
u(x)=?p(x) (1)
?·u(x)=0
whereuisthe?uidvelocityandpisthepressure.Byspecifyingthecombvelocity,
and then computing the surface pressure and the normal derivative of velocity,
onecandeterminethenetdragforceonthecomb. Onceagain,thisisaproblem
in which the known quantities and the quantities of interest are on the surface.
3
4.3 Electromagnetic Coupling in a Package
Slide 6
Picture Thanks to Coventor.
Note 4 Example 3: Electromagnetic coupling in a package
In the 40 lead electronic package pictured in the slide, it is important to de-
termine the extent to which signals on di?erent package leads interact. To
determine the magnetic interaction between signal currents ?owing on di?erent
wires, one must solve
?
2
H(x)=J(x) (2)
?·J(x)=0
whereH is the magnetic ?eld andJ is the signal current density. By specifying
the current, and then computing the magnetic ?eld at the surfaces of the leads,
one can determine the magnetic interaction. Again, this is a problem in which
the known quantities and the quantities of interest are on the surface.
4.4 Capacitance of Microprocessor Signal Lines
Slide 7
Note 5 Example 4: Capacitance of microprocessor signal lines
Thislastexampleintheaboveslideisapictureofthewiringonamicroprocessor
integratedcircuit. Atypicalmicroprocessorhasmillions ofwires,soweareonly
looking at a small piece of a processor. The critical problem in this example
4
is determining how long signals take to get from the output of a logical gate
to the input of the next gate. To compute that delay, one must determine the
capacitance on each of the wires given in the slide picture. To do so requires
computing charges given electrostatic potentials as noted above.
5 What is common about these problems?
Slide 8
Exterior Problems
MEMS device - fluid (air) creates drag
Package - Exterior fields create coupling
Signal Line - Exterior fields.
Quantities of interest are on surface
MEMS device - Just want surface traction force
Package - Just want coupling between conductors
Signal Line - Just want surface charge.
Exterior problem is linear and space-invariant
MEMS device - Exterior Stoke’s flow equation (linear)
Package - Maxwell’s equations in free space (linear)
Signal line - Laplace’s equation in free spce (linear)
But problems are geometrically very complex
6 Exterior Problems
6.1 Why not use FDM / FEM?
Slide 9
2-D Heat Flow Example
0 at T =∞
But, must
truncate t
mesh
Surface
Onlyneed
?T
?n
onthesurface,but T iscomputedeverywhere.
Musttruncatethemesh,? T(∞)=0becomes T(R)=0.
Note 6 Heat conduction in 2D
In this slide above, we consider a two dimensional exterior heat conduction
probleminwhichthetemperatureisknownontheedges,orsurface,ofasquare.
Here, the quantity of interest might be the total heat ?ow out of the square.
The temperatureT satis?es
?
2
T(x)=0x∈ ? (3)
5
T(x) givenx∈ Γ
lim
bardblxbardbl→∞
T(x)=0
where ? is the in?nite domain outside the square and Γ is the region formed by
the edges of the square.
Using ?nite-element or ?nite-di?erence methods to solve this problem requires
introducing an additional approximation beyond discretization error. It is not
possible to discretize all of ?, as it is in?nite, and therefore the domain must be
truncated with an arti?cial ?nite boundary. In the slide,the arti?cial boundary
is a large ellipse on which we assume the temperature is zero. Clearly, as the
radiusoftheellipseincreases,thetruncatedproblemmoreaccuratelyrepresents
thedomainproblem,butthenumberofunknownsinthediscretizationincreases.
6
7 Laplace’s Equation
7.1 Green’s Function
Slide 10
In 2D
Ifu=log
parenleftBig
radicalbig
(x?x
0
)
2
+(y?y
0
)
2
parenrightBig
then
?
2
u
?x
2
+
?
2
u
?y
2
=0 ? (x,y) negationslash=(x
0
,y
0
)
In 3D
Ifu=
1
√
(x?x
0
)
2
+(y?y
0
)
2
+(z?z
0
)
2
then
?
2
u
?x
2
+
?
2
u
?y
2
+
?
2
u
?z
2
=0 ? (x,y,z) negationslash=(x
0
,y
0
,z
0
)
Proof: Just di?erentiate and see!
Note 7 Green’s function for Laplace’s equation
In the next few slides, we will use an informal semi-numerical approach to
deriving the integral form of Laplace’s equation. We do this inpart because
such a derivation lends insight to the subsequent numerical procedures.
To start, recall from basic physics that the potential due to a point charge is
related only to the distance between the point charge and the evaluation point.
In 2-D the potential is givenby the log of the distance, and in 3-D the potential
is inversely proportion to the distance. The precise formulas are given on the
slide. A little more formally, direct di?erentiation reveals that
u(x,y)=log
radicalbig
(x?x
0
)
2
+(y?y
0
)
2
(4)
satis?es the 2-D Laplace’s equation everywhere except x=x
0
, y =y
0
and
u(x,y,z)=
1
radicalbig
(x?x
0
)
2
+(y?y
0
)
2
+(z?z
0
)
2
(5)
satis?es the 3-D Laplace’s equation everywhere except x = x
0
, y = y
0
, z = z
0
.
These functions are sometimes referred to as Green’s functions for Laplace’s
equation.
triangleright Exercise 1 Show by direct di?erentiation that the functions in (4) and (5)
satisfy ?
2
u=0, in the appropriate dimension almost everywhere.
7
8 Laplace’s Equation in 2D
8.1 Simple idea
Slide 11
22
22
+ = 0 outs
uu
xy
??
??
Surface
()
00
,xy
22
22
+ = 0 outside
uu
xy
??
??
Problem Solved
is given on surfaceu
()()
22
00
Let loguxxyy=?+?
Does not match boundary conditions!
Note 8 Simple idea for solving Laplace’s equation in 2D
Here is a simple idea for computing the solution of Laplace’s equation outside
the square. Simply let
u(x,y)=αlog
radicalbig
(x?x
0
)
2
+(y?y
0
)
2
wherex
0
, y
0
is a point inside the square. Clearly such auwill satisfy ?
2
u=0
outsidethesquare,butumaynotmatchtheboundaryconditions. Byadjusting
α, it is possible to make sure to match the boundary condition at at least one
point.
triangleright Exercise 2 Suppose the potential on the surface of the square is a constant.
Can you match that constant potential everywhere on the perimeter of the
square by judiciously selectingα.
8.1.1 "More points"
Slide 12
22
22
+ = 0 out
uu
xy
??
??
()
11
,xy
()
22
,x y
is given on surfaceu
(),
nn
xy
Letu=
summationtext
n
i=1
α
i
log
parenleftBig
radicalbig
(x?x
i
)
2
+(y?y
i
)
2
parenrightBig
=
summationtext
n
i=1
α
i
G(x?x
i
,y?y
i
)
Pick the α
i
’s to match the boundary conditions!
8
Note 9 ...contd
To construct a potential that satis?es Laplace’s equation and matches the
boundary conditions at more points, let u be represented by the potential due
to a sum of n weighted point charges in the square’s interior. As shown in the
slide, we can think of the potential due to a sum of charges as a sum of Green’s
functions. Of course, we have to determine the weights on then point charges,
and the weight on the i
th
charge is denoted herebyα
i
.
8.1.2 "More points equations"
Slide 13
()
11
,xy
()
22
,xy
(),
nn
xy
()
11
,
tt
xy
Source Strengths selec
to give correct potent
testpoints.
()( )
()( )
()
()
11 11 1111 1
11
,,,
,,,
nn nn nn
tt tntn tt
n
tt t tn tt
Gx xy y Gx xy y x y
Gx xy y Gx xy y x y
α
α
g170g186g170g186?? ?? Ψ
g170g186
g171g187g171g187
g171g187
g171g187
=
g171g187g171g187
g171g187
g171g187
g171g187g171g187
g171g187?? ?? Ψ
g172g188
g172g188g172g188
g22g22
g23g25 g23 g23g23
g23g25g23 g23g23
g22g22
Note 10 ...contd
To determine a system of n equations for the nα
i
’s,consider selecting a set of
ntest points, as shown in the slide above. Then, by superposition, for each test
pointx
t
i
,y
t
i
,
u(x
t
i
,y
t
i
)=
n
summationdisplay
i=1
α
i
log
radicalbig
(x
t
i
?x
0
)
2
+(y
t
i
?y
0
)
2
=
n
summationdisplay
i=1
α
i
G(x
t
i
?x
0
,y
t
i
?y
0
).
(6)
Writing an equation like (6) for each test point yields the matrix equation on
the slide.
The matrixA in the slide has some properties worth noting:
? Aisdense,thatisA
i,j
never equals zero. This is because every charge
contributes to every potential.
? If the test points and the charge points are ordered so that the i
th
test
point is nearest thei
th
charge, thenA
i,i
will be larger thanA
i,j
for allj.
Item 2 above seems to suggest that A is diagonally dominant, but this is not
thecase. Diagonaldominancerequiresthat theabsolute sum oftheo?-diagonal
entries is smaller than the magnitude of the diagonal. The matrix above easily
violates that condition.
9
triangleright Exercise 3 Determine a set of test points and charge locations for the 2-D
squareproblemthatgeneratesanAmatrixwherethemagnitudeofthediagonals
are bigger than the absolute value of the o?-diagonals, but the magnitude of
the diagonal is smaller than the absolute sum of the o?-diagonals.
8.1.3Computational Results
Slide 14
R=10
r
Circle with Chargesr=9.5
n=20 n=40
Potentials on the Circle
Note 11 results
It is possible to construct a numericalscheme for solving exterior Laplace prob-
lems by adding progressivelymore point charges so as to match more boundary
conditions. In the slide above, we show an example of using such a method to
compute the potential exterior to a circle of radius 10, where the potential on
the circle is given to be unity. In the example, charges are placed uniformly
on a circle of radius 9.5, and test points are placed uniformly on the radius 10
circle. If 20 point charges are placed in a circle of radius 9.5, then the potential
produced will be exactly one only at the 20 test points on the radius 10 circle.
The potential produced by the twenty point charges on the radius 10 circle is
plotted in the lower left corner of the slide above. As might be expected, the
potential produced on the radius 10 circle is exactly one at the 20 test points,
but then oscillates between 1 and 1.2 on the radius 10 circle. If 40 charges and
test points are used, the situation improves. The potential on the circle still
oscillates, as shown in the lower right hand corner, but now the amplitude is
only between 1 and 1.004.
10
8.2 Integral Formulation
8.2.1 Limiting Argument
Slide 15
Want to smear point charges to the surface
Results in an integral equation
Ψ(x)=
integraldisplay
surface
G(x,x
prime
)σ(x
prime
)dS
prime
Howdowesolvetheintegralequation?
Note 12 Single layer potential
The oscillating potential produced by the point charge method is due to the
rapid change in potential as the separation between evaluation point and point
chargeshrinks. If the point chargescould be smeared out, so that the produced
potential did not rise to in?nity with decreasing separation, then the resulting
computed potential would not have the oscillation noted on the previous slide.
Inaddition,itmakesthemostsensetosmearthepointchargesontothesurface,
as then the charge density and the known potential have the same associated
geometry. The result is the integral equation given on the bottom of the above
slide,wherenowtheunknownisachargedensity onthesurfaceandthepotential
due to that chargedensity is given by the well-knownsuperposition integral. In
the case of two or three dimensional Laplace problems,G(x,x
prime
) can be written
as
?
G(x?x
prime
), asthe potentialis onlya function of distance to the chargedensity
and not a function of absolute position. For such a Green’s function,
Ψ(x)=
integraldisplay
surface
?
G(x?x
prime
)σ(x
prime
)dS
prime
,
which one may recognize from system theory as a convolution. This connection
is quite precise. A space-invariant system has an impulse response, which is
usually referredto as a Green’s function. The output, in this case the potential,
is a convolution of the impulse response with the input, in this case the charge
density. Such an integral form of the potential is referred to as a single layer
potential.
Note 13 Types of integral equations
Thesinglelayerpotentialisanexampleofaclassofintegralequationsknownas
"Fredholm integral equation of the First Kind". A Fredholm integral equation
11
of the Second Kind results when the unknown charge density exists not only
under the integral sign but also outside it. An example of such an equation is
Ψ(x)=σ(x)+
integraldisplay
surface
K(x?x
prime
)σ(x
prime
)dS
prime
,
Fredholm integral equations, in which the domain of integration is ?xed, usu-
ally arise out of boundary value problems. Initial value problems typically give
rise to the so-called Volterra integral equations, where the domain of integra-
tion depends on the output of interest. For example, consider the initial value
problem
dx(t)
dt
=tx(t); t∈ [0,T],T>0.
x(t=0)=x
0
The "solution" of this equation is the following Volterra integral equation:
x(t)=x
0
+
integraldisplay
t
0
ξx(ξ)dξ
8.3 Basis Function Approach
8.3.1 Basic Idea
Slide 16
() ()
g44
1
Basis Functions
Represent
n
ii
i
xxσα?
=
=
g166
The basis functions are “on” the surfac
Example Basis
Represent circle with straight lines
Assume is constant along each lineσ
Can be used to approximate the density
May also approximate the geometry.
Note 14 Numerical solution of the single layer potential
As we have studied extensively in the ?nite-element section of the course, one
approach to numerically computing solutions to partial di?erential equations is
to represent the solution approximately as a weighted sum of basis functions.
Then, the original problem is replaced with the problem of determining the
basis function weights. In ?nite-element methods, the basis functions exist in
a volume, for integral equations they typically exist on a surface. For 2-D
problems that means the basis functions are restricted to curves and in 3-D the
basic functions are on physical surfaces.
Asanexample, considerthecircleintheslideabove. Onecouldtrytorepresent
the charge density on the circle by breaking the circle intonsub-arcs, and then
12
assumethechargedensityisaconstantoneachsub-arc. Suchanapproachisnot
commonly used. Instead the geometry is usually approximated along with the
charge density. In this example case, shown in the center right of the slide, the
sub-arcsofthecirclearereplacedwithstraightsections,thusformingapolygon.
The charge density is assumed constant on each edge of the polygon.The result
is a piecewise constant representation of the charge density on a polygon.
8.3.2 Geometric Approximation is Not New
Slide 17
Piecewise Straight surface basis
Functions approximate the circle
Triangles for 2-D F
approximate the circ
() ( ) ( )
1
,
n
ii
i
approx
surface
xGx xdSα?
=
′′′Ψ=
g166
g179
Note 15
The idea that both the geometry and the unknown charge density has been ap-
proximatedisnotactuallyanewissue. Asshowninthe?gureintheaboveslide,
if FEM methods are used to solve an interior problem, and triangular elements
are used, then the circle is approximated to exactly the same degree as when
straightsectionsreplacethesub-arcsforthesurfaceintegralequation. Asshown
at the bottom of the above slide, we can substitute the basis function represen-
tation into the integral equation, but then we should also note that the integral
is now over the approximated geometry. It is common, but not mathematically
justi?ed, to ignore the errors generated by the geometry approximation. We
will also ignore the error in the geometric approximation, and justify ourselves
by assuming that there are no circles, but only polygons, as then there is no
geometric approximation.
8.3.3 Piecewise Constant Straight Sections Example
Slide 18
1) Pick a set of n Point
surface
1
x
2) Define a new surface
connecting points with n
()
i
3) Define 1if is
i
xx l? =
()
i
otherwise, 0 x? =
2
x
n
x
1
l
2
l
n
l
13
Ψ(x)=
integraldisplay
approx
surface
G(x,x
prime
)
n
summationdisplay
i=1
α
i
?
i
(x
prime
)dS
prime
=
n
summationdisplay
i=1
α
i
integraldisplay
line l
i
G(x,x
prime
)dS
prime
How do we determine the α
i
’s?
Note 16
Intheaboveslide,wecompletethedescriptionofusingconstantchargedensities
on straight sections as the basis. If we substitute this example basis function
into the integral equation, as is done at the bottom of the slide, the result is to
replace the original integration of the product of the Green’s function and the
density with a weighted sum of integrals over straight lines of just the Green’s
function. The next step is then to develop an approach for determining the
weights, denoted here byα
i
’s.
8.3.4 Residual De?nition and Minimization
Slide 19
R(x)=Ψ(x)?
integraldisplay
approx
surface
G(x,x
prime
)
n
summationdisplay
i=1
α
i
?
i
(x
prime
)dS
prime
We pick the α
i
’s to make R(x) small
General Approach: Pick a set of test functionsφ
1
,...,φ
n
and forceR(x) to be
orthogonal to the set
integraldisplay
φ
i
(x)R(x)dS =0 foralli
Note 17
One way of assessing the accuracy of the basis function based approximation
of the charge density is to examine how well the approximation satis?es the
integral equation. To be more precise, we de?ne the residual associated with
the integralequation and an approximate solution at the top of the above slide.
Note thatR(x) is just the di?erence between the given potential on the surface
and the potential produced by the approximated chargedensity. Note also that
the equation is now over the approximate geometry and thereforex andx
prime
are
both on the approximated surface.
If the representation satis?es the integral equation exactly, then the residual
R(x) will be zero for all x and the approximate solution is equal to the exact
solution (provided the integral equation has a unique exact solution, more on
this later). However this will usually not be possible, and instead we will try
to pick the basis function weights, the α
i
’s, to somehow minimize R(x).One
14
approach to minimizing R(x) is to make it orthogonal to a collection of test
functions, which may or may not be related to the basis functions. As noted
on the bottom of the slide, enforcing orthogonality in this case means ensuring
that the integral of the product ofR(x) andφ(x) over the surface is zero.
15
8.3.5 Residual Minimization Using Test Functions
Slide 20
integraltext
φ
i
(x)R(x)dS =0 ?
integraldisplay
φ
i
(x)Ψ(x)dS?
integraldisplayintegraldisplay
approx
surface
φ
i
(x)G(x,x
prime
)
n
summationdisplay
j=1
α
j
?
j
(x
prime
)dS
prime
dS =0
Wewillgeneratedi?erentmethodsbychoosingthe φ
1
,...,φ
n
Collocation : φ
i
(x)=δ(x?x
t
i
) (pointmatching)
Galerkin Method : φ
i
(x)=?
i
(x) (basis=test)
Weighted Residual Method : φ
i
(x)=1if ?
i
(x) negationslash=0(averages)
Note 18
As noted in the equation on the top of the above slide, by substituting the
de?nition of the residual into the orthogonality equation given in the previous
slide, it is possible to generate n equations, one for each test function. The
generated equation has two integrals. The ?rst is a surface integral of the
product of the given potential with the test function. The second integral is
a double integral over the surface. The integrand of the double integral is
a product of the test function, the Green’s function, and the charge density
representation.
Asnotedintheslide,andaswewilldescribeinmoredetailbelow,threedi?erent
numerical techniques can be derived by altering the test functions.
8.3.6 Collocation
Slide 21
Collocation: φ
i
(x)=δ(x?x
t
i
) (pointmatching)
integraltext
δ(x?x
t
i
)R(x)dS=R(x
t
i
)=0 ?
summationtext
n
j=1
α
j
A
i,j
bracehtipdownleft bracehtipuprightbracehtipupleft bracehtipdownright
integraldisplay
approx
surface
G(x
t
i
,x
prime
)?
j
(x
prime
)dS
prime
=Ψ(x
t
i
)
?
?
?
?
?
A
1,1
··· ··· A
1,n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
A
n,1
··· ··· A
n,n
?
?
?
?
?
?
?
?
?
?
α
1
.
.
.
.
.
.
α
n
?
?
?
?
?
=
?
?
?
?
?
Ψ(x
t
1
)
.
.
.
.
.
.
Ψ(x
t
n
)
?
?
?
?
?
Note 19 Collocation
The collocationmethod, described in the aboveslide, uses shifted impulse func-
tions as test functions, φ
i
(x)=δ(x?x
i
). Impulse functions, as called “delta”
16
functions, havea sifting propertywhen integratedwith a smooth functionf(x),
integraldisplay
f(x)δ(x?x
i
)dx=f(x
i
).
Impulse functions are also referred to as generalized functions, and they are
speci?ed only by their behaviorwhen integratedwith a smoothfunction. In the
case of the impulse function, one can think of the function as being zero except
for a very narrow interval around x
i
, and then being so large in that narrow
interval that
integraltext
δ(x?x
i
)dx=1.
As the summation equation in the middle of the above slide indicates, testing
with impulse functions is equivalent to insisting that R(x
i
)=0,orinwords,
that the potential produced by the approximated charge density should match
the given potential at n test points. That the potentials match at the test
points gives rise to the method’s name, the point where the potential is exactly
matched is “co-located” with a set of test points.
The n×n matrix equation at the bottom of the above slide has as its right-
handsidethe potentialsatthe testpoints. The unknownsarethe basisfunction
weights. The j
th
matrix element for the i
th
row is the potential produced at
test pointx
i
by a charge density equal to basis function ?
j
.
8.3.7 Centroid Collocation for Piecewise Constant Bases
Slide 22
() ( )()
1
,
ii
n
tjtj
j approx
surface
xGxxdSα?
=
′′′Ψ=
g166 g179
2
x
n
x
1
l
2
l
n
l
Collocation point in
line center
1t
x
?
?
?
?
?
?
A
1,1
··· ··· A
1,n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
A
n,1
··· ··· A
n,n
?
?
?
?
?
?
?
?
?
?
?
?
α
1
.
.
.
.
.
.
α
n
?
?
?
?
?
?
=
?
?
?
?
?
?
Ψ(x
t
1
)
.
.
.
.
.
.
Ψ(x
t
1
)
?
?
?
?
?
?
Ψ(x
t
i
)=
summationtext
n
j=1
α
j
integraldisplay
linej
G(x
t
i
,x
prime
)dS
prime
bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipupright
A
i,j
Note 20
In the above slide, a speci?c collocation algorithm is described. First, the basis
being used is the constant charge density on n straight sections or lines, as de-
scribedabove. Notethatthereforethegeometryisbeingapproximated. Second,
the collocation points being selected are the centroids of the basis functions, in
this case just the center of each straight line. Note that the collocation point is
on the approximated geometry, not the original geometry. So, one can think of
the problem as having been restated to be on a polygon instead of the original
circle. One could also have selected the collocation points on the original circle,
but then the replacement interpretation does not hold.
17
In collocation, or point- matching, the charge densities on each of the straight
lines are selected so that the resulting potential at the line centers matches the
given potential. As the equations on this slide make clear, the matrix element
A
i,j
is the potential at the center of line i due to a unit charge density along
line j.
It should be noted that the matrix A is dense, the charge on line j contributes
to the potentialeverywhere. Also note that if linej is far awayfromlinei,then
A
i,j
≈length(line
j
)?G(x
t
i
,x
t
j
) (7)
triangleright Exercise 4 Suppose we are using piecewise constant centroid collocation
to solve a 2-D Laplace problem, so G(x,y,x
prime
,y
prime
)=log
radicalbig
(x?x
prime
)
2
+(y?y
prime
)
2
.
Roughlyhowfarapartdolinesectionsiandjhavetobefor(7)tobeaccurateto
within one percent? Assume linej has length of one. Does your answer depend
on the orientationof linej? Does your answer depend on the orientationof line
i? (You shouldansweryesto one ofthese and no to the other,do yousee why?)
8.3.8 Centroid Collocation Generates Nonsymmetric A
Slide 23
Ψ(x
t
i
)=
summationtext
n
j=1
α
j
A
i,j
bracehtipdownleft bracehtipuprightbracehtipupleft bracehtipdownright
integraldisplay
linej
G(x
t
i
,x
prime
)dS
prime
1
l
2
l
1
t
x
2
t
x
A
1,2
=
integraltext
line2
G(x
t
1
,x
prime
)dS
prime
negationslash=
integraltext
line1
G(x
t
2
,x
prime
)dS
prime
=A
2,1
Note 21
Consider the two line sections, l
1
and l
2
given in the above slide. For Laplace
problems, G(x,x
prime
)=G(x
prime
,x), and this suggests a symmetry in the underlying
integralequationwhichisnotrepresentedinthecollocationdiscretization. This
asymmetry is shownon the slide above by noting thatA
1,2
negationslash=A
2,1
.Thatis,the
potential at the center of l
2
due to a unit charge density on l
1
is not equal to
the potential at the center ofl
1
due to a unit charge onl
2
.
It is possible to scale the variables to improve the symmetry, consider a change
of variables
?α
i
=α
i
?length(line
j
).
18
In this changeof variables,the unknowns ?α
i
arenowthe net line chargesrather
than the line charge densities. In this new system,
?
A?α=Ψ,wheretheelements
of the matrix
?
A are given by
?
A
i,j
=
1
length(line
j
)
integraldisplay
line
j
G(x
t
i
,x
prime
)dS
prime
.
Under the change of variables, if line j is far away from line i,then
?
A
i,j
≈G(x
t
i
,x
t
j
) ≈
?
A
j,i
. (8)
In other words, the elements of
?
A corresponding to distant terms are approxi-
mately symmetric.
triangleright Exercise 5 Give an example which shows that the scaled entries of
?
Acan be
far from symmetric. Assume we are using piecewise constant straight sections
with centroid collocation and the 2-D Laplace’s equation Green’s function.
8.3.9 Galerkin
Slide 24
Galerkin: φ
i
(x)=?
i
(x) (test=basis)
( ) ( ) () () () ( ) ( )
1
,0
n
ii ij
japprox
surface
x R x dS x x dS x G x x x dSdS?? ?α?
=
′′′=Ψ? =
g166g179g179g179g179
() () ( ) () ( )
1
,
,
n
iji
japprox approx approx
surface surfacesurface
iji
xxdS Gxx xxdSdS
Ab
?α?
=
′′′Ψ=
g166g179g179g179
g9g12g12g12g10g12g12g12g11 g9g12g12g12g12g12g12g12g10g12g12g12g12g12g12g12g11
?
?
?
?
?
A
1,1
··· ··· A
1,n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
A
n,1
··· ··· A
n,n
?
?
?
?
?
?
?
?
?
?
α
1
.
.
.
.
.
.
α
n
?
?
?
?
?
=
?
?
?
?
?
b
1
.
.
.
.
.
.
b
n
?
?
?
?
?
If G(x,x
prime
)=G(x
prime
,x) then A
i,j
= A
j,i
? A is symmetric
Note 22 Galerkin method
In the slide above, we give the equations for the Galerkin method, in which the
test functions are equal to the basis functions. In particular, one generates n
equations for the basis function weights by insisting thatR(x) is orthogonal to
each of the basis functions. Enforcing orthogonality corresponds to setting
integraldisplay
?(x)R(x)dS =0
and substituting the de?nition of R(x) into the orthogonality condition yields
the equation in the center of the above slide.
Note that the Galerkin method yields a system of n equations, one for each
19
orthogonality condition, and n unknowns, one for each basis function weight.
Also, the system does not have the potential explicitly as the right hand side.
Instead, the i
th
right- hand side entry is the average of the product of the
potential and the i
th
basis function.
8.3.10 Galerkin for Piecewise Constant Bases
Slide 25
2
x
n
x
1
l
2
l
n
l
() ( )
1
,
,
iij
n
j
jline line line
i ij
x dS G x x dSdS
b A
α
=
′′Ψ=
g166g179g179g179
g9g12g10g12g11 g9g12g12g12g10g12g12g12g11
?
?
?
?
?
?
A
1,1
··· ··· A
1,n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
A
n,1
··· ··· A
n,n
?
?
?
?
?
?
?
?
?
?
?
?
α
1
.
.
.
.
.
.
α
n
?
?
?
?
?
?
=
?
?
?
?
?
?
b
1
.
.
.
.
.
.
b
n
?
?
?
?
?
?
Note 23
In the above slide, a speci?c Galerkin algorithm is described using the same
example basis introduced above, constant charge density onn straight sections
or lines. Once again, we will think of the problem as having been restated to
be on a polygon instead of the original circle.
In the Galerkin method, the charge densities on each of the straight lines are
selected so that the resulting line averaged potential matches the line averaged
given potential. As the equations on this slide make clear, the matrix element
A
i,j
is the average potential over line i, scaled by the length of line
i
, due to a
unit charge density along line j.
As with the collocation method, the matrix A is dense because the the charge
on line j contributes to the averaged potentials everywhere. Also note that if
line j is far away from line i,then
A
i,j
≈length(line
j
)?length(line
i
)?G(x
t
i
,x
t
j
) (9)
triangleright Exercise 6 Suppose we are using piecewise constant centroid collocation
to solve a 2-D Laplace problem, so G(x,y,x
prime
,y
prime
)=log
radicalbig
(x?x
prime
)
2
+(y?y
prime
)
2
.
Roughlyhowfarapartdolinesectionsiandjhavetobefor(9)tobeaccurateto
within one percent? Assume linej has length of one. Does your answer depend
on the orientationof linej? Does your answer depend on the orientationof line
i? (Youranswershouldbedi?erentthantheansweryougaveforthecollocation
method. Do you see why?)
20
9 Laplace’s Equation in 3D
9.1 Electrostatics Example
9.1.1 Dirichlet Problem
Slide 26
v
+
-
2
0 Outsi?Ψ=
is given on SΨ
potential
First kind integral equation for charge:
Ψ(x)
bracehtipupleftbracehtipdownrightbracehtipdownleftbracehtipupright
Potential
=
integraldisplay
surface
1
||x?x
prime
||
bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipupright
Green
prime
sfunction
σ(x
prime
)
bracehtipupleftbracehtipdownrightbracehtipdownleftbracehtipupright
Charge density
dS
prime
Note 24
Intheaboveslide,wegiveatypicalapplicationexample,determiningthecharge
density on the surface of conducting plates given an applied voltage. As noted
next to the parallel plate ?gure in the above slide, in the domain exterior to
the plates the potential, Ψ, satis?es Laplace’s equation, and the potential is
given on the surface of the two plates. In this particular example, the top plate
potential is Ψ=0.5V and the bottom plate potential is Ψ=?0.5V,whereV
is the voltage noted in the ?gure. As mentioned in previous lectures, it is also
important to note that
lim
|x|→∞
Ψ(x)=0.
That is, far away from the plates the potential decays to zero.
For this exterior Dirichlet problem, one can write an integral equation that
relates the surface chargedensity on the platesσ to the potential on the plates.
This integral equation, given at the bottom of the above slide, is often referred
to by physicists as the superposition integral. In the integral equation,xis any
point on the plate surfaces and the surface being integrated over is the union
of the top and bottom plate surfaces. Note that the integration surface is not a
connected domain, but this presents no di?culties.
9.2 Basis Function Approach
9.2.1 Piecewise Constant Basis
Slide 27
Integral Equation : Ψ(x)=
integraltext
surface
1
||x?x
prime
||
σ(x
prime
)dS
prime
21
()1 if is on
j
xx? =
() 0 otherwise
j
x? =
Discretize Surface into
Panels
Panel j
() ()
g44
1
Basis Functio
Represent
n
ii
i
xxσα?
=
≈
g166
Note 25
Consider solving the integral on the top of the above slide where the surface is
thesurfaceofthecube showninthelowerleft-handcorner. The ?rststep,aswe
have mentioned in previous lectures, is to develop a basis in which to represent
the surface charge densityσ.
The cube pictured in the slide has had its surface divided into panels, and
a basis is derived from the panels. In particular, one can associate a basis
function ?
j
with each panel j by assigning ?
j
(x) the value one when x is a
point on panel j,andsetting?
j
(x)=0otherwise. If σ is approximated by
a weighted combination of these basis functions, then the approximation is a
piecewise constant representation of the charge density on the surface of the
cube.
A few aspects of this basis set should be noted.
? The basis functions are orthogonal, that is ifinegationslash=j,
integraldisplay
?
j
(x)?
i
(x)dx=0.
? These basis functions are normalized with respect to l
∞
,notl
2
.Thatis,
bardbl?bardbl
∞
=1but
bardbl?
j
bardbl
2
2
=
integraldisplay
?
i
(x)?
i
(x)dx=panel area.
Finally, many papers in the literature on solving integral equations refer to
“panel methods”. The name is derived from the idea of breaking a surface into
?at panels.
In the application area of analyzing ocean wave forces on ship hulls, panel
methods are commonly used. However, it is not possible to represent a curved
hull with quadrilateral ?at panels. Researchers in the area often create a best
?t panelled surface in which there are gaps between the edges of the panels.
Such a discretization technique is often referred to as using “leaky panels”, a
very compelling image.
22
9.2.2 Centroid Collocation
Slide 28
Put collocation points at panel centroids
( )
1
,
1
i
i
n
cj
j
c
ij
panel j
xdS
xx
A
α
=
′Ψ=
′?
g166
g179
g9g12g12g12g10g12g12g12g11
()
()
11,1 1, 1
,1 ,
n
c
n
nnnn
c
x
AA
AA
x
α
α
g170g186Ψ
g170g186g170g186
g171g187
g171g187g171g187
g171g187
g171g187g171g187
=
g171g187
g171g187g171g187
g171g187
g171g187g171g187
g171g187
g171g187g171g187Ψ
g172g188g172g188
g172g188
g22g22
g23g25 g23 g23g23
g23g25g23 g23g23
g22g22
i
c
x
Collocation
point
Note 26
After one has decide on a basis with which to approximately represent the
surface charge density, the next step is to develop a system of equations from
whichtodeterminethebasisweights,denotedintheslideaboveastheα
i
’s. The
most commonly used approach to forming such a system is to use collocation,
though Galerkin methods are also quite widely used. Recall that in collocation,
the basis function weights are determined by insuring the the integral equation
is exactly satis?ed at a collection of “collocation” points. For panel methods,
the most common choice for the position of the collocation points are the panel
centroids, as shown in the cube diagram in the above slide.
The equation in the top right of the above slide relates the potential at collo-
cation pointx
c
i
to the weights for the panel-based basis functions. To see how
the equation was derived, consider evaluting the potential at thei
th
collocation
point using the original integral equation
Φ(x
c
i
)=
integraldisplay
surface
1
bardblx
c
i
?x
prime
bardbl
σ(x
prime
)dS
prime
, (10)
where Φ is the know potential on the problem surface and σ is the unknown
charge density. Substituting the approximate representation for σ,
σ(x) ≈
n
summationdisplay
j=1
α
j
?
j
(x)
into the integral equation results in
Φ(x
c
i
)=
integraldisplay
surface
G(x
c
i
,x
prime
)
n
summationdisplay
j=1
α
j
?
j
(x
prime
)dS
prime
, (11)
where G(x,x
prime
) ≡
1
bardblx?x
prime
bardbl
is used to simplify the formula. Exploiting the fact
that ?
j
(x)=1if x is on panel j, and zero otherwise, results in the formula at
the top right of the slide.
23
The system of equations from which to determine the basis function weights is
given in the lower right corner of the slide. The right hand side of the system
is the vector of known potentials at the collocation points. The i,jth element
of the matrix A is the potential produced at collocation point i due to a unit
charge density on panel j. The vector of α’s are the unknown panel charge
densities.
triangleright Exercise 7 Determine a scaling of the α’s (?α
i
= c
i
α
i
) such that the scaled
matrix
?
A has the property
A
i,j
≈
1
bardblx
c
i
?x
c
j
bardbl
when bardblx
c
i
?x
c
j
bardbl is much larger than a panel diameter.
9.2.3Calculating Matrix Elements
Slide 29
Panel j
i
c
x Collocation
point
,
1
i
ij
pa cnel j
xx
AdS′
′?
=
g179
,
ij
c centr
j
id
i
o
Panel Area
xx
A
?
≈
One point
quadrature
Approximation
x
y
z
t
4
,
1
in
0.25*
ij
co
ij
j
p
Ar a
xx
A
e
= ?
≈
g166
Four point
quadrature
Approximation
Note 27
In order to calculate the matrix entries for the system of equations described in
the previous slide, recallthatA
i,j
is the potential produced at collocationpoint
i due to a unit charge density on panel j. The formula for A
i,j
is given on the
top right side of the above slide.
The ?gure on the left of the above slide is a diagram of how one typically
computes the panel integral given on the top right. First, consider a shift and
rotatation of the coordinate system so that the panel lies in the x?y plane
at z =0, with the panel’s center at x =0, y =0. The ?gure in the top left
shows the panel in the shifted and rotated coordinate system. Note that the
collocation point must also be placed in the new coordinate system.
If panelj is reasonably well seperated from collocation pointi, it is possible to
approximate the integral given in the top right by a single point quadrature.
More speci?cally, one could approximates the integral of
1
bardblx
c
i
?x
prime
bardbl
by a product
of
1
bardblx
c
i
?x
centroid
j
bardbl
and the panel area. As show in the middle ?gure, a single
24
point quadrature is like treating the panel as if it were a point charge at the
panel’s centroid, where the point charge’s strength is equal to the panel area.
If the collocation point is close to the panel, then a single point quadrature
will be insu?ciently accurate. Instead, a more accurate four point quadrature
scheme would be to break the panel into four subpanels, and then treat each of
the subpanels as point charges at their respective centers. This simple idea is
shown in the ?gure at the bottom of the above slide. This four point scheme is
equivalent to
integraldisplay
panel
j
1
bardblx
c
i
?x
prime
bardbl
dS
prime
≈
4
summationdisplay
j=1
0.25?Area
bardblx
c
i
?x
point
j
bardbl
.
If the panel is a unit square in the x-y plane whose center is at the coordinate
system origin, then the four x
point
j
’s are (x,y,z)=(0.25,0.25,0), (x,y,z)=
(?0.25,0.25,0), (x,y,z)=(?0.25,?0.25,0),and(x,y,z)=(0.25,?0.25,0).
9.2.4 Calculating "Self Term"
Slide 30
Panel i
i
c
x Collocation
point
,
0
ii
ii
cc
Panel Area
A
xx?
≈
g9g12g10g12g11
One point
quadrature
Approximation
x
y
z
,
is an integrable sin
1
i
ii
paneli c
xx
AdS
′?
′=
g179
,
1
i
ii
pa cneli
xx
A dS′
′?
=
g179
Note 28
The diagonal termsA
i,i
can not be computed using the quadrature approxima-
tiongivenonthepreviousslide. Toseethis, considerthe?gureatthetopleftof
the above slide, where a panel has been shifted and rotated into the x-y plane,
and the collocation point is the center of the panel. The integral that must be
computed is given on the right side of the top of the above slide.
As shownin the middle of the slide, using a single point quadraturescheme will
fail, because the distance between the point charge approximation to the panel
and the collocation point will be zero. Therefore, the single point formula will
require computing the reciprocal of zero, which is in?nite. The problem is that
the integrand in
integraldisplay
panel
i
1
bardblx
c
i
?x
prime
bardbl
dS
prime
(12)
is singular. That is, the integrand approaches in?nity at a point x
prime
which is in
the domain of integration. What is not so obvious is that (12) is an integrable
25
singularity. Therefore, even though the integrand approaches in?nity at some
point, the “area under the curve” is ?nite.
9.2.5 Calculating "Self Term"
Tricks of the Trade
Slide 31
Panel i
i
c
x Collocation
point
,
1
i
ii
pa cneli
xx
AdS′
′?
=
g179
x
y
z
Disk of radius R
surrounding
collocation point
,
11
ii
cc
ii
disk restof panel
AdS dS
xx xx′′??
=+
g179g179
Disk Integral has
singularity but has
analytic formula
Integrate in two
pieces
2
00
1
2
1
i
R
dk cis
dS rdrd R
rxx
π
θπ
′
′==
?
g179g179g179
Note 29
In the above slide, we both show that
integraldisplay
panel
i
1
bardblx
c
i
?x
prime
bardbl
dS
prime
is integrable, and also give an idea about how to compute the integral.
As shown in the slide, ?rst rotate and shift the coordinate system so that the
panel is in the x-y plane at z =0, and so that the collocation point (or equiv-
alently the panel centroid) is at the origin. In this new coordinate system, the
integral can be written as
A
i,i
=
integraldisplay
panel(rs)
i
1
bardblx
prime
bardbl
dS
prime
where the notation panel(rs) is used to indicate that the integral is over the
rotated and shifted panel.
On the top left of the above slide, a circular disk of radius R and center at the
collocation point is inscribed in the rotated panel. In the equations that follow
the ?gure, it is noted that the panel integral can be recast as the sum of an
integral over the disk plus an integral over the rest of the panel. The integrand
in the integralover the rest of the panel is no longer singular, but the integrand
in the integral over the disk is still singular.
The integral over the disk can be computed analytically by using a change of
variables. After rotating and shifting the panel, the disk is in the x-y plane and
its center, equal to the collocation point, is at zero. Therefore,
integraldisplay
disk
1
bardblx
c
i
?x
prime
bardbl
dS
prime
=
integraldisplay
disk
1
bardblx
prime
bardbl
dS
prime
=
integraldisplay
R
0
integraldisplay
2π
0
1
r
rdrdθ
26
where r =
radicalbig
x
2
+y
2
, rsinθ = x and rcosθ = y. As given in the bottom of the
above slide, the integral over the disk is 2πR.
9.2.6 Calculating "Self Term"
Other Tricks of the Trade
Slide 32
Panel i
i
c
x Collocation
point
Integrand is singul
,
1
i
ii
paneli c
xx
AdS
′
′=
?
g179
g9g12g10g12g11
x
y
z
1. If panel is a ?at polygon, analytical formulas exist.
2. Curved panels can be handled with projection.
10 Summary
Slide 33
Integral Equation Methods
Exterior versus interior problems
Start with using point sources
Standard Solution Methods
Collocation Method
Galerkin Method
Integrals for 3D Problems
Singular Integrals
Wewillexaminecomputingintegralsnexttime,andthenexamineintegralequa-
tion convergence theory.
27