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