Numerical methods for pdes Integral Equation Methods, Lecture I Discretization of boundary Integral equation Notes by Suvranu De and J. White April 23, 2003
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 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 1st and 2nd 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 Standard solution methods Collocation Method Galerkin Method Some issues in 3D 3 Interior vs Exterior Problems SLIde 3 Interior Exterior yn on surtace known on surfac Temperature in a tank""Ice cube in a bath What is the heat Heat flow =Thermal conductivity sur face an Note 1 Why use integral equation methods? For both of the heat conduction examples in the above figure, the temperature, T, is a function of the spatial coordinate, z, and satisfies V2T(a)=0. In both
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 1st and 2nd 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 ∇ = T 0 inside 2 ∇ = T 0 outside Temperature known on surfac "Temperature in a tank" "Ice cube in a bath" What is the heat flow? Heat flow = Thermal conductivity surface ∂T ∂n Note 1 Why use integral equation methods? For both of the heat conduction examples in the above figure, the temperature, T , is a function of the spatial coordinate, x, and satisfies ∇2T (x)=0. In both 1
problems T(r)is given on the surface, defined by T, and therefore both problems are Dirichlet problems. For the"temperature in a tank"problem, the problem domain, Q2 is the interior of the cube, and for the "ice cube in a bath"problem, the problem domain is the infinitely extending region exterior to the cube. For such an exterior problem, one needs an additional boundary condition to specify what happens sufficiently far away from the cube. Typically, it is assumed there are no heat sources exterior to the cube and therefore limT(x)→0. For the cube problem, we might only be interested in the net heat flow from he surface. That How is given by an integral over the cube surface of the normal derivative of temperature, scaled by a thermal conductivity. It might eem inefficient to use the finite-element or finite-difference methods discussed in previous sections to solve this problem, as such methods will need to compute the temperature everywhere in Q. Indeed, it is possible to write an integral quation 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 n applications to have exterior problems where the known quantities and the quantities of interest are all on the surface 4 Examples 4.1 Computation of Capacitance SLIDE 4 potential given on s What is the capac Capacitance= Dielectric Permittivity J Example 1: Capacitance problem In the example in the slide, the yellow plates form a parallel-plate capacitor with an applied voltage V. In this 3-D electrostatics problem, the electrostatic potential y satisfies V-y(r)=0 in the region exterior to the plates, and the otential is known on the surface of the plates. In addition, far from the plates. 4-0. What is of interest is the capacitance, C, which satisfies
problems T (x) is given on the surface, defined by Γ, and therefore both problems 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 infinitely extending region exterior to the cube. For such an exterior problem, one needs an additional boundary condition to specify what happens sufficiently far away from the cube. Typically, it is assumed there are no heat sources exterior to the cube and therefore lim x→∞ T (x) → 0. For the cube problem, we might only be interested in the net heat flow from the surface. That flow is given by an integral over the cube surface of the normal derivative of temperature, scaled by a thermal conductivity. It might seem inefficient to use the finite-element or finite-difference methods discussed in 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. 4 Examples 4.1 Computation of Capacitance Slide 4 v + - 2 ∇Ψ= 0 Outsi Ψ is given on S potential What is the capacitance? Capacitance = Dielectric Permittivity ∂Ψ ∂n Note 2 Example 1: Capacitance problem In the example in the slide, the yellow plates form a parallel-plate capacitor with an applied voltage V . In this 3-D electrostatics problem, the electrostatic potential Ψ satisfies ∇2Ψ(x)=0 in 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 satisfies 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 permittivit 4.2 Drag Force in a Microresonator SLIDE 5 tom v⊥ew 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 vibrate using electrostatic forces. The changing character of those vibra tions can be used to sense rotation. The particulars of how the microresonato 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 fluid drag force on comb structures shown in the bottom part of theslide. The Huid is the air sur rounding the structure, and at the micron-scale of these devices, air satisfies the compressible Stokes equation Vu(a)= Vp(a) (1) V·u(x)=0 the fuid velocity and p is the pressure By specifying the comb velocity, and then computing the surface pressure and the normal derivative of velocity one can determine the net drag force on the comb. Once again, this is a problem in which the known quantities and the quantities of interest are on the surface
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 vibrations 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 fluid drag force on comb structures shown in the bottom part of theslide. The fluid is the air surrounding the structure, and at the micron-scale of these devices, air satisfies the incompressible Stokes equation, ∇2u(x) = ∇p(x) (1) ∇ · u(x)=0 where u is the fluid velocity and p is the pressure.By specifying the comb velocity, and then computing the surface pressure and the normal derivative of velocity, one can determine the net drag force on the comb. Once again, this is a problem 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 different package leads interact. To determine the magnetic interaction between signal currents flowing on different wires, one must solve H(r)=J(a) V·J(x)=0 where h is the magnetic field and J is the signal current density. By specifying he current, and then computing the magnetic field 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 LIDE i Note 5 Example 4: Capacitance of microprocessor signal lines This last example in the above slide is a picture of the wiring on a microprocessor integrated circuit. A typical microprocessor has millions of wires, so we are only looking at a small piece of a processor. The critical problem in this example
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 determine the extent to which signals on different package leads interact. To determine the magnetic interaction between signal currents flowing on different wires, one must solve ∇2H(x) = J(x) (2) ∇ · J(x)=0 where H is the magnetic field and J is the signal current density. By specifying the current, and then computing the magnetic field 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 This last example in the above slide is a picture of the wiring on a microprocessor integrated circuit. A typical microprocessor has millions of wires, so we are only 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 IEMS 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 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 Only need a on the surface, but T is computed everywhere Must truncate the mesh, = T(oo)=0 becomes T(R)=0 Note 6 Heat conduction in 2D In this slide above. we consider a two dimensional exterior heat conduction problem in which the temperature is known on the edges, or surface, of a square Here, the quantity of interest might be the total heat flow out of the square The temperature T satisfies V2T(x)=0x∈9
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 T = ∞ 0 at But, must truncate t mesh Surface Only need ∂T ∂n on the surface, but T is computed everywhere. Must truncate the mesh, ⇒ T(∞)=0 becomes T(R)=0. Note 6 Heat conduction in 2D In this slide above, we consider a two dimensional exterior heat conduction problem in which the temperature is known on the edges, or surface, of a square. Here, the quantity of interest might be the total heat flow out of the square. The temperature T satisfies ∇2T (x)=0 x ∈ Ω (3) 5
T(x) giuen a∈r 1imx-∞T(x)=0 where Q is the infinite domain outside the square and r is the region formed by he edges of the square. Using finite-element or finite-difference methods to solve this problem requires introducing an additional approximation beyond discretization error. It is not possible to discretize all of Q, as it is infinite, and therefore the domain must be truncated with an artificial finite boundary. In the slide, the artificial boundary a large ellipse on which we assume the temperature is zero. Clearly, as the radius of the ellipse increases, the truncated problem more accurately represents the domain problem, but the number of unknowns in the discretization increases
T (x) given x ∈ Γ limx→∞T (x)=0 where Ω is the infinite domain outside the square and Γ is the region formed by the edges of the square. Using finite-element or finite-difference methods to solve this problem requires introducing an additional approximation beyond discretization error. It is not possible to discretize all of Ω, as it is infinite, and therefore the domain must be truncated with an artificial finite boundary. In the slide,the artificial boundary is a large ellipse on which we assume the temperature is zero. Clearly, as the radius of the ellipse increases, the truncated problem more accurately represents the domain problem, but the number of unknowns in the discretization increases. 6
7 Laplace’ s Equation 7.1 Green’ s Function SLIDE 10 血n2D If u=log(V(-Io)2+(y-yo)2 hen岩+=0(x,y)≠(xo,) In 3D If hen+分+=0(x,y,2)≠(x0,,20) Proof: Just differentiate and see! Note 7 Greens function for Laplaces 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 given by 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 differentiation reveals that satisfies the 2-D Laplace's equation everywhere except = To, y= yo and u(, y, 2) (x-x0)2+(y-y0)2+(2-20)2 satisfies the 3-D Laplace's equation everywhere except I=Io, y= yo, 2=20 These functions are sometimes referred to as Greens functions for Laplace's equation b Exercise 1 Show by direct differentiation that the functions in(4)and(5) satisfy V-u=0, in the appropriate dimension almost everywhere
7 Laplace’s Equation 7.1 Green’s Function Slide 10 In 2D If u = log (x − x0)2 + (y − y0)2 then ∂2u ∂x2 + ∂2u ∂y2 = 0 ∀ (x, y) = (x0, y0) In 3D If u = √ 1 (x−x0)2+(y−y0)2+(z−z0)2 then ∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 = 0 ∀ (x, y, z) = (x0, y0, z0) Proof: Just differentiate 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 given by 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 differentiation reveals that u(x, y) = log (x − x0)2 + (y − y0)2 (4) satisfies the 2-D Laplace’s equation everywhere except x = x0, y = y0 and u(x, y, z) = 1 (x − x0)2 + (y − y0)2 + (z − z0)2 (5) satisfies the 3-D Laplace’s equation everywhere except x = x0, y = y0, z = z0. These functions are sometimes referred to as Green’s functions for Laplace’s equation. Exercise 1 Show by direct differentiation that the functions in (4) and (5) satisfy ∇2u = 0, in the appropriate dimension almost everywhere. 7
8 Laplaces Equation in 2D 8.1 Simple idea SLIDE 11 u is given on surface out fdex-*)+(ry Does not match boundary conditions! Note 8 Simple idea for solving Laplace's equation Here is a simple idea for computing the solution of La s equation o u(r, y)=alog V(a-xo)2+(y-yo where ro, yo is a point inside the square. Clearly such a u will satisfy Vu=0 outside the square, but u may not match the boundary conditions. By adjusting a, it is possible to make sure to match the boundary condition at at least one b 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 a. I 8.1.1 More points SLIDE 12 u is given on surface a2u du 0 out (x,y)● (., y.) Letu=∑11lg(√(-)2+(-)=∑aG(x-r;y-) Pick the ais to match the boundary conditions
8 Laplace’s Equation in 2D 8.1 Simple idea Slide 11 2 2 2 2 + = 0 outs u u x y ∂ ∂ ∂ ∂ Surface ( ) x y 0 0 , 2 2 2 2 + = 0 outside u u x y ∂ ∂ ∂ ∂ Problem Solved u is given on surface ( ) ( )( ) 2 2 Let log u xx yy = − +− 0 0 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 (x − x0)2 + (y − y0)2 where x0, y0 is a point inside the square. Clearly such a u will satisfy ∇2u = 0 outside the square, but u may not match the boundary conditions. By adjusting α, it is possible to make sure to match the boundary condition at at least one point. 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 2 2 2 2 + = 0 out u u x y ∂ ∂ ∂ ∂ ( ) x y 1 1 , ( ) x2 2 ,y u is given on surface ( ) x y n n , Let u = n i=1 αi log (x − xi)2 + (y − yi)2 = n i=1 αiG(x − xi, y − yi) Pick the αi’s to match the boundary conditions! 8
ote 9 To construct a potential that satisfies Laplace's equation and matches the undary 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 Greens functions. Of course, we have to determine the weights on the n point charges and the weight on the i n charge is denoted hereby a 8.1.2 More points equations SLIDE 13 Source Strengths selec to give correct potent - testpoints a(x-x,￥-y G(x-x,x-Y Note 10 contd To determine a system of n equations for the n ais, consider selecting a set of test points, as shown in the slide above. Then, by superposition, for each test point ti, yt,, u(xr4,v,)=∑alog√(xt1-xo)2+(vt4-yo) Writing an equation like(6) for each test point yields the matrix equation or the slide The matrix A in the slide has some properties worth notin . A is dense, that is Ai, 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 ith test point is nearest the ith charge, then Ai. i will be larger than Ai i for all j Item 2 above seems to suggest that A is diagonally dominant, but this is not the case. Diagonal dominance requires that the absolute sum of the off-diagonal entries is smaller than the magnitude of the diagonal. The matrix above easil violates that condition
Note 9 ...contd To construct a potential that satisfies 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 the n point charges, and the weight on the i th charge is denoted hereby αi. 8.1.2 "More points equations" Slide 13 ( ) x y 1 1 , ( ) x y 2 2 , ( ) x y n n , ( ) 1 1 x y t t , Source Strengths selec to give correct potent testpoints. ( )( ) ( )( ) ( ) ( ) 1 1 1 1 11 1 1 1 1 1 , ,, , ,, n n n n nn t t t nt n t t n t t t nt n t t Gx x y y Gx x y y x y Gx x y y Gx x y y x y α α −− − − Ψ = −− − − Ψ Note 10 ...contd To determine a system of n equations for the n αi’s,consider selecting a set of n test points, as shown in the slide above. Then, by superposition, for each test point xti , yti , u(xti , yti ) = n i=1 αi log (xti − x0)2 + (yti − y0)2 = n i=1 αiG(xti − x0, yti − y0). (6) Writing an equation like (6) for each test point yields the matrix equation on the slide. The matrix A in the slide has some properties worth noting: • A is dense, that is Ai,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 the i th charge, thenAi,i will be larger than Ai,j for all j. Item 2 above seems to suggest that A is diagonally dominant, but this is not the case. Diagonal dominance requires that the absolute sum of the off-diagonal entries is smaller than the magnitude of the diagonal. The matrix above easily violates that condition. 9