v . .fih. .mmwm. .50."... z 4’, ,9 ‘ I ‘ N in... J .3539}. o '3'- . .‘fivio c . . $5.30.... t Inmmflvuflmuaar 1.. “1...?! u 2 lex‘ a. tau. - . RH... .. .11 .r .r .x. Z. . 1 J in! a . .t.(.. I A}. (3.5. zhwaufivduqu 2:...» Mar—i 1‘ 2&6... . «mild»... WI. .Iurfi-miii... 3.1.] : "at. u v d? it: 21.: C. \ 4.43.1. s. a: it... I! ”H 3.4. ‘1... Eh». I , \ 3...; a m i. 33’ < V 1‘ “art... Jr \r - .; 1" 5?. kmusuMmm é. Hun . «dine... .I‘x‘nr: u «Wang .8 \ 5 ~‘llivi . rfi. km F1 . .i .. 1;. :u I a . .I g: .w,€..,... ‘ 8w. .‘. an: a: 1'13 3.5» 3 .2 53?. 12:32:... .i i . nu: , I sc 4.. #55. .3 .v .. .d. a 53¢.“qu .1 I: n :19... 44 3...: :- .,. Q :: . t . z)“:- u .t. .1... «It. I: 12.... r51?) 9.“ v . I l a l’.' 1! i l. e... is), o .- ‘1 (3:, z .1: Sin 53.3.... 4 .| .2. A. . .1... v11 .339 xv n 1.3. 2.3,}; I. kins...- _. . .3 :3! ,‘r.oaumLxK .11).? . rh irah . .. ,v \.. . ..I.o5. 3007- This is to certify that the dissertation entitled MATCHED INTERFACE AND BOUNDARY (MIB) METHOD AND ITS APPLICATIONS TO IMPLICIT SOLVENT MODELING OF BIOMOLECULES presented by YONGCHENG ZHOU b has been accepted towards fulfillment of the requirements for the PhD degree in Applied Mathematics Major Professor’s Signature Twat/1 2ST], 2005, Date MSU is an Affirmative Action/Equal Opportunity Institution LIBRARY Michigan State University 4 _.-;a.-.---—»—--,-I-.--o—.- —- - -.-.-4--—n-w~'-.--.-.-.—.-.--.-.-.-.-.- - PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 2/05 pLIClRC/DateDuejnddpJ MATCHED INTERFACE AND BOUNDARY(MIB) METHOD AND ITS APPLICATIONS TO IMPLICIT SOLVENT MODELING OF BIOMOLECULES By YONGCHENG ZHOU A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics 2006 © Copyright by Yongcheng Zhou 2006 All Rights Reserved ABSTRACT MATCHED INTERFACE AND BOUNDARY(MIB) METHOD AND ITS APPLICATIONS TO IMPLICIT MODELING OF BIOMOLECULES By Yongcheng Zhou This dissertation describes the matched interface and boundary (MIB) method [94, 96] for elliptic interface problems, in which the equations have discontinuous coefficients and possible singular source on the interface. This low regularity leads to the slow convergence or divergence of most traditional numerical methods for smooth problems. The complexity of the interface makes it more difficult to develop efficient numerical methods. The matched interface and boundary method is closely related to the methods with ghost points while it differs from these methods in the implicit enforcing of the interface conditions. A uniform Cartesian grid is used in the formulation of the MIB to take advantages of the conventional high order central finite difference schemes. Attentions are only paid to the irregular grid points near the interface where the difference schemes involve the fictitious values instead of solely solution values. A local coordinate transformation is used to project the interface conditions defined in the gradient direction on the interface to the coordinate directions, which allows the MIB method to handle irregular interfaces. By iterative application of the interface conditions, fourth and sixth order numerical methods for general elliptic interface problems are successfully derived for the first time. Both immersed interface method (HM) and MIB methods are used for the accu- rate solution of the Poisson-Boltzmann equation for electrostatic potentials. Molec- ular surface generated with MSMS is chosen as the dielectric interface and special techniques are developed to implement this triangulized molecular surface into the finite difference method with Cartesian grid. This is the first application of analyti- cal molecular surface in the finite—difference-based Poisson-Boltzmann solvers, and is also the first method which conserves the continutity Of potential flux at the molecu- lar surface. Substantial improvements in the surface potentials and the electrostatic solvation energies are found as a result of this work. It is also found that, by mesh re- finements, the solutions of all other traditional Poisson-Boltzmann solvers essentially converge to results of interface methods attained with a coarse mesh. Thus interface methods are relatively more efficient in achieving solutions of the same accuracy. The convergent solvation energies calculated from interface methods at coarse meshes not only permits an accurate energetic analysis of the molecular interactions but also pro- vide an accurate calibration for the Generalized Born method. The accurate potential itself also has significant implications in the analysis of protein/ligand binding and association. Acknowledgements This dissertation is the conclusion of what has been an extremely long journey, one that I could not have completed alone. First of all, I would like to thank my thesis advisor, Dr. Guowei Wei, for giving me the opportunity to come to MSU, and for his invaluable assistance, guidance and encouragement throughout my Ph.D study and various research projects. He led me into the fascinating world of the scientific research, gave me the opportunity to take on such a stimulating challenge, and has shown me many things that have aided in my scientific growth. Our innumerable times discussing and questioning the problems, and cheering any achievements together will leave a lasting impression. It is always an enjoyable experience working with him. I am also very grateful to Dr. Michael Feig for his help in my understanding of molecular biology, the applications of my mathematics in computational biology, and for providing me the computational resources. I also much appreciated my advisory committee members, Dr. Keith S. Promislow, Dr. Tien—Yien Li and Dr. Beisheng Yan for their directions, comments and suggestions. Sincere thanks are also extended to Dr. Shan Zhao for many productive discus- sions, and Dr. Seiichiro Tanizaki for helping me on CHARMM, and Ms. Sining Yu and Mr. Weihua Geng, for their critical reading of my papers as well as source codes, and Mr. Matthew Ryan Otoole for proofreading this dissertation, and Dr. Zhengfeng Zhou, Director of Graduate Studies and Ms. Barbara Miller, Graduate Secretary in the Department of Mathematics, for their generous help throughout my graduate study at Michigan State University. I would like to thank my family members back home in China, especially my ' parents, Xiaonv Ye and Chenhuo Zhou, and my sisters, Yonghong Zhou and Xuemei Zhou, for their long-time sacrifices, support, confidence and encouragement. Finally and most importantly, I would like to express my utmost appreciation to my wife, Yuhui Sun, for her constant love, understanding and support. vi Table of Contents Chapter 1 Introduction 1.1 Elliptic interface problem and its numerical methods: a brief review . . . . 1.2 Illustration of MIB method: 1-D example ................... Chapter 2 MIB Method: 2-D Formulation and Numerical Experiments 2.1 Second-order 2-D MIB scheme for irregular interface ............. 2.2 Higher order 2-D MIB scheme for irregular interfaces ............ 2.3 Subtleness in domain extensions ........................ 2.4 Interpolation formulation ............................ 2.4.1 One dimensional formalism ....................... 2.4.2 Two dimensional formalism ...................... 2.4.3 Comparison of two formulations .................... 2.5 Numerical experiments on irregular interfaces ................. Chapter 3 MIB Method: 3-D Formulation and Numerical Experiments 3.1 Second-Order 3—D MIB Scheme for Irregular Interfaces ........... 3.2 Numerical Experiments ............................. 3.3 Comparison of MIB and IIM .......................... Chapter 4 Poisson-Boltzmann Equation and Interface Methods 4.1 Implicit Solvent Modeling and the Poisson-Boltzmann Equation ...... 4.2 Analytical Solutions and Numerical Techniques for the Poisson-Boltzmann Equation: A Brief Review ............................ 4.3 Dielectric Interface ................................ Chapter 5 Application of Interface Methods to the Poisson-Boltzmann Equation 5.1 Implementation of molecular surface in MIB and IIM ............ 1 2 7 14 15 31 33 37 38 42 45 47 64 64 78 79 81 82 88 91 94 95 5.2 Numerical validation of interface method for the Poisson-Boltzmann equation101 5.3 Application of interface method for biomolecules ............... Chapter 6 Thesis Contribution and Eiture Research Plan 6.1 Thesis Contribution ............................... 6.1.1 The Matched Interface and Boundary Method ............ 6.1.2 Interface Methods for the Poisson-Boltzmann equation ....... 6.2 Future research plan ............................... 6.2.1 Theoretical analysis and fast algebraic solvers for the MIB method 6.2.2 Investigation of biological significance of accurate surface potential APPENDIX A Kirkwood Solution to Non-Centered Charged Sphere References vii 106 119 119 119 120 121 121 122 125 128 List of Tables 1.1 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 3.1 3.2 4.1 5.1 5.2 5.3 Numerical convergence test for the equation (18113)]; = f (2:) The computational domain is [0,2] and the interface is located at x = a = 0.74174243050441595. 13 = l,11.(.7:) = Sin(3:13) if .17 g a and fl = 10, u(x) = cos(3a:) if a: > a. ............ , ......... 12 Numerical efficiency test of the 2-D Poisson equation in Example 2.5.1. b = 10, [u] = 0, [flan] = —0.75. ..................... 48 Numerical efficiency test of the 2-D Laplace equation in Example 2.5.3. B = 1. ................................... 51 Numerical convergence test and accuracy test for Example 2.5.4. . . . 53 Robustness test of high order MIB scheme with Example 2.5.4. ,8‘ = 1. n; = 72y: 40. .......................... _ ..... 53 Numerical convergence test and accuracy test for Example 2.5.5. . . . 60 Numerical convergence test for Example 4. ............... 62 Numerical convergence test for Example 5. ............... 62 Comparison of the second-, fourth— and sixth-order MIB schemes for the Poisson equation (Case 6). The numbers in the parentheses are the CPU times in seconds for the generation of local difference schemes. 63 Numerical convergence test for Example(3.2.1) ............. 78 Numerical convergence test for example(3.2.2) ............. 79 Major Poisson-Boltzmann solvers. Notation: Poisson-Boltzmann equa- tion(PBE), multigrid(MG), algebraic multigrid(AMG) finite difference(FD), finite element(FE), Gauss-Seidel(GS), conjugate gradient(CG), succes- sive over relaxation(SOR), stochastic dynamics(SD), quantum mechan- ics(QM), molecular mechanics(MM), molecular dynamics(MD). . . . 93 Positions and radii of atoms of cyclohexane. .............. 100 Numerical convergence test for Example (5.1.1). ............ 100 Electrostatic solvation energy AG in kcal/mol and the error in the surface potential for a sphere with a centered unit charge. The exact solvation energy is -81.98 kcal/mol .................... 110 viii 5.4 5.5 5.6 5.7 5.8 5.9 5.10 Electrostatic solvation energy AG in kcal/mol and error in the surface potential for a sphere with a single noncentered charge. ....... Electrostatic solvation energy AG in kcal/mol and error in the surface potential for a sphere with a single noncentered charge. ....... Electrostatic solvation energy AG in kcal/mol and the error in the surface potential for a sphere with two noncentered charges. ..... Electrostatic solvation energy AG in kcal/mol and the error in the surface potential for a sphere with two noncentered charges. ..... Electrostatic solvation energy in kcal/mol of two positively charged spheres at four different mesh sizes h = 0.5, 0.2, 0.1 and 0.05. The distance between the centers of two spheres D vary from 2A to 5A. Summary of test proteins and the computed electrostatic solvation energies with MIB and PBEQ at different mesh sizes. ........ Summary of test proteins and the computed electrostatic solvation energies with MIB and PBEQ at different mesh sizes. ........ ix 111 112 114 115 116 List of Figures 1.1 2.1 2.2 2.3 2.4 Illustration of MIB method for 1-D problem. The dash line indicates the position of interface a: = a, where u(:r) has a discontinuity. Real solutions 22(2) at grid points are marked in black and the fictitious values f (:r) are marked in green. Fictitious values [HI and fi+2 are the smooth extension of the solution in left subdomain; fi_1 and f,- are the smooth extension of the solution in right subdomain. These fictitious values are solved by matching two interface conditions simultaneously. Three typical situations for an interface crossing the mesh lines. In the left chart, the interface passes the 33- mesh line at point A and y- mesh line at point B. The interface conditions are approximated at point A and point B to generate the difference scheme for (82255)]; and (Bug/)3” respectively. In the middle chart, the interface conditions are approximated at point A to obtain the difference scheme only for (flux); and the regular central difference is used for (Bug/)3, In the right chart, the difference scheme for ([3223])3, is set up from the approximation of the interface conditions at point B while term (1311;); is treated with the regular central difference. ...................... Irregular point (II, j) and the interface. The interface crosses the :1:- mesh line at (rmyo). The fictitious values are fw- and fi+1j (green). The vertical dash line is the auxiliary line on which three auxiliary points (in empty circle) are defined: (0,j + 2), (0,j+1) and (0,3) right at (1:0, yo). The jumps [72], [flun] and [”7] are evaluated at (1:0, yo). . . Irregular point (2', j) and the interface. The interface crosses the y- mesh line at (rmyo). The dash horizontal line is the auxiliary line on which three auxiliary points (in empty circle) are defined: (2 -— 2, o), (2' — 1,0) and (2', 0) right at (rmyo). The jumps [22], [flan] and [227] are evaluated at the intersect point (:ro, yo). .............. Irregular point (2', j) and the interface. The interface passes the grid point (2', j ) No auxiliary lines and auxiliary points are needed. The jumps [12], [min] and [227] are evaluated at the point (2, j) ........ Irregular point (2, j) and the interface. The interface crosses the :r- mesh line at ($0,310). The fictitious values are fm' and fi+1,j (green). The vertical dashed line is the auxiliary line on which three auxiliary points (in empty circle) are defined: (o,j - 1), (0,j + l) and (0, j) right at (1:0, yo). The jumps [22], [flan] and [227] are evaluated at ($0,110). . . 16 17 26 30 2.6 2.7 2.8 2.9 2.10 2.11 2.12 2.13 2.14 2.15 2.16 2.17 A typical stencil used in constructing a fourth-order scheme for um and um. There are two pairs of fictitious points in this case: fig, fi+l,j and f2—1,ja fi+2,j .............................. A typical stencil used in constructing a sixth order scheme for um and 223;. Now there are three pairs of fictitious points in this case: fi,j)fi+l,j and f2—1,jafi+2,j as well as fi—2,j,f2'+3,j ........... Fictitious value fig-+2 shall be used for the discretization of uy,uyy at grid point (2', j + 1). However, fi,j+2 can not be solved along the y-direction but it can be solved along x-direction, i.e., solved together With fictitious values fi—2,j+2:f2'-—l,j+2 and f2+l,j+2- ......... The distribution of fictitious points (stars) for a 40 x 40 mesh with irregular interface. Each fictitious domain roughly has two layers of fictitious points, supporting a fourth-order central difference scheme. The green stars represent the smooth continuation of the interior sub— domain and the red stars denote the smooth continuation of the outer subdomain. ................................ (Same as Figure 2.2) Irregular point (2', j ) and the interface. The in- terface crosses the 33- mesh line at ($0. yo). The vertical dashed line is the auxiliary line on which three auxiliary points (in empty circle) are defined: (0, j + 2), (0, j + 1) and (0, 3') right at (170,y0). The jumps [22], [Gun] and [227] are evaluated at ($0,310) ................ Computed solution (upper) and the error (lower) for the 2-D Poisson equation in Example 2.5.1. ,6 = 10, [u] = 0, [Bun] = -0.75. ...... The computed solution (upper) and the error (lower) for Example 2.5.2. Computed solution (upper) and the error (lower) for the 2-D Laplace equation in Example 2.5.3. 13 E 1 ..................... The computed solution (upper) and the error (lower) for Example 2.5.4 with second-order method and n3; = my: 40 ............... The computed solution (upper) and the error (lower) for Example2.5.5. Upper: the computed solution of Example(2.5.6); Lower: computed solution of Example(2.5.7). Both are computed with mesh 72,; .—. ny = 100. .................................... The solution of the high order MIB method for the 2-D Poisson equa- tion in Example 2.5.8. —u(:r, y) is plotted. ............... xi 31 32 34 36 43 49 50 62 3.1 3.2 5.1 5.2 Local topology around irregular point (2', j, k). The interface crosses the x-mesh line at the point S = (1:0, ya, 20) between (i, j, k) and (2+1, j, k), which are two irregular points, on which two fictitious values, fi,j,k and fi+1, j,k are defined (marked with green dots). Two auxiliary lines (dashed line) are sketched passing through 8', one on the x — y plane and the other on the x — 2 plane. Two auxiliary points (in empty circle) (0, j, k + 1) and (0, j, k + 2) are placed on the auxiliary line on the a: — 2 plane to facilitate the discretization of 22;; Also, two auxiliary points (0, j, 12+ 1) and (0, j, k +2) are placed on the auxiliary line on the :r — y plane to facilitate the discretization of 11;. ............... Local topology with interface intersecting at a grid point (2’, j, k). Two 73 fictitious values, fig-,1,7 and fi+1,j,k are defined (marked with green dots). 77 The illustration of mesh lines and triangle as well as their intersect points. If the intersect point happens to be a vertex, or passes through a side of the triangle, or the sum of three angles 01,02 and 03 is 27r, then the intersect point is inside the triangle ............... 2-D diagram of the accessibility determined by chasing intersec points of mesh lines and molecular surface. Accessible grid points are marked . by solid red circles and the inaccessible grid points are marked by empty black circles. For the upper mesh line (y = y4) which has 6 intersection points with the molecular surface(the right most one is counted twice as the mesh line is tangent to the molecular surface at this point), the grid points (1,4) through (3,4) are accessible because they are located before the first intersection point; grid points (4,4) through (5,4) are inaccessible because they are located between the first and the second intersection point; grid points (6,4) through (8,4) are accessible again because they are located between the second and the third intersection point; grid points (9,4) and (10,4) are inaccessible again because they are located between the third and the fourth intersection point; grid points (11,4) through (15,4) are accessible because they are located between the fourth and the fifth intersection point. Between the fifth and the sixth intersection points there are no grid points. After the sixth intersection point which is the last one on this mesh line, all the grid points are accessible. ........................ xii 97 5.3 The grid points and interface point needed for the formulation of the C51 Cf! CH IIM scheme [13] at the irregular point(2, j, k). The IIM method man- ages to find the Taylor expansions of the solution accurate to 0(h3) at both sides of the interface, for which a total of 16 expansion co- efficients are to be determined. Original IIM requires 10 grid points and 6 interface conditions at an interface point(p’ or p) to solve for the 16 polynomial coefficients. These 10 points(rcd) are chosen to be (2', j, k), other 6 grid points directly connecting to (2', j, k) and 3 more grid points inside the cube, preferably at the other side of the inter- face. IIM with optimization employs all the 27 grid points (red and green) inside the cube to produce constrained quadratic optimization problems from which 16 polynomial coefficients are solved. See [13] for the detailed formulation of 3-D IIM scheme. It is not necessary to choose the interface point to be the nearest one as p’. For simplicity we choose 7/, the intersect point of the interface with the mesh line, as the interface point for the irregular point (2', j, k). ........... Maximum absolute error 81 (upper left), Maximum relative error 62(upper right), average absolute errors(lower left) and average relative error 63 (lower right) in the surface potential of two positively charged spheres. The errors of MIB and PBEQ are plotted in black and blue, respec- tively. The average absolute error is the total of the absolute errors at all irregular points divided by the number of irregular points. The absolute error is in kcal/mol/e and the relative error is the percent- age. The :r-axis is the distance between centers of two spheres. Circle: 1:. = 0.5; Diamond: )2 = 0.2; Square: 12 = 0.1 ............... Left y-axis: the deviation of electrostatic solvation energy calculated with PBEQ at h. = 0.5 (circle), 12 = 0.25(diamond) and h = 0.2(square) from that calculated with MIB at h = 0.5. Right y-axis: the solid line plots the total absolute charge of the proteins. x-axis from left to right: 1ajj, 1vii, 2erl, 1bbl, lcbn, 2pde, lshl, lfca, 1ptq, lbor, luxc, lvjw, lfxd, 1hpt, liiibg, lbpi, 1r69, 451e, 1neq, 1a2s, lsvr, lfrd, 1a63 and 1a7m. ................................... The map of electrostatic potential on the surface of protein 1a63. Up- per: potential computed with MIB; Lower: Difference of potentials computed with PBEQ and with MIB. ................. xiii 107 Chapter 1 Introduction Many problems in fluid dynamics, material science and molecular biology can be modeled by elliptic equations with singular sources or discontinuous diffusion coeffi- cients at fixed or moving internal interfaces, such as viscous incompressible flow with an elastic weightless immersed boundary (Stokes problem of blood flow [63]), stable or unstable crystal growth (Stefan problem [22, 25]) or implicit solvent modeling in biomolecular simulation (Poisson-Boltzmann equation [28]). The singularities in the source and / or the diffusion coefficient will cause the reduction of the global regularity of the solution, and consequently, lead to poor performances of most standard nu- merical methods designed for smooth solutions, usually seen through degradation of the convergence rate or even leadint to divergence of the numerical scheme. This dissertation is about the design, analysis and applications of an accurate and robust numerical algorithm, the matched interface and boundary(MIB) method, for elliptic interface problems on Cartesian grids. The MIB method was first proposed for the numerical solution of the Maxwells equation for electromagnetic fields in media with straight material interfaces [94]. It will be demonstrated that the MIB method can be generalized to 2-D or 3—D elliptic problems with arbitrary interfaces given an appropriate local coordinate transformation, a sophisticated organization of interface conditions and a deliberated selection of grid points for the solution of fictitious values. Substantial improvement on accuracy and resolution. as well as savings of the computational cost are achieved with this new method. Among the variety of potential applications of interface methods, we are particularly interested in the application of MIB to the numerical solution of the Poisson-Boltzmann equation derived from an implicit solvent model of biomolecular salty, aqueous environment. This application possesses physical elegance, mathematical complexity and biological importance. 1.1 Elliptic interface problem and its numerical methods: a brief review An elliptic interface problem is given by the equation V- (,13Vu)— K(X)U(X) = [)(x), x 6 fl = (1+ U (2— (1.1) with an appropriate boundary condition on 80. For simplicity, Q is assumed to be a regular domain, such as a rectangle in two dimensions (2D) or a cuboid in three dimensions (3D). 82“ and 12+ are two complementary subdomains of Q which are separated by an internal interface P. On this interface either the diffusion coefficient 13 is discontinuous or the source term p(.r) is singular, or both have singularities. Supplemented by the underlying physics of the equation, we have two jump conditions at the interface: [it] = 'u+(X(s)) — u—(X(s)) = 95(3). [322"] = [3+227:(X(s)) — ,B—u;(X(s)) = 5(3). (1.2) where X(.s-) is a point on the interface F, s is the arc-length parametrization of the interface F, and n is the outer normal direction of the interface. The superscript :1: denotes the limiting value of a function from one side or the other of the interface. In either subdomain Q" or (2+ the function u(x) is smooth, but global regularity of the solution is usually low, due to the finite jumps in the solution and/or its gradient on the interface. Most of the standard numerical methods (excluding dis- continuous Galerkin or similar methods) assume. that the solution is smooth inside its discretization stencil, which will not hold for the case at hand if the discretization stencil includes grid points located in different subdomains. Special treatments or modifications are therefore necessary to extend the application of many well devel- oped numerical methods such as finite difference method and finite element method to these elliptic interface problems. One of the earliest investigations on the elliptic interface problems was conducted by Babuska in 1970 using the finite element method [3], whereby he formulated the interface problem as an equivalent minimization problem and incorporated all the jump conditions into the cost functions. This method and a technique developed by Xu [91], require the exact evaluation of integrals on the boundary and the internal interface, which is usually difficult to obtain for an arbitrary interfaces or boundaries. In the 1970-19808, Peskin [50, 63, 64] developed an immersed boundary method to model blood flow in the heart. In the immersed boundary method, the complicated time-varying geometric boundary is regarded as being immersed in the fluid, and the Navier-Stokes equation is solved on a simple Cartesian grid. The presence of the embedded boundary is modeled via a singular source localized on the interface. This singular source is then approximated and distributed over a small neighborhood near the interface, usually through some discrete delta functions. The immersed boundary method has had success in modeling flows in complicated time—dependent geometry [17, 92] due to its flexibility, efficiency, and robustness. Recently, Tornberg and Engquist have proposed and analyzed a class of regularization methods for treating the singular source terms [79]. An extensive review of the immersed boundary method for turbulent flow simulations can be found in [33]. However, the immersed boundary method is typically first order accurate in higher dimensions and the smearing of the jumps along the interface might prevent its application to circumstances where sharp jumps and/or the resolution of the solution at the interface are of concern. The success of these pioneering works and the identification of more interface prob- lems have motivated the emergence of a variety of interface methods in the last decade, among them a major advance in this field was due to LeVequc and Li [42]. T hcsc authors proposed the immersed interface method (IIM) for solving elliptic equations with discontinuous coefficients and singular sources. In the IIM, local corrections of finite difference schemes on irregular grid points where the discretized Laplacian operator involves nodes from both sides of the interface, are pursued throughout the domain. This is achieved by incorporating interface jump conditions into local Taylor expansions of the operator on irregular points, from which FD schemes with genuinely first-order accurate truncation error can be derived. The resulting scheme is of second order accuracy and preserves the jumps at the interface. For 2-D or 3—D problems, a local coordinate is typically required to offer a better representation of the jump conditions since they are given in the direction normal to the interface. The IIM is robust and efficient, and has been successfully applied to a variety of interface-related problems, such as moving interface problems [30], elliptic irregular domain problems [14] and 3—D interface problems [13], to name a few. The reader is referred to a re- cent review [45] or Li’s monograph [47] on IIM method for more details about these applications. Various extensions and further improvements of the HM have been considered in the literature [2, 6, 7, 14, 29—32, 35, 38, 41, 48, 49, 52, 69, 71, 81, 83, 90], includ- ing the formulations in polar coordinates [49], the formulation based on the finite element method [44] and the formulation proposed by Dumett and Keener [14] for solving anisotropic elliptic boundary value problems. On the other hand, the prob- lem of convergence and efficiency of the IIM has attracted many research interests because the original IIM [42] typically leads to a matrix of non—symmetric coeffi— cients though the original problem is self—adjoint and strictly elliptic. This reduces the number of standard fast solvers that can be utilized with IIM, and convergence may not be rigorously guaranteed [31]. To address this issue, a maximum principle preserving IIM was proposed by Li and Ito [46] to attain a diagonally dominant albeit still non-symmetric coefficient matrix. Specially designed multigrid solvers can then be employed to speed up the convergence of the maximum principle preserving 11M [1, 40]. For interface problems with pi(_'.cewise constant coefficients, a fast lIM was constructed [43] through introduction of an unknown jump condition for [22”]. To- gether with the elliptic equation, this unknown will also be solved numerically. The success of the fast IIM lies in the fact that it models the jump conditions in 22(2) and its gradient as correction terms for the standard finite difference approximation to the equation. The latter is solved with fast Poisson solvers whilst the correction term is updated using the solution of 22(1). A satisfactory convergence was found from this iterative procedure. Motivated by the fast IIM [43], an explicit jump IIM was developed by Wiegmann and Bube [90]. Recently, a decomposed IIM for elliptic equations with variable coefficients was proposed by Berthelsen [6]. The linear sys- tems in two IIMs [89] and [6] are all symmetric and diagonally dominant, allowing the use of conventional fast Poisson solvers. Another popular sharp interface scheme using the Cartesian grid is the ghost fluid method (GFM) originally developed for treating contact discontinuities in the inviscid Euler equations by Osher and his coworkers [18]; The GFM is typically first order accurate for interface problems, including the elliptic one [51] and could be of second order accuracy for elliptic irregular domain problems. In the flavor of the level set method which gives an implicit representation of the interface, the interface jump conditions are captured implicitly by extending values across the interface into a ghost fluid. On irregular grid points, when the FD discretized Laplacian refers to a node from the other side of the interface, a ghost fluid value instead of the real one will be supplied. Such a jump condition capturing procedure is directly incorporated into the numerical discretization in a way that the symmetry of the FD coefficient matrix is maintained, allowing the use of standard fast solvers. In higher dimensions, the jump in the normal derivative is correctly captured through a projection to Cartesian coordinate directions in the GFM, while the jump in the tangential derivatives is neglected [51]. In this way the GFM can be applied dimension by dimension. The GFM is very simple and robust, and its practical extensions to complex interface problems such as 3D moving interfaces or the multiphase Navier-Stokes equations are promising. Recently, an interesting jump condition capturing FD scheme was constructed by W'ang [84] by using a body-fitting curvilinear coordinate system. It is noted that there also exist numerous studies in the literature about the quite relevant elliptic irregular domain problem [36, 59, 60]. One way to solve this prob- lem is to embed the irregular domain to a larger regular computational domain, and reformulate the original boundary conditions as interface jump conditions. A simple Cartesian grid can then be adopted, so that various fast algebraic solvers developed in the literature can be utilized. This essentially converts an elliptic irregular domain problem into an elliptic interface problem. Due to this close relationship, the meth— ods originally designed for one type of problems may be extended for another one. However, in the immersed boundary problems, no solution is sought outside the do- main boundary, whereas, in immersed interface problems, interface jump conditions couple the solutions on both sides of the interface. Moreover, an immersed inter- face problem has two interface conditions by definition while an immersed boundary problem has only one ‘interface condition’ coming from the boundary condition of the original problem. In this dissertation, we will primarily focus on solving the elliptic immersed interface problems. The application of our method will be extended to a few immersed boundary problems. To solve an elliptic problem with an irregular domain or interface, a body-fitting grid can be employed [3, 10, 11]. Nevertheless, for certain geometrically complex domains, the construction of a good body-fitting mesh remains a nontrivial and time- consuming task, even though considerable progress has been made. Furthermore, a considerable increase of computational difficulties will be encountered for moving in- terface problems, where a moving mesh method is required to regenerate or deform the grid during the simulation. Therefore, numerous modified finite difference methods that are based only on a simple Cartesian grid have been chosen to solve the elliptic interface—related problems in the literature [1, 6, 17, 18, 30, 31, 33, 36, 40, 42, 44- 46, 50, 51, 59, 60, 63, 64, 84, 89, 92]. One obvious advantage of Cartesian grid methods is that there is no computational cost for grid generation. A Cartesian grid also allows the use of simple data structures and a standard FD stencil over a ma- jority of the domain. Moreover, many contemporary software packages, such as fast Poisson solvers, multigrid, level set method, etc. are mainly developed for a Cartesian grid, and thus could be taken advantage of. On the other hand, in order to properly maintain the accuracy at the interface, some extra numerical work needs to be done near the interface in a Cartesian FD method. In general the growth of computational cost due to these extra operations is very small compared to the solution of the linear system as the modification of numerical methods occurs only in a small region near the interface. The matched interface and boundary (MIB) method to be formulated in this dis- sertation is based on fictitious points, which are regarded as the smooth extension of the subdomains separated by an internal interface. The solution values at these fictitious points are determined by hierarchically matching two given interface con- ditions. The idea of fictitious points or ghost domain was also used in the literature for solving elliptic irregular domain problems by Mayo [59, 60] in the 19803. The ghost cell is introduced outside of the. (inner) domain as a fictitious domain. Similar ideas have been explored in Refs. [36, 92], but the fictitious values are smoothly extended or extrapolated according to the jump conditions. Our MIB formulation differs from these in the sense that fictitious values on the different subdomains are solved pairwisely and simultaneously. Most importantly, the determination of ficti- tious values are independent of the coordinate direction and the discretization of the elliptic equation. The resultant flexibility makes it possible for the MIB method to handle elliptic problems with irregular interfaces. 1.2 Illustration of MIB method: 1-D example First proposed by Zhao and Wei [94] for the solution of the Maxwell’s equation with straight material interface, the MIB method manages to implicitly enforce both inter- face conditions by modifying the discretization schemes near the dielectric interface. A major ingredient of the MIB method can be illustrated with the following l-D example. Consider a uniform mesh with grid spacing h. and let the interface be lo- cated at a: = a where at, g a S :r,+1 for some 2', as in Figure 1.1. Suppose the diffusion coeflicient is [3“ in the left subdomain and [3+ in the right subdomain. For a one-dimensional elliptic equation (1.3) with interface conditions (1.4) (Hi/1):; : p(.1:), (1.3) [u] = u+(a) - ”ll-((1), [Bum] = B+ur(a) — ,[3_uI(a), (1.4) a straightforward application of the finite difference scheme says at 1:,- 541 2712+1— 203'.“ 2+,5v:_1 zlui +54 2ui—1 l / I / h2 I / l / _—__ p($i), (1.5) where ,8i+1/2 and fl,_1/2 are the values of diffusion coefficient at positions 1,41” and 17,-_1/2, respectively. This discretization, based on Taylor expansions of smooth func- tions, renders a second—order convergence only if both u(x) and [3(1) are continuous, which is not the case for the elliptic interface problems. ”Ur—3 Uzi-2 l l l i-3 i—2 1—1 i ' i+1 i+2 i+3 i+4 lax Figure 1.1: Illustration of MIB method for 1-D problem. The dash line indicates the position of interface 1‘ = a, where 22(1) has a discontinuity. Real solutions 22(1) at grid points are marked in black and the fictitious values f (:13) are marked in green. Fictitious values f,“ and fi+2 are the smooth extension of the solution in left subdo— main; f,-_1 and f, are the smooth extension of the solution in right subdomain. These fictitious values are solved by matching two interface conditions simultaneously. The degradation of the convergence of standard finite difference on elliptic inter- face problems is attributed to failure to ensure the interface conditions. \Vith the MIB 8 method, one modifies the finite difference scheme near the interface by first smoothly extending the subdomains separated by the interface and then applying the standard finite difference scheme on the extended subdomains. As seen from Figure 1.1, a second-order central difference scheme at 1:, involves grid points :r,_1,:r,- and 1,-4.1, but n+1 is located on the other side of the interface. Since on the left subdomain, only the difference scheme at 3:,- involves the grid point on the other subdomain, it is sufficient to smoothly extend the left subdomain up to xi“ to accommodate the discretization at 12,-. Similarly, the right subdomain would be left-extended up to 13,-. The solution values on these two extended domains are called fictitious values f (:r) to distinguish them from the real solution values 22(x) at the same place. These two fictitious values f,- and fi+1 can be solved from an approximation of the interface conditions (1.4) + + + — — — (“10.26 + wo.2+1”2’+1 + w0,i+2“i+2) _ (w0,2—1“i—1 + wO,-iui + “"0,zt+1fi+1) = [22] (1.6) .++ + ,.+ .. —,— .. .—...— . [3 ((1)1in + ll,l,2’+luf+1 + ”"1,i+2“1+2) —- {3 (u21,i_1u.1_1 + 201,242,, + wl,i+lfl+1) = [22.2,], (1.7) where (waif,- + wfl,2‘+1“i+1 + wai+2ui+2) is the interpolation of 22+(a) on ui+1,u,¢+2 . | . . ‘ . '. . . . . + + and fictitious value I, mth corresponding interpolation weights wO,i+1’"’0,i+2 and 222+ (222+ f- + 222+ 22- + 222+ 22- ) is the a roximation of 22+(a) at the same 0,2: 1,2 2- 1,i+1 2+1 1,2+2 2+2 W r set of nodes via finite difference. Similar approximations to u"(3:) and u; (a) are conducted on the nodes 22,-1, 221 and fictitious value f,“ with respective interpolation weights or finite difference coefficients, which are calculated with Fornberg’s algorithm [21]. The equations (1.6 and 1.7 are linear in f, and f,+1: + — _ + + — — + + — — __ , , + + + [3 “11.4% _ B w1.2+1fi+1 — [13223;] - 3 (wl.2‘+lu‘i+1 + wl.2+2ui+2) +fi—(urii_lu,_1+luii22,). The fictitious values f,- and f,“ are solved to be the linear combinations of real solution values 22,-_1,u,;, 22,-+1,22.,-+2 and the prescribed jumps [22], [23223;]: ft = Cite—1 + C22": + C§W+1 + 03111242 + Cisl’ul + Célfiuxl, fm = c',+lu.,-_1 + 02“”, + c5212,-+1 + c:,+lu,-+2 + cg+1[u] + CEHUJux], where the vectors are “-57U1,2+1K1 + ""0,i+1K2 ci = D _ + + + CH1 _ (3 “’1,iK1+w0,,fiK2 D _ - — _ ,+ _ + K1 — (”’0,2-—1vwo,2v “Je.2+1~ “’0,2+2v120) _ —— —,— __,+,+ ___-++ K2 - (L3 “mi—1’13 “1,2, *3 “0,241, ’3 t”0.2+2’021) and D = —,23—w;j, +1203: 1 + [31126, i +1wf:i' These two fictitious values are supplied to the standard central difference schemes at .r, and 13+] such that 23,,_1/2u,-_1 —— 2(B,-_1/2 + [3,-+1/2)ui+.32+1/2f2+1 I22 {3+1 2ft — 2(2341 2 + 5+3 2W2“ + 54-3 2“i+2 I / I / h; / I / = 20(1‘2—1)» at Ti+l P012), at 3172, are actually conducted on smooth subdomains, and thus have a better convergence property than (1.5), as illustrated by Table 1.1. More fictitious values can be solved to support higher order discretization of the equation near the interface via iterative enforcing of the interface conditions as only 10 two fictitious values can be solved from two interface conditions. Following the above procedure, we first solve f,- and f,“ from 2+4 ' (u'01f2‘l' Z “'Okukk (: u’Okfl’k+w0i+1fi+l) = in] (1'8) 12: 2+1 k: 2-— 3 2+4 3+ 2011f, Z w :kuk) —/3 (Z w [221 +201 2+1f1+1) = [2322;] (1.9) k:2+1 k: 2— 3 which approximate the interface conditions(1.4) at fifth-order accuracy. We then solve fi—l and f2+2 from 2+4 (“’fii—Ifi-l + “6113+: 1" 3.2.122) - k:2+1 2' ( Z ”’(Tkuk + u'(f,z‘+1fi+1 + uiéi+2fi+2l = [H]. (1.10) k=2—3 2+4 23+(1Uifii_lfi_l + lvlfl'fi + Z uytkuk) _ k=2+1 2'. 213-( Z lL’ikltk+u’i—,i+lfi+l +u,1—.i+2fi+2) 2 “3111]. (111) k=2—3 Using two fictitious values at either side of the interface is consistant with a fourth- order central finite difference in the vicinity of the interface. for example, f2+1 4'12 “2—1 4111—2 “vi—3 . z __.+___ + , — 1 .--. 1.12 u” 12212 3212 4h? 32.2 12h2 a :2, 1' ( ) [2+2 4f2+1 22, 4712—1 "2—2 = ____+ _ __,____ 2 1.13 "$1 12h2 3112 4112 + 3212 12h2 a I’ ( ) Ui+3 4U2+2 "1+1 4fzi—1 fi—Q : _ + _. + — t ‘ , 1.14 "fl 12212 3h2 4h? 322.2 12212 a “73’“ ( ) “H4 411143 “2+2 4“2+1 f1 ,. . : ———‘- | — ‘ ‘ '— t . 7' , 1.1!, "ff 12212 321.2 4112 + 311.2 1211.2 a "+2 ( J) to achieve fourth-order convergence, as shown in Table 1.1. In addition to the convergence rate of the scheme, we are also interested in the relation between the accuracy of the MIB scheme and the magnitude of the jump in the diffusion coefficient /3. Considering a. secorid-order MIB scheme for which only 11 Table 1.1: Numerical convergence test for the equation (13221.); = f(.r). The compu- tational domain is [0, 2] and the interface is located at :2: = a = 0.74174243050441595. 23 = l,22(;r) = sin(3:r) if :2: S a and fl = 10, 22(1) = cos(31‘) if I > a. 211d MIB 4th MIB Loo Order Loo Order 20 2.80E—2 1.18E—4 40 6.89E—3 2.02 8.80E-5 3.75 80 1.64E—3 2.07 5.37E—6 4.03 160 3.98E-4 2.04 2.88E—7 4.22 320 1.02E—4 1.96 1.79E—8 4.01 one fictitious value is defined on either side of the interface, i.e., f, and f2“, the interface conditions ( 1.4) can be approximated by + + + (woyz‘fi + 2120,14]qu + u’o,i+2“2+2) - (zu&i_122,_1+ 211(ii22,+ 221(ii+1f,1+1) = [22], (1.16) + + + + ,13_(w1—,,-_1Uzi—1 + “Ti"“i + IL’;i+lfz‘+1) = [23221.]. (1.17) On the other hand, interface conditions (1.4) can also be approximated by using the exactly extended values 22;f and u- 2+1’ respectively defined at 1:, and at,“ + + + + , + 3 (WU-,i—lm‘l + 2163i“,- +u"’0—,2+1"2—+l + (30—123) 2 [22], (1.18) ++,+,.+ ,. ,+.. +2 {3 (“[1,1"1 + “'1,i+1"1+1 + ”’1,,j+2”2+2 +Cl h. ) [Bur], (1.19) — - — — — — 2 f3 (222191._122,j_1+w1122,+2vlfi+122li+l+ C1 12 ) where C7168” , 3'1— and (31+ are the coefficients of the leading truncation errors. Note that fictitious values f,- and f,+1 are the approximation to the exact values 22;L and 221.11 on fictitious domains, respectively. We denote the difference between these 12 fictitious values and the extended exact values as c,- and €i+1, i.e, 2.,- = at - f,- (1.20) 2 62+1 = U;+1‘f2+1- (121) By subtracting Eq. (1.16) from Eq. (1.18), and Eq. (1.17) from Eq. (1.19), we notice that e,- and 624.1 satisfy the following two equations: waif,- — mafia,“ = (00- - 05‘)h3, (1.22) B+wtici—fi_wii+lci+1 —_- (fi’Cf— (#03212. (1.23) Solving for e,- we obtain , _ ._ ,+ ’ ,_ 3 , — — — _ + + 2 _ —. 9+ _ + + . — . (3 “'0,2“’1,2+1+/3 “"1,2‘”0,2+1 (1.24) c,- 2 Since wit.“ ~ O(1),2ua:i ~ 0(1) ~ ()(1/12) and 221:2. ~ 0(1/h), it follows ’wl,2'+1 that 2,: ~ ()(113). Moreover, since — + - 3 ,— — +B+ 2 ‘(Co ‘ Co )wl,i+lh‘ “ “‘0,2+1(Ci ‘ C1 37)" , (1.25) c,- = + - 6+ ,+ — ‘1"0,2w1,2+1 + ,T—“Liwaz'H as 6+ [13" increases, 6,: asymptotically approaches C (123), where the constant C de- pends only on (11+ and 212;. In other words, for large jumps in the diffusion coefficient, the approximation error depends only on the solution of the problem. This proves that the present method is robust to a large contrast of ,{3 at the interface. 13 Chapter 2 MIB Method: 2-D Formulation and Numerical Experiments This chapter is devoted to the formulation and test of the MIB scheme for solving 2-D elliptic equations with irregular (curved) interfaces. Compared with the treatment of interface in 1-D or a straight interface in 2-D, it is much more intricate to treat irregu- lar interfaces because the local topology of irregular points varies from point to point. Therefore, the procedure outlined for regular interfaces cannot be directly applied. A second-order 2-D MIB scheme will be described first followed by its generalization to higher order schemes. We will then present an alternative interpolation formulation of the MIB method without the explicit use of fictitious domains. We show that the new formulation is essentially equivalent to the improved fictitious domain formulation. It is believed that this alternative formulation offers not only a better understanding of the MIB method, but also provides insights into other similar high order interface methods that might be constructed. Extensive numerical experiments are performed to demonstrate the high accuracy and robustness of the MIB method, as well as its flexibility in handling irregular interfaces. 14 2.1 Second-order 2-D MIB scheme for irregular in- terface Before proceeding to the construction of the MIB scheme for the discretization of (Bug); or (23223,)y, we have to identify the irregular grid points since only these points necessitate special care when the standard central difference scheme is applied to the whole domain. We define a grid point to be irregular if all the discretization points participating a standard central finite difference scheme at this point are not on the same side of the interface. For example, in a second order 2D scheme, an irregular point has at least one of its 4 nearest neighbor grid points lying on the other side of the interface. Note that the number of irregular points increases when a higher order FD scheme is employed. Assume as before that there are two given conditions associated with the interface, i.e., [22] 2 22+ — 22' : o(;r. y) (2.1) [Bun] 2 8+2”: — B—u; = 2/1(.r,y) (2.2) and also assume that both ¢(r,y) and tb(r,y) are C 1 continuous along the inter- face I‘. When considering the interface which is not always aligned with the 17- or y- mesh lines, as what shown by Figure 2.10, one more interface condition can be attained by differentiating Eq. (2.1) along the tangential direction of the interface, i.e. [227] = (6701,21). Hence for a point (20,210) on the interface, we have three jump conditions, [22] 2 22+ — 22.— 2 (5(20, yo) (2.3) [22¢] = 22: — 22; = 957(2'0, 210) (2.4) [137122] = .8+U+ — .8_'U'_ : ill/(1504101 (2-5) where the normal vector of the interface is 22' = (cos 6. sinfi) and 227 is the derivative 15 in the tangential direction 2' = (— sin 6. cos 6), while 0 S 0 < 27r is the angle between the positive .2:- direction and the vector 22’. Considering these relations, these three local interface conditions can be reformulated as [22.] 2 22+ — 22.— = cb(;170, yo) (2.6) [227] (—22;]" sin6 + 22; cos 6) — (—22; sin6 + 22,"; cos 6) = p(:2:0, yo) (2.7) [13227,] = fi+(22: c086 + 22;” si116) — 6’02; c050 + 22;; si116) = 212(20, yo). (2.8) In the MIB approach, the implementation of jump conditions is disassociated with (b) (C) Figure 2.1: Three typical situations for an interface crossing the mesh lines. In the left chart, the interface passes the :r- mesh line at point A and y- mesh line at point B. The interface conditions are approximated at point A and point B to generate the difference scheme for (,1'322JE)Jr and (Buy)y, respectively. In the middle chart, the interface conditions are approximated at point A to obtain the difference scheme only for (6221):; and the regular central difference is used for (fiuy)y. In the right chart, the difference scheme for (Buy)y is set up from the approximation of the interface conditions at point B while term (Burk. is treated with the regular central difference. the discretization of the elliptic equation. Also (23221.)1 and (822y)y will be treated separately. Therefore, we only need to illustrate how to locally recover the second order accuracy of the standard 3—point finite difference scheme for ((3221);. The mod- eling for (2322“,, can be achieved similarly. Considering an interface point 0 = (2:0, go) which is located at the intersection of interface I‘ and the :2: mesh line, see Figure 2.2, one fictitious point on each side of (IO, yo) is required in the present MIB modeling. +| ., To estimate these two fictitious values we will discretize 22 2122; and 22; involved in the jump conditions by using a 1-D grid partition as in the regular interface cases. However. it remains a problem to deal with the two derivatives 22,] and 22,7 in the 16 jump conditions (2.7) and (2.8) at the interface point (330,210). We overcome this difficulty essentially via two steps. First, by using jump conditions (2.7) and (2.8), we eliminate one 3; derivative that is more difficult to discretize near the interface. Second, we carry out the discretization for the remaining y derivative by using a one- sided FD scheme on the auxiliary points whose values are obtained by interpolation. This is usually possible in practice. If 22;] is easier to be evaluated, we will cancel 22,: <[(3.2%) F y $9.1“) % 3i (3701310 Figure 2.2: Irregular point (2, j) and the interface. The interface crosses the .2:- mesh line at (130,310). The fictitious values are ftj and f,+1,j (green). The vertical dash line is the auxiliary line on which three auxiliary points (in empty circle) are defined: (0,j + 2), (0,j +1) and (0,3) right at (1:0, ya). The jumps [22], [(3227,] and [227] are evaluated at (#150,110)- from (2.7) and (2.8) to attain [22] : 22+ — 22—, and [(3227,] — /3_ tan6[22T] : C1732}L — C1722; + CJUJ, (2.9) where C: = /3+ c056 + 23‘ tanfisinfl, C; = 13‘ cost? + ,r3_tanc93i119 and C; 2 (3+ sinH — 23" sinO. If otherwise 22,7 is easier to be calculated, one shall use the following jump conditions derived from Eqs. (2.6), (2.7) and (2.8) by canceling 22;: [2,2] : 22‘L — 22-, and [23227,] —,1’3+ ta116[227] = (7:2ij — C7722; - C" ” J y 22,]. (2.10) 17 where C; : [3+ cos6 + [3+ tan6sin 6, C; = ,6“ cos6 + [3+ tan6sin6 and C; = ,[3‘ sin6 — /3+ sin6. Cancellation of 225,]L or 22; from Eqs. (2.7) and (2.8) is required when modeling (23223))?” which leads to two other combinations of jump conditions. One combination IS [22] 2 22+ — 22—, and [13227,] + 6‘ cot 6[uT] = Cgu: + C322,; - CW2- y y, (2.11) where C; : (,6+—,t3‘) cos 6, C; = /3‘ cos6cot6+6+ sin6 and C; : 13_(cos6cot 6+ sin 6); and the other is [22] 2 22+ — 22— and [Hun] + 6+ cot 6[uT] = Cgu; + C322; — Cgug. (2.12) It is easy to check that here C; 2 (6+ — 23’) cos 6, C; = 13+ cos6cot6 + 6‘ sin6 and C; = [3+(cos6cot6 + sin 6). Either Eq. (2.11) or Eq. (2.12) can be chosen to formulate the discretization scheme for (62%),], depending on whether 22; or 22; is easier to be evaluated. In the ewe that either tan 6 or cot 6 is undefined in choosing jump conditions, i.e., whenever the outer normal direction is aligned to the x— or y- direction, the interface is locally perpendicular to the mesh line so that it can be treated as a straight interface locally. The procedure introduced in Section (1.2) for the 1-D problem can be directly employed to handle these irregular points. It is seen that for irregular interfaces each irregular point should have its own local MIB representation because of different local topology, while for regular interface cases, in contrast, only one global representation is needed throughout the whole domain. Some typical topology near irregular points is depicted in Figure 2.1. We first consider the detailed discretization of jump conditions (2.9) at the inter- face point (10,210) for the topology given in Figure 2.2, where the interface I‘ passes through point (:20, yo) between two irregular points (2', j) and (2 + 1, j ) Four nodes along the fit mesh line, i.e., (2 — 1.)), (2.)), (2 + l,j) and (2' + 2,3), are required 18 fat to approximate r-derivatives in jump conditions (2.9). In order to evaluate 22 J (20,210) = (:ro,;2/j), we add an auxiliary y- mesh line (dash line) through (20,210). Three auxiliary grid points on this auxiliary line will be employed to approximate 22;] . Two of these three points, (0, j + 1) and (0, j + 2), are on the positive side of the interface, while the third auxiliary point (0, j) is right on the interface. Refer again to Figure 2.2, where point (2 -~ 1, j ) and (2, j) are on the same side of the interface while (2 + 1,3), (2' + 2,3), (0,j + 1) and (0,j + 2) are on the other side. We deploy two fictitious points at (2, j) and (2 + l, j), the corresponding fictitious unknowns are denoted as fi,j and fi+1,j, respectively. The two conditions in Eq. (2.9) are then approximated to be ’ll’&i_l'fli_l’j + 2226;222:31- + “[0:2i+1f77+11]. + [22] : ”(ii/id +Iu'0fii+1"i+lvj + U’ai+2ll.i+2,j (2.13) .C; (201—,i_122,-__1,j + "T121113 + '2221-‘i+1f,j+1,j) + [6227,] - (3— tan 6[u,.] = Cff'll’iifid + 222:,+1U2+1,j + “IL-HUM”) v+ + +C’y (pf-Judj + Ptj+1paj+l + P1yj+2”o,j+2) (2-14) where superscripts, - and +, signify that the FD approximation is on the - and + side of the interface, respectively. Here 222 and p are FD weights for approximation along the 1?— and y- directions, respectively. Their first subscript (0 or 1) represents that it is for interpolation or for the first-order differentiation, while their second subscript is the node index. For example, 22.2(ii_1, 222612. and 222671.+1 are the interpolation weights of 22(ro,y0)- at grid points (2 — 1,j), (2,j) and (2+ 1,j). As aforementioned, 22,] shall be approximated on a one-sided finite difference stencil comprised of three auxiliary points, i.e., 22:~ -, 220.141 and nag-+2, which need J +. with function values at real grid points further numerical treatments. To relate 22 0 J we adopt the relation 113:]. : ”(:j + [11.] : u’fTJ—luf-IJ + u’fijuij + 4H,(I.i+lff+l~j + [It]. (215) 19 The other two auxiliary values, 220, j+1 and 220,342 are interpolated on three grid points around them to ensure an accuracy of 0(123). These interpolation points are chosen 22 1222022, and should be on the same side of the interface. For the geometry in Figure 2.2 we choose interpolation points to be (2,j + 1), (2 + 1,j + 1), (2' + 2,j + 1) and (2,j + 2), (2 + 1,j + 2), (2 + 2,j + 2) for 220,341 and 220,342, respectively. As a result of these interpolations, nag-+1 and 2200312 become known variables in jump condition (2.14). There are only two unknowns in Eqs. (2.13) and (2.14), i.e., fictitious values f2",j and [H.102 To solve for the representation of these fictitious values in terms of the real function values and the known jumps, we introduce two expansions f,- - = Ci - U, (2.16) 1] 2121.) = C111,, (2.17) where Cf = ( 1, .. ,C6) and CHI: (0,“, CH1, ,Cé“) are the expansion coefficients of the two fictitious unknowns with respect. to 6 real unknowns and 3 jumps which are also given by the vector U = (""2—1,j1“2,j1"‘2+1,j~""‘2+2.j1"’-0,j+11”0,j+21[7"l1l1dunl-lnTll- Inserting these expansions of ./2‘,j and f2+1,j to Eqs. (2.13) and (2.14), we end up with a small linear system for vectors Ci and CH1: Ci+1U—211+C"U Kl-U (2.18) || “0,2+1 (C; 111 —C;p:j221,1,+,)cf'+l1U—waflCi-U = Kg-U, (2.19) 12+l or equivalently 111,,,+,C'+1— 110,0 = K] (2.20) (C;11’i,-+1 - ijij222a,+l)C’+l —— CjwiiC' : K2. (2.21) 20 where _ — - + + _ K1 — (—u10,,_,,—u10’,.,u10,,+1,11)0’i+2,0,0, 1,0,0), _ + + K2 — ( ‘10 “’11—1+Cyp1,j“’01— 11 Crw1,1+Cypljw0,1’Cxw1,1+10110111121 + + + + — _ Cy p1j+11Cy p,,J-+2,Cy p14" —1,13 tanO). The solution of this 2 x 2 linear system is readily available with Cramer’s formula: _ "+.+ _ “’0,1+1K2‘7r((_ “’1,1+1 (’y1’1,j"’0,1+1)K1 Ci = (2.22) 1+, + ,._. _ ‘+ + _ _(t’.’L‘ 11,11,231! 0,i+l + U)O,i(C; ”Iii-Fl (1y P1,jw0,,+1) _ + + + Ci+1_ 03” w, iKl + 2110,,K2 (2 23) + + — ' ' ‘xC “’1,1“’0, 3 + “’0 1(C “’1 ,1+1 Cy 7’1,j“’0.i+1) Before finishing the solution of the fictitious values, we need to distribute the expansion weights of fw- and f,:+1,j at auxiliary points “OJ“ and 110,142 to their respective interpolation nodes defined above. Suppose '+1 'U'o.,j+1 1’ ° (111,142,U1+1,j+11u1+2,j+1) . '+2 . . . . ”0,14?- — I] '(”1.3+11“1+1,J+11U-1+2,J+1)1 where I].+1 = (IJ,+1,I%+1,I?,+1) and I] +2(IJ +2 I] +2 13+?) are the interpolation coefficients in the vector form. The expansion weights of fz',j and f1+l,j could then be distributed onto these interpolation points and the final expressions of fij and f1+1,j are: f1,j = (31111—14 + (3321111 + (3511111)“ + CZU1+2J 1 j+1 j+1 1+1 + 05(11 “11.141 + 12 U1+1.j+1 + 13 "1+2.j+1) 1 1+2 j+2 1+2 + 06(11 “11,142 +12 u1+1,j+2 +13 “11+2.j+2) + CHM] + C§[/111n]+ C§[2I,T] (2.24) 21 _ ‘+1 '.+1 '+1 ' 1 fi+l,j — CI] 11,-_1,j + C12 21“ + cs 21,11,3' '1" 02+ 111+2J 1+1 j+l j+1 j+l + C5 (11 u1‘,j+1 + 12 “i+l,j+l + 13 u1f+2,j+1) 1+1 2+2 2+2 2+2 + C6 (11 u1i,j+2 + 12 ui+l,j+2 + 13 U1+2,j+2) + C",+1 [11] + Cg+1[,1111,,] + Cg+1111,]. (2.25) One can therefore discretize (1811 11);: at irregular points (2', j) and (2 + 1, j ) as at the regular point: 2.“ 1 #1—11' [3;1j+/j1;1j 1+ j -—21 _ 1 a a . . (fiu1)1 = Tia—1,1 — 2 h2 2 “1,3‘ + —h22 f1+1,j at (21]) ,3+ 3+ + + + 1111.1 ‘ 115,1 5113.1 fii+§1 . . (1311101: = ——‘h2 (J — ’12 122-+1.,- + —h2 ui+2J at (2 +1,]) by substituting the above expansions to f1,j and [111, j. The known terms involving given jumps [u], [227] and [1312;] should be collected and combined into the right—hand vector B in the ultimate linear system Ax = B from which the elliptic equation is L. 1,-1.1 Dj+1 / (1:0, yo) (11211131,3 (111:6),9 ’ ' ' ’ :\ eventually solved. 1 1 ,1—1 Figure 2.3: Irregular point (2, j) and the interface. The interface crosses the y- mesh line at (2:0, yo). The dash horizontal line is the auxiliary line 011 which three auxiliary points (in empty circle) are defined: (2' - 2,0), (2 — 1,0) and (2,0) right at ($0,110). The jumps [11], [1311"] and [21.7] are. evaluated at the intersect point (170,310). The difference scheme for ([3 my) y can be generated f(.,1llowing the procedure similar 22 to that for (,Buxh. Let (2', j ) and (2, j + 1) be a pair of irregular points and between them the interface I‘ intersects the 2th mesh line at ($0,210), see Figure 2.3. We need two fictitious points, fi1j at point (2, j) and fig-+1 at (2', j + 1), to facilitate the discretization of (Bug/)3, at (2', j) and (2', j +1). Jump conditions, Eq. (2.12), are chosen and approximated at the intersect point (20, yo). Considering the fact that we do not have grid points to directly approximate 22; , a horizontal auxiliary mesh line passing through (2:0, go) is introduced and three auxiliary points, (2' —- 2, 0), (2' —- 1, 0) and (2', 0), are deployed on which a one-sided finite difference approximation for u; is formulated. The jump conditions (2.12) are then replaced by the following approximate equations: P55411114 + P(_):ju2',j +P6J+1f1yj+1+lul ’- pajfz'g' + Paj+1utj+1 +Paj+2utj+21 (226) C;(p1—‘j_lu,j,j_1 + pijui1j+ pij+lfi1j+1) + [Hun] + ,6+ cot 0(117] = Cy+(Pi*-in,j +Pitj+1u2§j+l +Pij+2u2,j+2) +C; (ml—,1-—2ui-210 + mil-421F130 + wiiugo). (2.27) The notations used here follow the same nomenclature as that in Eqs. (2.13) and (2.14). Since 21:0 is defined at an auxiliary point on the interface rather than at a grid point, we relate it with the unknowns at grid points by an interpolation: — + + + + “1,0 = “1,0 _ [11] : 1’0,ij + P0,j+1“j+1 + Pod-+210” — lU-l- (228) Nevertheless the interpolation on points 21,1, j_2, 11,, j and fi,j+2 is also applicable. With this relation, Eq. (2.27) is changed into - - — V — ,1 -+ , Cy (p13j_122,3j_1 + p1,j“2.j + p1,j+1fi1j+1)+ [13111,]«1—13 cot ()[ILT] ._ + + . . .+ . . + , . . — Cy (Pnjfw +P1.j+1“w+1 +1’1,j+2“z.2+2) +C;[u1,—,,_222,j_210 + 101,412,110 +'1111_‘,(])0fjf112j + ]);j+121,j,j+1 +7)fl,j+2"l'i1j+2 — [11D] (2.29) 23 Also, we are interested in the representations of the two fictitious values fiJ and f1, j+1 With respect to the known jumps and real unknowns, i.e., f2,j = Cj.U (2.30) f1,j+1 = Cj+l'U (231) Where U = (“2,j—11ui,j1Ui,j+l1u1',j+21”0,2—21uo,1’—11lul1lfiunl1lurl) and the compo- nents of vector Cj , Cj+1 are the expansion coefficients of f1, ,- and fig-+1 with respect to U, respectively. A 2 x 2 linear system for Ci and Cj+1 can be obtained by replacing fi,j and fin in Eqs. (2.26) and (2.29) with above representations, as follows 11,;,+,Ci+1-U—11,;jjcj.u = Kl-U (2.32) C inj+1cj+l ‘U ' (Cg-ptj + ngiipaj)cj ' U = K2 - U. (2.33) It is easy to verify that __ — — + + K1 ‘ (‘1’01—1"Po..1*P0.1+1’P0,1+21010"1-0’0) __—.— _-—— +.+ —,—+ K2 — ( Cyp1,j—1’ C1/p1,j’Cyp1,j+1+C$ l“1,1P0,j+11 + + — - + 1— — — - — + Cy p13j+2 + C1: u’l’ipo‘j+2, CI 101,?1—‘2‘ Cx w1._i_1, —C.’L' u’l’i, —1, —B C0126). After dropping the vector U from the above system, we obtain the solution of the system again by Cramer’s formula: Ci : -—(C,'fpil:j + (";;7"'1_,1P0-,j)K1 + 723:ng (2.34) CJPIJ'HPEJ' ’ p&,+,(Cthj + 0510172651) 0111 = ’05 21,11,119 + ”51+1K2 . (2.35) (’fy—IIZj+17)(l.j_p1T.j+l(331+p-1fij + (":17“’1—,17’(l,j) The expansions of C] and C1+1 solved from this linear system involve the un- knowns at auxiliary points 12091-2 and 120311, which are not regular points, so we also need to distribute the weights 011 these two points to the regular grid points. 24 Considering the local topology we choose grid points u.,~_1,j_1,u,-_1,j and Mum-+1 for the interpolation at auxiliary points 2204-4, and ""2'—2,j—l1 21,41]- and 22,12,141 for nay--2. With these two interpolation relations: '—1 ”—1 ”—1 'll.0,,'_1 : It1 '11.i._1,j_1+lz2 '11.,'_1,j+I§ 7"1—1,j+11 (2.36) '—2 '—2 ‘—2 u01—2 == 13 lH—23—1-fIE lH—zg-tlé 2t1-2,j+12 (237) it turns out that both fi1j and fig-+1 have the expansions in terms of unknowns at regular grid points and known jumps, as below _ j j ' ' f1.1‘ — Chum—1 + 02111.1 + 0311111“ + (331111.112 cj Ii—l . , Ii—l _ _ Ii—I , , + 5( 1 ul—I,j-I + 2 uz—l,] + 3 “2-1,J+1) j '—2 '—2 :_2 1' delj 1M—2J—1-tlé 1H—2J-tlé 1t1-2,j+1) + Cflu] + Cg,[,1311n] + C{,[uT] ' (2.38) '+1 '+1 "+1 “+1 f1.1‘+1 = Ci ”1.1—1+ 0% ”1.1 + Cf; ”1,111+ C31 “11,1412 '+1 7—1 '—1 7—1 + C"; (171 'Il.i_1,j_1 + I122 “'2—1,j + I; 'Il.i_1,j+1) '+1 '—2 :_2 '—2 + C‘é (Ill 11,-_2,j_1 + 122 22,;24- + If} 11,92,341) + C1,“ [21] + cg,“ 11:111..) + C§+1111T1 (2.39) With these two fictitious values, the discretization of uyy at point (23') and (2',j + 1) becomes straightforward (3.— 1 .1311 +137. 1 (3.11 . _ J—Q .. .12 J—Q... .72. , (Hugh, — (A—ypul’J—l — (Ag/)2 u’l-J + mfj+l at (l,}) 01.1 611-+fi11 811 (13“) : i2”. .__:?_l—_Z,,. .+_1:?_f. 1 at (1j+1) ” 3’ (A1)? “J (A112 1" (A1) ’+ ’ The known terms involving [11], [217] and [1311"] should also be collected and combined into the vector B in the final linear system AI = B. The finite difference schemes constructed above are applicable regardless whether the intersect point of the interface with the mesh line is at the grid point. However, significant simplification can be obtained if the intersect point is a grid point. When it comes to such a case, the difference scheme for (11113;)3 + ([312y)y can be generated in a single run rather than separately. lj+2 x 01‘” 1+1 i+2 22 Tj_1 Figure 2.4: Irregular point (2', j) and the interface. The interface passes the grid point (2', 3'). No auxiliary lines and auxiliary points are needed. The jumps [11], [Gun] and [UT] are evaluated at the point (2', j) Figure 2.4 shows such a situation where grid point (2, j) is an intersection point. Note that when [u] 75 0 the function on the interface is not well defined but from the computational point of view the interface itself can be treated either as in the interior of the interface or outside it. Here we regard it as in the interior of the interface, i.e. 11(1, j) = 21_(1?,_j) if the grid point (2', j) is on the interface. In such a case, one only needs to take care of the difference scheme at point (2', j) since both (2 — 1, j) and (2',j - 1) are now regular points. For (2 + 1,3) and (2,j + 1), the jump in [u] can be directly incorporated to generate the difference schemes 13: 3 , 13:3 113: , . 13: 1 , 2 ..J 1 ,j 1 .j 2 ,j ((311,), = figuflzr 2032:)? 2 1‘i+1*j+(T.12)7(2‘i:j+[1‘])(2'40) 26 at (2+ 1,3) and fi:j+§ 823% + 822% fi:j+% (”My”) 3 mum-+2 — (Ag/)2 Ui,j+1 + (ATP—(aid. + [12]) (2.41) at (2,,2-1-1). For point (2', j) on the interface, however, two fictitious values, fig-+1 and fi+1,j, are required in order to formulate the difference scheme for ([3121);- + (fill/51),, at the point (2', j). The expansion coefficients of these two fictitious values with respect to the unknowns at surrounding grid points and the known jumps would also be solved from the approximation equations of the jump conditions at point (2', j ) However, since the grid point is now on the interface, the interface conditions adopted are slightly different from what we used earlier. On the one hand, the Lagrange interpolation representation of the interface relation for 12, Eq. (2.6), does not involve either f,+1,j or f1, j+1 at all. Thus this jump condition can not provide an approximate equation for the fictitious values. On the other hand, two targeted +and fictitious values make the approximation possible for all derivatives, 22;, 22;,12,’ 21;. Therefore the cancellation of one of Eqs. (2.7) and (2.8) is no longer necessary and these two conditions exactly provide two approximate equations for f1.j+1 and [14.1,]; as follows —(1111—,,_111,_1‘j + 11.11—91.11.” + 1111—7,+1./',j+1,j)sin() + (P;j_1“2.j—l + Fiji-‘24 + P1_.J'+1f-1T,J‘+l)COS0 + [“7”] = —[w:,-(u,j,j +[11])+111:2-+,u,+1,j+ wii+2ui+2fl8i119 + [1%,]le + [11]) + pij+122w+1 + pij+221,-,j+2] cos 6, (2.42) {3-(111Ii_,11,_13j+ 211E.,21.i,j + u»’1-,,j+1f1f+l.j)(’056 +13_(12ij_112.,,j_1 + 121—Jun + pl—.j+lfi1j+1) sinl9 + [1311"] : 15+[wii(u,j,j + [11]) +111ii+111,+1 + '111:,+2u,-+2) c086 / + + ‘ , + 1 . . , + .1 . . .' / +13 ([21,).(12Lj + [11])+p,,j+,u,,1+1 + p].j+2”111+2) sin 1‘). (2.43) These two equations construct a 2 x 2 linear system for vectors C"+1 and CjH, 27 which are the respective expansion coefficients of fictitious values fi+1,j and fi,j+1a i.e. [24.10“ = Ci+l ' U, fihj+1 = Cj+1‘ U With U = (111—1,j,u2‘,j,lti+1,j=uz‘+2,j,U1,j—1,U1,j+1.uzf,j+2= lul~ [(31172], [MD- The solution of this linear system is _ — - 1+1 — j+1 _ fl’w‘ 0086Ci+1 +fi’ _- sinOCj+1 = K (2 45) 1,1’+1 p1,]+1 2’ ' where _. . _ + . . + _ K1 2 (tum;l 51116, (101,1. — w1,z.)sm6 + (pm- -— pm.) cos 0, + ,+ ,o _ — + + pi]- COS6’ -1111+J.si110,0,—1) —— — . + + - — . ,+ + — — - + + . . + + __ — — ~ 1+ + - + + . L3 w13i+1cos6hfi 11213242 c036, fl pl,j-l 81116, ,(3 191,].+1 sm6,,6 P14429110: 5+(wiicosfl +ptj sin6), —1,0), and U_pij+lsi116K1— pin cos 6K2 C"+1 (2.46) finwi,i+11’1,j+1 ”—2”1,-i+1 cos 9K1 + mm.“ smflKg Cj+1 (2.47) v‘i_“’1,1'+17’1,j+1 28 Thus we obtain the full representations of fi+1,j and f1; j+1 [my = C]+1u.,-_1,,- + C§+1um + C§+lui+1j+ C§+111i+2d~ + Ogle-J“ + cg+1u,;, M + Ci,“ [u] + Cg+1[uT] + cg‘,+1[pu,,] (2.48) fi,j+l = CJj+1Ui—l,j + CQHUM + Cg+lui+hj + Cit-“Urns + CgHUiJH + CéHum-Jrg + C;+I[11] + C?“ [217] + Cg+l[fiun] (2.49) and the difference scheme for (flux); + (,Buy)y at the point (i, j ) is 81,—+1, 5;1.+5,__1- [31:1 (fluxlx'l'wuyly = fifflld‘f I 2:122 2‘] z,j+ (A152): l—l,] 55% [323+] +flij—5 ‘ng— Wfi,j+la+ (Ay)2 Ui,j+ (Ag)? 1,1—1 (2-50) It is noted that the representation scheme described above has a local truncation error of 0(h3) so that the central finite difference scheme for the second-order deriva- tives have a local truncation error of 0(h) at irregular points. Because the number of irregular points is of 0(N) while the number of the total grid points is of 0(N2), the global approximation error of the proposed MIB method is dominated by the error of 0(h2) at regular points, which has been verified in our numerical experiments. Comments 2.1.1. The placement of anzrilarty points is flerrihle in the MIB method and it is not necessary to use a one-side difference scheme to approximate the partial derivative that is not in the primary direction. For example, considering the topology of the interface and the distribution of grid points in Figure 2. 5, we can place two auxiliary points on the opposite side of the sir-mesh line y = yj hence y: can be approximated with a central finite difference on these auxiliary points, in contrast to the situation in Figure 2.2 where both. aurihary points are placed on the same side of the mesh line. In this case, the expansion weights of fictitious values on ”OJ—l and Uo.j+l will be distributed to u(i,j - 1),11(i + 1,j — 1), 11(i + 2,j — 1) and u(i,j + 1). u(i + 1,] + 1), 21(i + 2,)" + 1), respecti-Uely. 29 Figure 2.5: Irregular point (i, j) and the interface. The interface crosses the x- mesh line at ($0,110). The fictitious values are fi,j and fi+l,j (green). The vertical dashed line is the auxiliary line on which three auxiliary points (in empty circle) are defined: (0, j - 1). (0, j + 1) and (0, j) right at ($0,310). The jumps [u], [flun] and [117] are evaluated at ($0, yo). Comments 2.1.2. In case that the interface crosses the mesh line at a grid point, it is possible to solve the fictitious values rather than to formulate the difference scheme at that grid point directly as did in Eq. (2.42) to Eq. (2.50). For three fictitious values fi‘j, fi+1,j and fig-+1 in Figure 2.4, fi’j can be directly represented as aw = [u] since it is on the grid point (i, j) which is regarded as inside of the interface by default. Therefore the fictitious value fi+1,j is be determined using the second relation in Eq.(2.9), i.e., [Hun] — ,13" tan()[uT] : (1:11; — (I; '11.; + ($713]. The approximation to this relation presents an equation for fi+1J as follows [Hun] — li— ta110[uT] — (If (pijmij + [11]) + pij+luid+1 + P:j+2'ui,j+2) = C';(?lli+:i('lti,j + [11]) +"-’i’:i+1"’i+l,j + 11:17:20.) -— — + + + Car (wl,i-lui-1J + wl,iuiaj + 'u’1.1+1fi+1,j)- (251) The third fictitious value f ,j, j+1 is solved from the approximation to the second relation 30 in Eq.(2.11) [dun] + fi_ COt’ 6[u7] — C: (wig: (”i-J + lull + wii+lui+l,j + wlff‘ljui'l'zj) 1‘ C'flPIJ-(Um + ["l) + PLHWJH +1113”) '- CJ(P:j_1U1,j—1+ ptjuij +pij+1fi,j+l)° (2.52) The central finite difference scheme at grid point (i,j) can then be formulated using these fictitious values. 2.2 Higher order 2-D MIB scheme for irregular interfaces '(0,j+4) t r o—é—o———o— 1(0J+3) & r I ‘+ a k c 53(0’102) .— [y KOJH) x o——+—e—5 + o— ——c c ' A a i—3 i—2 2 1+2 i+3 1+4 Figure 2.6: A typical stencil used in constructing a fourth-order scheme for um and ux. There are two pairs of fictitious points in this case: ff}, f1+l.j and fi_1‘j, fi+2‘j. The generalization of the proposed MIB method to higher order convergence is quite straightforward. In the present approach, a high order scheme means the use of standard high order finite difference discretization in the whole computational domain. In the vicinity of the interface, irregular grid points are first identified 31 according to the finite difference discretization. Fictitious points are then created at the same locations of the irregular points. The function values at fictitious points are determined by enforcing the interface jump conditions. For high order schemes, the same set of interface conditions is repeatedly used to determine the required set of fictitious values. At each irregular grid point, this procedure is exactly the same as what was described for the straight interfaces. (0J+6) f % G f t O—Q—Cl(o’j+5) .— c c Filo—’1'?) o— koi+3> G & koi+2> a—o t v—é—c c o— [y koi+li x +—+¢—o e t c 4 Figure 2.7: A typical stencil used in constructing a sixth order scheme for um and ur. Now there are three pairs of fictitious points in this case: fij, fi+19j and f1—1,j,fzt+2,j as well as fi—2,j= fi+3,j~ Let us consider a fourth-order case. We start with the interface condition (2.9) and the stencil shown in Figure 2.6 to derive a fourth-order scheme for um and u,r by using following steps: 1. Use vi_3‘j, - . . 311,33]- and fi+1,j as a stencil to approximate u- and 21;. For u+ and u; we choose the stencil f1, j. "i+l,j- - - - ,i1z-+4,j. Here, '11,] is discretized on . . . . . V aux1hary pomts (0,3). - - - . (0,] + 4), to ensure an accuracy of ()(h‘-’). 2. Solve the 2 x 2 linear system resulting from Step 1 for the representations of 32 fi,j and fi+1,j- 3. Use ui_3,j,--- ,uiJ and fi+1,j,f,:+2‘j to approximate 11‘ and a}, and use f,_1,j,f,',j and ”UH-1,93 - ~- aui+4,j to approximate u+ and a; Use the same approximation for 11;] as Step 1. Note that at this moment, both fio' and fi+1,j are known. 4. Solve the 2 x 2 linear system resulting from Step 3 for the representations of fi—1,j and fi+2,j- 5. Substitute appropriate terms of fi—1,j, f 1', j» f i+l,j, f i+2, j for the values at irreg— ular points when the standard fourth-order central discretization of am or an crosses the interface. This iterative procedure is systematic and is of arbitrarily high order provided there are enough number of regular points available on both sides of the interface. For example, a sixth order discretization of an and am can be constructed by using the stencil illustrated in Figure 2.7. In this case, we successively solve the representations for three pairs of fictitious values (figj,fi+1,j), (fi_1,j,fi+2,j) and (ft-2,jvfi+3,j)- It is noted that in constructing high order MIB methods, one needs to ensure that the local truncation error at the irregular point is one order lower than the designed global order of accuracy. The nature of high order accuracy of these schemes can be appreciated from the numerical examples presented in Section 2.5. 2.3 Subtleness in domain extensions We have not yet specified in sufficient detail how to determine fictitious values for a given topology of irregular points, which is a subtle issue indeed. At the first glance, fictitious values have to be solved along the direction of discretization, as they are needed in this discretization. In this regard, a fictitious value at a given point varies when it is obtained via smooth extension from different directions, because it is solved from the smooth continuation only in the direction associated with the original dis- cretization. According to this association between the domain extension and the 33 discretization, fi+1,j should only be solved along the x-direction, see Figure 2.8. Similarly, fig-4 should be determined twice, along the .r- and y-directions, respec- tively. Unfortunately, f1} j_1 cannot be solved along the y-direction up to fourth-order accuracy due to the large curvature in Figure 2.8. It is this difficulty that has severely constrained the applicability of the MIB scheme in Ref. [96], where only one generic elliptic interface problem with relatively low curvatures was solved up to fourth-order accuracy. \l 1 1+3 i+2 >———<5—— j+1 >———<>——j IR \l \r j-l j—2 — [ i-3 i-2 i-l T i+1 i+2 3:: Figure 2.8: Fictitious value fi’j+2 shall be used for the discretization of uy,uyy at grid point (i, j + 1). However, fig-+2 can not be solved along the y-direction but it can be solved along air-direction, i.e., solved together with fictitious values fr—2,j+2, fi—1,j+2 and fi+1,j+2- A close look at the fictitious values in Figure 2.9 reveals that the fictitious values must fall into one of the following five categories: 1. Those can be solved only along one of the :1:- and y-directions and will be used for the discretization along the same direction. 2. Those can be solved only along one of the x- and y-directions, but will be used for the discretization along the other direction. 3. Those can be solved only along one of the :r- and y-directions, and will be used for the discretization along both directions. 34 4. Those can be solved along both the x- and y-directions, and will be used for the discretization along either one, or both directions. 5. Those can not be solved along any direction, but will be used for the discretiza- tion along one or both directions. It is seen that a fictitious value in either the second, the third or the fifth category would be deemed as unsolvable. As seen in Figure 2.8, to find four fictitious val- ues fig, - - - , fig-+3 along the y-direction, grid points (:ri, yj_2) through (517,-, yj+1) are supposed to be on same side of the interface, whereas the other four grid points, (xi,yj+2) through (r,,yj+5), must be on the other side. Similarly, to find four fic- titious values fig-4, - - - , fig-+1, by considering the extension in the y-direction, we need points (23,-, yj_4) through (xi,yj_1) on one side and (:ri,yj) through (2:,, yj+3) on the other side. It is obvious that the distribution of grid points with respect to the interface in Figure 2.8 does not satisfy such an assumption. In particular, it is unable to find fictitious values fi,j_2 through fi9j+3 by the extension in the y-direction along the ith mesh line since there are only two grid points in one subdomain on that mesh line. According to the above discussion, however, at least four grid points in either subdomain are needed in order to define a fourth-order extension. These fictitious values can be classified into the third type, and were unsolvable by using the practice of Ref. [96]. To enhance the flexibility of MIB method in handling irregular interfaces, we intro- duce a new concept, the disassociation between domain extension and discretization. This concept is established based on the following error analysis. If the approximation errors of fictitious values at a given point (xi. yj) obtained from the z- and y-directions are both of 0(h”) for a positive n, the difference between these two fictitious values must be of 0(hn). Therefore, a fictitious value at an irregular point (16,-. yj), regard- less the direction in its calculation, can also be used for the discretization in any direction involving the grid point. For example, fictitious values fiJ through fi’j+3 are essentially the extension of the inner or the outer subdomain. Although they cannot be calculated through the extension in the y—direction, they indeed can be Figure 2.9: The distribution of fictitious points (stars) for a 40 x 40 mesh with irregular interface. Each fictitious domain roughly has two layers of fictitious points, supporting a fourth-order central difference scheme. The green stars represent the smooth continuation of the interior subdomain and the red stars denote the smooth continuation of the outer subdomain. obtained by the extension in the x-direction. Particularly, one can show that ft,j+2 can be calculated considering the extension along the (j + 2)th horizontal mesh line. The known fig-+2 can be used to facilitate the discretization, not only for um and an at (r,_2, yj+2) and (:ri_1,yj+2), resepctively, but also for uy and uyy at (:ri,yJ-) and (iii, y J41), respectively. This new understanding makes it possible to solve the ficti- tious values of the second and third types, and significantly broadens the applicability of the MIB method to general interface geometry. For the fictitious values of the fourth type, i.e., those that can be extended in both the :r- and y-directions, such as fi_11j_2, fi_1,j+3 in Figure 2.8, we may attain their values in the most convenient manner and use them for necessary discretization. Finally, fictitious values of the fifth type may be made available in most cases by refining the mesh. For example, if there are only four grid points inside a circular interface, fictitious values cannot be solved for all higher order schemes, i.e., orders higher than two, on such a grid. However, when the grid size is doubled in both directions, the fourth-order fictitious values can be obtained. For this reason, the present MIB method is a robust high order approach for all curved interfaces. 36 It was claimed in Ref. [96] that the proposed MIB scheme is of arbitrarily high or- der in principle. Indeed, the MIB procedure is systematic and it encompasses a variety of higher order schemes. However, the higher order convergence was demonstrated only for an interface problem on a special geometry due to the association between the domain extension and the discretization. It is believed that with the present understanding and extension technique, we are able to fully realize the higher-order potential of the MIB method for arbitrarily curved geometry. It is noted that for a general interface the difference scheme generated by the MIB involves more nodes than the HM and of course, more than the standard central dif- ference scheme. The number of grid points and the distribution of the MIB stencils vary with the underlying mesh and the local extension of the interface relative to the mesh. Also, the resulting linear system is no longer symmetric and diagonally dominant, as in the case of the HM. Nevertheless, all the iterative solvers we tested, including the successive over relaxation (SOR) or the preconditioned biconjugate gra- dient (PBCG), always yield the solution with a favorable convergence rate. Detailed comparison in terms of accuracy and CPU time is given in Section 2.5. 2.4 Interpolation formulation The MIB method described in the preceding sections makes use of fictitious domains and values for the standard high order finite difference discretization of the governing equation near the interface. The fictitious values give a smooth continuation of the solution across the interface. Interface jump conditions are iteratively used to deter- mine fietitious values. In this section, the possibility of the smooth continuation in terms of polynomial expansions is investigated. This interpolation formulation does not require the use of fictitious domains and values. To illustrate the idea, we start our discussion with a 1-D problem. The general principle is then applied to 2-D prob- lems with curved interfaces. Particular attentions are paid to the relation between these formulations. 37 2.4.1 One dimensional formalism Consider the 1-D elliptic problem given in Eq. (1.3). The interface is located at r,- g a S 23,-“. In a second-order central finite difference scheme, as,- is the irregular point to the left of the interface, and 1,-4.1 the irregular point to the right. On each side of the interface, we define a second-order polynomial '11—(17) = (16 + 97?” —— 1172') + $122.0: -— :ri)2, (2.53) + + “f a; 2 u (I) = 00 +T($—Ii+1)+h_2(x_$i+1) , (2-54) where the polynomial coefficients are scaled by h to reduce the condition number of the coefficient matrix. Obviously, a6 = u,- and a3" = u,“ and the remaining four coefficients can be determined from the expansion of the polynomials at two grid points, and from two interface jump conditions —al_ + a; = u,j_1 — u,- (2.55) a? + a; = 21,42 — Iti+1 (2.56) (-aTI—,~ + up?) — (0.1—11+ 02-112) 2 [u] — 11,41 +11,- (2.57) ,I3+(ail' — 2a.:arr) — 8"(a1- + 202—171) = h[,1’3uI] (2.58) where r, = a 1117i ,1‘,- = PILL—J. The first. two equations, (2.55) and (2.56), are the realizations of u‘ (at) at 13,1, and u+ (1:) at $142, respectively. The last two equations, (2.57) and (2.58), are the approximations of the interface jump conditions with u‘ (r) and 11+ (1"). Eqs.(2.55«2.58) essentially provide an algebraic system for the coefficients of u" (.r) and 11+(r), whose solutions are the representations of those four polynomial coefficients in terms of ‘11.,-_1, 11,-.iii+1.i1.,+2.e5 and 111. These. polynomial coefiicients 38 are solved by inverting the coefficient matrix of equation ( —1 1 0 0\ {01—\ fatal—u,- W 0 0 1 1 (12— : ui+2 -- 1124.1 , (2.59) —.r1 "112 —.r,» I; of" [u] — ”11,311 + u,- ( —,13- —2s-.r, 5+ —2B+.rr ) ( a; f ( hlfiurl f It can be seen from Eqs. (2.53) and (2.54) that h2' _ a; _ 2.,- + at + “a: (331') = T’ Ila-Afr) = 717.1% (33m) = 7;. and unfit”) = The l-D elliptic equation (1.3) can therefore be approximated at irregular points by ,(3— ’12 = p(.l‘,i) at 12', (2.60) 20+ {ff—1;;— = p(.ri+1) at :I:,j+1, (2.61) if fl is piecewise constant. Otherwise the approximation equations would be , a- 2a— firfirill—i'tdfxi) I13 = 0051') at In (262) 11+ 20+ .1317(I'1’+1)‘}_:‘ +/3(Ii+1)—h§— = p(Ii+1) at $i+1- (2-63) The constants in the representations of (1.1—.a2_,af and a; in the above equations, i.e., the terms involving given jumps ab and w, should be moved to the right hand side. This finishes the establishment of a second-order interpolation scheme at both irregular points. The above procedure can be easily generalized to construct higher order schemes. For example, to attain a fourth-order scheme at irregular points near the interface, one can start the formulation of the scheme by defining a fourth-order polynomial on 39 each side of the interface _ _ a a— a a— u (r) = a0 + -#(a: — xi) + h—22(r — :13)2 + h—%(:1: — 10,-)3 + "$4“ — Liz-)4, (2.64) + + “f “i 2 a; 3 “if 4 U (I) = 00 + 72—(31? — Ii+1) + [7(3 — n+1) + fit? — l‘i+1) + hjff — 33i+1l - (2.65) The first two coefficients a6 and a6" again have to be u,- and u,“ respectively. The remaining eight coefficients are to be determined by using two interface jump condi- tions, and by expanding these two polynomials at six grid points, i.e., points 13.], ri_2 and 22,--3 for Eq. (2.64), and points 13242, 35,-+3 and 17,-+4 for Eq. (2.65). We could therefore end up with eight linear algebraic equations ( —3 9 —27 81 0 0 0 Oi _2 4 -8 16 0 0 0 0 -1 1 —1 1 0 0 0 0 0 0 0 0 1 1 1 1 [ 0 0 0 0 2 4 8 16 0 0 0 0 3 9 27 81 —-.r1 —I12 ”I? ‘33? —x’ I? fit; $3 (43‘ ”fl-“ _3,,_,,.12 41,—3,23 (3+ 4pm,. 315+]; —4li+_rg f — — — — + + + +T (al,a2,a3,a4,al,a2,a3,a4) : (full-3 _ uii ui—2 — uii u'l—l — ui: “i+2 _ lll+13ui+3 — ui-i-l: ui+4 — ui-i-l: [u] — ui+l + “i: h'lfl'llxllT where the first three equations are generated from 11’ (it), while the next three equa- tions are generated from u+(:1:). The last two equations are the approximation of two interface jump conditions with polynomials u“(.r) and u+(r) at a: = a. By inverting the coefficient matrix, the representations of all the eight polynomial coefficients are solved, and expressions for u‘(.r) and u+(;r) are fully determined. 40 Unlike the second-order case, where each polynomial is used only once in for- mulating the difference scheme at either 1:,- or 272;“, for the fourth-order case, each polynomial is applied twice since there are two irregular points on each side of the interface, i.e., 23,3-1 and x,- on the left, and 23,41 and 1,42 on the right. These points are identified since the regular central difference scheme at these points involves the grid point(s) on the opposite side. In particular, u‘(a:) would be used to formulate the difference scheme at the points $17-1 and mi, while u+ (2:) is used at mi.” and 33,-+2. Assuming a piecewise continuous ,6, the difference schemes at the two closest irregular points are a— 2a7 aaafwaah—g = per.) at x. (2.66) a+ 2a+ .13x(1‘i+1)7:++fl(17i+1)722' P(Ii+1) at gr1+1- (2.67) At the other two irregular points, the difference schemes are more complicated since the derivatives of the polynomials involve more terms 10 (I'—1) _ _ _ _ BLT—1) _ _ — ifml _ 25,2 + 3a3 — 4a4 ) + $90.2 — 6a3 +12a4) = Mitt—1) at ilk—1, (2-68) B~ 8 . J,(:12+2) (a? + 20.3- + 3113+ + 40.1) + —’ (7:;2) (2a; + 60.;- + 120$) = Pf$i+2l at $i+2~ (2-69) It is concluded from the above discussion that at irregular points near the interface, a difference scheme of arbitrary order can be formulated by following these steps: 1. Define a polynomial on each side of the interface u_(;r) = 21,: + Z %(:r —— cm), (2.70) k=l l 71 0+ u+(.‘z:) -— 'Ui-l-l + 2 ET“? -- $14.1), (2.71) k=l ' where n, the order of the polyncm‘iial, depends on the global accuracy required. 41 These polynomials are essentially the approximations of the solution in the vicinity of the interface. 2. Expand each polynomial at the nearest n — 1 grid points to the interface on the same side. The total number of expansions is therefore 2n — 2. 3. Approximate two interface jump conditions with these polynomials. Differenti- ation is needed in approximating [flux], the jump in flux. 4. Compose a 2n x 2n linear algebraic system by combining 2n — 2 expansion poly- nomials and two discretized interface jump conditions. Solve this linear system by inversion, and the solution is the representation of polynomials’ coefficients in terms of the approximation solution at involved grid points and the given jumps. Approximate the original elliptic equation at all the irregular points by differ- C31 entiating the polynomials. Replace the polynomial coefficients with their rep- resentations. This gives rise to a finite difference scheme at the corresponding grid point. 2.4.2 Two dimensional formalism The new formulation of the MIB method for 2-D problems with curved interfaces can be accomplished similarly. The essential idea is to construct l-D interpolation polynomials, instead of 2-D ones. Referring to Figure 2.10, we have two irregular points, x,- and “+1, along the r-axis for a second-order scheme. Two l—D polynomials, in exactly the same form as Eqs. (2.53) and (2.54), are defined 11'(.r, yj) = a6 + -h—l(.1: — 21,) + [fake — 1,)2, (2.72) + + ”f “'2 2 u (1‘. yj) 2 a0 72—(1: - 17,41) + h—2(r — 17,-+1) . (2.73) These polynomials can be regarded as the 1-D approximation of the solution in the vicinit of the interface alon ,2 : i - mesh line. Here, one has a" = u- - and 1 JJ 0 I.) 42 Figure 2.10: (Same as Figure 2.2) Irregular point (i, j) and the interface. The interface crosses the 15- mesh line at (1:0, yo). The vertical dashed line is the auxiliary line on which three auxiliary points (in empty circle) are defined: (0, j + 2), (0, j + 1) and (0, j) right at (130,310). The jumps [u], [Bun] and [ur] are evaluated at (120,310). a; = u,+1,j by definition. To solve for the rest four coefficients, these two polynomials are to be evaluated at (17,--1411) and (13,42, yj), which gives two equations. The other two equations are obtained from the approximations of interface jump conditions (2.74) and (2.75) [96] I[u] = u+ — 21—, (2.74) [Bun] — L3— ta116[uT] 2 C311: + (71:11.; + CIT-11;, (2.75) where 6 is the angle between the normal vector and the :r—axis, C; = fi+ c086 + 5“ tanOsinB, C; = —,(3"(cos6 + tan65i116), and C; 2 (8+ — ,8_)sin 0. Although the approximation to Eq. (2.74) is trivial, special care has to be taken in treating Eq. (2.75), which couples two directions. Here we approximate Eq. (2.75) as C',+ ([v— —l‘::—(n.;L — 20.3.1.3.) + f—(rlf + 2072—11) + (2.76) CJUTJUJJ' +Pij+1llo.j+l + Pij+guaj+2l = [£31172] — (3— ta116[uT]. 43 . . . . . 1' ’ 13‘ 17' — I where p 18 the FD we1ght 1n the y-direction, and I, = —0T—5, :rr = —flh——O-. Note 2. that here auxiliary values ”(if nag-+1 and nag-+2 are also introduced to calculate u; as done in the fictitious domain formulation [96], where -+_ — _ — —2 uo,j — “10,310 + [u] — (”233' + (11 z, + a2 1:1) + [u]. (2.77) The final algebraic system for four polynomial coefficients is f —1 1 0 0) (q) 0 0 1 1 a; —$1 "$12 "Ir $3 ail. 1 — + 2 — + 1 2 + {MAJ-W4 ) WHJ—WHJ ¢’WHJ+WJ \ [fiun] — (3“ tanfilur] — C; {pa-(um- + [u]) + ptj+1uo,j+l + ptj+2u0J+2l ) Desirable difference schemes for ur or "u.” at two irregular points can be obtained from the direct differentiation of these two polynomials after all of their coefficients are determined. Following a similar procedure, one can determine polynomials u- (13,-, y) and u+(:ri,y), which are then used to approximate 11;, ugy, u; and 21.33, at corre- sponding irregular points. High order 2—D interpolation schemes can be established in a manner similar to what was described for 1-D high order schemes. However, for general interfaces with large curvatures, the construction of a high order interpolation MIB scheme is subject to the same difficulty as that in the fictitious domain formulation, i.e., there are not sufficient grid points to support high order polynomials along all the J:- and y-directions on one side of an interface, see at = 1:,- mesh line in Figure 2.8. It normally takes four grid points to determine a fifth-order polynomial u.— (.rl-, y), which is required for calculating u], and uyy at four-irregular points, .ri‘j_2, mid-4, 17-i,j+2 44 and mid-+3, in a fourth-order scheme. One solution to this problem is to make use of two interface intersection points, which would provide four jump conditions. However, a simple remedy is to make use of two other polynomials found in the .r-direction, namely, u‘(r, yj+2) and u_(.r, yj_1) to supply function values at Ii,j-—1 and 332‘,j+2- These two computed function values, together with the original function values “id and aid-+1, as well as two interface jump conditions, can determine u‘(a:,:,y). As usual, this determination is sought in coupling with the determination of u'l’ (xi, y), whose solution has no additional problem. Moreover, we can choose either one of the two interface intersection points along the :1: = .12,- mesh line. The two formulations also differ from each other in dealing with large curva- tures. The fictitious domain approach avoids determining a fictitious value along a discretization direction whenever there are insufficient number of grid values along the direction inside an interface. It makes use of grid values in the other direction to determine the fictitious values. However, the interpolation formulation cannot avoid determining polynomials in all directions at an irregular point because it has to calculate both the .r- and i/-derivatives at the irregular point. It therefore makes use of other polynomials to provide function values outside the interface to complete the determination of each required polynomial. These two approaches might involve different sets of grid values in dealing with a given situation. 2.4.3 Comparison of two formulations In this subsection, we discuss the similarities and difference of the fictitious domain formulation and the interpolation formulation mainly based on 1—D cases. In terms of similarities, both formulations share the same set of irregular points for a given interface geometry and given order of the scheme. Both approaches utilize only the standard (high order) central finite difference (.liscretization and the lowest order inter- face jump conditions. The essential equivalence between the two formulations would be obvious if the fictitious value formulation was also casted in polynomial expres- sions. Referring to Figure 1.1, bridged by the fictitious value f,+1, the solution in the left vicinity of the interface is actually approximated by a Lagrange interpolation polynomial 1 u'tr) = Emit-saws.“ (2.78) k=0 to an accuracy of 0(h3) or 3 fix) = ZLku._3+k+L4f.-+1 (2.79) k=0 to an accuracy of 0(h5). These two polynomials, although not explicitly constructed, are actually approximated in seeking the fictitious values. Moreover, a very important common feature is that, in either 1D or higher dimensions, both approaches make use of only I—D polynomials. This treatment significantly simplifies the scheme. Comparing Eq. (2.78) with Eq. (2.53), or Eq. (2.78) with Eq. (2.64) it can be found that with the fictitious domain approach one only needs to determine two parameters, i.e, the fictitious values f, and f5“, to resolve the approximate solutions in both the left and the right vicinities of the interface. With the new formulation, however, one has to solve for all the six polynomial coefficients to determine the approximate solution. Since approximate polynomials are explicitly determined in the new formulation, they are then directly differentiated to provide the approximation to the first and second derivatives in the elliptic equation. In the fictitious domain approach, instead of differentiating approximate polynomials implicitly determined with the first pair of fictitious values f,- and fi+1e a formal central difference scheme involving f, or fi+1 is adopted at each irregular point to approximate the partial derivative. Moreover, to support a fourth or higher order scheme, more fictitious values are needed, and they have to be solved progressively after determining f,- and fi+1- With the new formulation, however, all the coefficients in the high order polynomials are solved simultaneously. 46 2.5 Numerical experiments on irregular interfaces The performance of the MIB scheme for 2-D elliptic problems with irregular interfaces is to be examined in this section with a couple of examples. We will primarily focus on the second-order and fourth-order MIB method in our numerical experiments, while the last case is devoted to the validation of high order MIB schemes. Numerical results are compared to the analytical solutions of the equations, in terms of both numerical accuracy and computational efficiency. The IIM of LeVeque and Li [42] is regenerated for a comparison in some test cases. The performance of our IIM code has been verified with that in the literature [42] and is found to be similar to a later version of the IIM given by Li and Ito [46]. The PBCG solver [65] is adopted to solve the linear system due to its efficiency and the simplicity in implementation. The standard LOO norm error measurement is employed in this section. Example 2.5.1. We consider a 2-D Poisson equation (fiu'I)I + (Wyn, = ftr, 31) (2-80) defined in a square {—1, 1] x {—1, 1] with a circular interface r2 E .172 + 312 = 3 inside. Following [42], the exact solution is designed to be :r:2+y2 r3053 4 711(1— 81?} ’11;)+(I2‘ + T2)/b Otherwise b other-wise 47 Table 2.1: Numerical efficiency test of the 2-D Poisson equation in Example 2.5.1. b = 10, [u] = 0, [Gun] = —0.75. n x n Second-order MIB IIM x 9 L00 Order Loo Order 20 x 20 2.852 (-4) 2.167 (—3) 40 x 40 7.707 (-5) 1.9 5.000 (4) 2.1 80 X 80 2.069 (~5) 1.9 1.131 (-4) 2.1 160 x 160 5.131 (-6) 2.0 2.748 (-5) 2.0 320 x 320 1.257 (6) 2.0 6.781 (-6) 2.0 Such a designated solution forces the discontinuous inhomogeneous term f (3:, y) to be 8.0 1‘ S 0.5 ./(:r, y) = - 8(1):2 + 312) + 4.0 otherwise Let b = 10 such that u(:r,y) is continuous throughout the domain and [Bun] = —0.75 on the interface. The computed result with a 20 x 20 mesh is plotted in Figure 2.11. Table 2.1 lists the computed error of the second-order MIB scheme in a comparison with the results of the IIM. Both methods have very clear second-order accuracy. The MIB delivers a more accurate result than that of the IIM. It is found that the CPU times used for generating the local immersed grids are almost the same for both methods although the underlying algorithms are different. The CPU times used for solving the linear algebraic equation systems are not compared because the IIM method has its own optimal solver [46], while no such solver is available for the MIB yet. It is expected that the IIM method is faster since it involves fewer irregular nodes. Example 2.5.2. We want to use our forth-order MIB method to solve the same Pois— son equation in erample(2.5.1) with a designated finite jump in [u] at the interface: 2:2 + y2 — 1 r g 0.5 um 11) = 4 211(1— glg — %) + (12- + 7‘2)/b otherwise It can be checked that on the interface [u] = 1 and [Han] =2 —0.75 with b = 10. 48 -1 -1 Figure 2.11: Computed solution (upper) and the error (lower) for the 2-D Poisson equation in Example 2.5.1. ,6 = 10. [u] = 0, [Ban] = —0.75. 49 Because the analytical solution is a fourth-order polynomial, it is expected that a fourth-order method shall render a numerical solution of the machine error. The error plot in Figure 2.12 exactly illustrates this prediction, where the maximum error is around 10‘”. 3. O, T‘VYZI‘,‘ ; ‘04“)? ‘ . "K‘VA;"'; 0 Figure 2.12: The computed solution (upper) and the error (lower) for Example 2.5.2. 50 Example 2.5.3. In this example we solve the Laplace equation u“; + uyy = 0, which is defined in the square [—1,1] x [—1,l] and has the following analytical solution: (,2: cos(y) r S 0.5 "(13 y) = 0 otherwise The jumps in u and on along the interface r = % can be evaluated from this solution. Note that here we have unit diffusion coefficient throughout the whole domain. Table 2.2 gives the computational results of the second-order MIB and the IIM for a comparison. Again, as expected, a steady second-order convergence is validated for both methods. Moreover, the MIB scheme is favored over IIM due to the smaller numerical error on all the 5 successively refined meshes. The computed result in Figure 2.13 sharply features designated discontinuity along the circular interface. Table 2.2: Numerical efficiency test of the 2-D Laplace equation in Example 2.5.3. (3 = 1. nr x ny Second—order MIB IIM ‘ Loo Order Loo Order 20 x 20 1.015 (-4) 4.389 (-4) 40 x 40 2.511 (-5) 2.0 1.079 (-4) 2.0 80 x 80 6.369 (-6) 2.0 2.778 (—5) 2.0 160 x 160 1.608 (-6) 2.0 7.500 (-6) 1.9 320 x 320 3.714 (-7) 2.1 1.740 (-6) 2.1 Example 2.5.4. We consider the Poisson equation with a computational domain [—1,1] x {—1.1} and an elliptical interface x 2 y 2 (IS/27) +(10/27) :1 9' ' .0 o , . 0:9 0 .' 0.. 0.. 0d ‘.' ‘0 O... O \I 63" a? ‘0 0.. O Q‘. 9“... . O... Figure 2.13: Computed solution (upper) and the error (lower) for the 2-D Laplace equation in Example 2.5.3. ,8 E 1. The analytical solution and the coefficient ,8 are given as the follows ex cos(y) inside P MI, 11) = 2 5exp (—x2 — 32—) otherwise b inside P fi($,y) == 1 otherwise. Two cases are considered, one with b = 10 and the other with b = 1000. The latter shows a strong discontinuity in the coefficient ,6 and demands more iterations in solving the linear system as it is ill-conditioned due to the large jump in 6. The lower accuracy for the case with b = 1000 can be attributed to the larger jump in the coefficient. Table 2.3: Numerical convergence test and accuracy test for Example 2.5.4. n : b = 10 b = 1000 ”I 4th MIB 2nd MIB 4th MIB 2“ MIB y Loo Order Loo Order Loo Order Loo Order 40 39213-5 5.21E—3 6.19E—3 2.76E—2 80 29213-6 3.75 1.49E—3 1.8 2.65E—4 4.55 7.52E—3 1.9 160 1.70E—7 4.10 3.75E—4 2.0 1.33E—5 4.32 2.17E—3 2.0 320 8.57E—9 4.31 7.80E—5 2.3 6.73E-7 4.30 4.84E—4 2.2 To validate the asymptotic behavior of the approximation error with the increasing of the jump in ,6, a series of 6+ are chosen for numerical tests. The results are collected in Table 2.4. It can be seen that for moderate magnitude the numerical error is increasing while for large jumps the error is almost a constant. Example 2.5.5. To examine the resolution of the proposed method, a Helmholtz-like Table 2.4: Robustness test of high order MIB scheme with Example 2.5.4. 0‘ = 1. e+ 10 20 100 500 103 104 105 108 Loo 3.92E-5 7.10E—5 1.89E—4 3.04E—4 3.28E—4 3.53E—4 3.56E—4 3.56E—4 O O '0 "3:90 0"}! II ’I/ I’ll] I], I’ll/Illllllll I III/ gill/11W?” ////,/,” ,, equation . \h Ir .1 as N 0.015 lllllllll‘llllllllll‘llllllll‘ 1‘ ll 1 11‘ um ‘\l\‘ 0.01 \\\ \\\ t l/////’lllllllllifll‘fllllllllllll °~°°5 I I”/// Illjlblthlll‘l‘fl‘ll \\\\\\\\‘l\“\‘~ 0 l \\\\\\\\\ wk 1 \ \\ \\\\\\\\i\\\\._ I I,” ” "0N“ ““\“ “\‘““ \\\\\\\\\\‘-;;. I III 0,1561%“\\\“\\\\\\\\\\“\\“\“ 1 _0_5 blgg::\\‘ y ‘1 '1 x Figure 2.14: The computed solution (upper) and the error (lower) for Example 2.5.4 with second-order method and n; = 72y: 40. v - (Wutr. y» + W. one. y) = «m. y) is considered, where k(.r, y) = no(x,y) is the dielectric function describing macro- 54 scopically the properties of the medium in which the wave propagates. Both 6 and k are discontinuous at the interface. The analytical solution to the equation is designed to be highly oscillatory x2 + y2 r S 0.5 14$. y) = sin(nx) cos(ny) otherwise, where n can be tuned to produce solutions of desired frequency. The computational domain and the interface are the same as those in Example(2. 5. 2), while 10 r S 0.5 .3 = 1 otherwise, and 1 r g 0.5 0(1. y) = \/1—0 otherwise. It is generally believed that a high order method usually comes with high resolu- tion, which is absolutely needed in solving the Maxwell’s equation or the Helmholtz equation for high frequency wave propagation and scattering. In this example two different frequencies, H. = 2 and K. = 12, are considered as in Table 2.5 where the solution of the present fourth-order scheme is compared with that of the second—order MIB scheme. It is observed that the order of convergence agrees with the theoretical analysis for both schemes. Also observed is that the low frequency solution can be well approximated by both schemes. However, for the high frequency case, considerable differences can be found in the two schemes. In particular, for the low frequency case, a sparse mesh (20 x 20) is sufficient for both schemes to produce a result of moderate 1;}! CJ'! accuracy, i.e., around 10‘3. For the high frequency problem, since the solution ad- mits large gradients and is anisotropic, one has to use a very dense mesh (320 x 320) to resolve the fast variation in the solution if one sticks to a second-order method. It is noted that by using the fourth-order MIB scheme, nevertheless, an 80 x 80 mesh can provide a sufficient resolution. It is therefore anticipated that with the proposed high-order interface method, significant saving on computing time can be achieved for problems that involve both material interface and high frequency oscillations. Example 2.5.6. This example was introduced in Ref. [43], and is adopted here to examine the flexibility of MIB in dealing with complex interfaces. The analytical solution to the equation, the coefficient 5 and a jigsaw puzzle-like interface F are given below (217(y2 + x2 sin(y)) inside P u(x, y) = —(x2 + y2) otherwise, 1 inside P 15(1'. y) = 10 otherwise, 17(0) 2 0.6 cos(6) - 0.3 cos(36) P : y(6) = 1.5 + 0.7 sin(6) — 0.07 sin(36) + 0.2 sin(76). A discretization of the interface is plotted in, Figure 2.9. The computed solution and the error for a 100 x 100 mesh are plotted in Figure 2.16 and the numerical error in terms of the L2 norm are collected in Table 2.6, together with the data of the second-order counterpart. We note that the predicted convergence rate for both methods are verified, whereas the fourth-order method gives a much more accurate result. The maximum error occurs at the irregular points near the interface where the local truncation error is one-order lower than that. at regular points. Example 2.5.7. This is another standard test case for testing numerical methods designed for solving elliptic interface problems. The interface is pararnc—ztrizcd with 56 \ ‘ ¢ ‘ ' 7?. ‘\‘ . . ;::_:::«. , . ; v.‘- “f“ ~i..., .. (l; . 0.03 H" 0.02. it, 0.01~ 13f? _. - Ill/[’6 '00,,” ,4, 7; 3,3;13’00’1’“; “\ Os J;1:..,y;"oI/II$\/“ W‘I‘. \\\ 1f;;.-.,o\\;;,m"o:,’o 4"“ ” fill/2:. %¢‘“““ .: \\\“ ‘63-“ . ., \\\/0‘ " H38“ 1 1233”“. :3 \‘\ u (I; “v": , 39". Law Y '1 '1 x Figure 2.15: The computed solution (upper) and the error (lower) for Example2.5.5. the polar angle 6. as 1 sin(6) — §+—7—-. The exact solution to the problem can be arbitrarily designed. Here we choose exp(x2 + y2) inside P “(13 21) = 0.1(x2 + 312)2 — 0.011n(2 172 + y2) otherwise, 1 inside P 8(x, y) = 10 otherwise This solution also prescribes the Dirichlet boundary condition of the problem. The non-homogeneous term of the Poisson equation can also be derived from the exact solution. Figure 2.16 plots the computed solution with a mesh of 100 x 100. Table 2.7 shows the results of the numerical accuracy tests on three successively refined meshes, in comparison with the second-order MIB method. In the last example, we will scrutinize the convergence of the fourth-order and sixth-order MIB methods. We note that the present method is particularly favorable if the immersed inner boundary is convex. In such cases, there is always an adequate number of grid points to support a given high order approximation of the interface conditions. Example 2.5.8. Here we choose the immersed inner boundary as a circle of radius % with its center located at the origin. The computational domain is {—1, 1] x {—1, 1], excluding the inner circle. The exact solution of the Poisson equation is set to be u(x, y) = 5exp(—x2 — y2/2) on the computational domain. We plot the numerical solution, computed on a 40 x 40 grid, in Figure 2.17 with a negative sign (i.e., —u(x, y)) to illustrate the immersed inner boundary. The accuracy and convergence order of our MIB method are presented in Table 2.8. For high order methods, it takes more time in the generation of the local difference schemes because of the iterative nature of the current algorithm, as given in Table 2.8. However, this time is still negligible compared to the time spent on the solution of the main linear system. It is seen that the result obtained by using the sixth order scheme at the grid of 40 x 40 is about 30 times more accurate than that obtained by using the second-order scheme at the grid of 640 x 640, while consuming only a fraction ( 1 / 473) of CPU time. In order for the second method to reach the same accuracy of the sixth order scheme at the grid of 40 x 40, a grid larger than 2600 x 2600 has to been used, which approximately costs 300000 seconds CPU time. Therefore, the sixth order scheme is about 28000 times more efficient than the second-order scheme in terms of CPU time. It is noted that the CPU time used for the generation of the local finite difference scheme (i.e, the immersed grids) is very small compared to that for solving the linear algebraic equation system. Q-‘DQW "' -‘ flr;‘113.1' as was? as $3 2s seems as saw: can a: 38.2 as «was a? 3:3 as mass 2: was was...” as was? as also: as 32.3 cm as Hag? 2: seems a? mass :3 mags 2. i. as: Haas 33.0. 33.... cm a 85 8a $20 8a $20 8a $20 8a 925% E: a... me: as E: it s: u a: an: an: .m.m.m oEEexm .28 $3 >035an use $8 oozowamZSQ 1.288: Z 6d @738 Figure 2.16: Upper: the computed solution of Example(2.5.6): Lower: computed solution of Example(2.5.7). Both are computed with mesh ul- : ny : 100. 61 Table 2.6: Numerical convergence test for Example 4. n _ n 4th MIB 29d MIB 3 — y L2 Order L2 Order 100 2.10E—8 3.78E—5 200 1.54119 3.77 1.06E—5 1.83 400 5.111111 3.95 2.49E—6 2.09 Table 2.7: Numerical convergence test for Example 5. n _ n 4th MIB 2nd MIB 3 _ y L2 Order L2 Order 100 1.29E—7 ‘ 5.03E-5 200 1.17E-8 3.46 1.35E—5 1.90 400 7.09E—10 4.04 3.41E—6 1.99 -14 2’”: i ' ’I ’I '1 II ’II I ’I ”III, -2~ fifty/”Mgr “If/001W; '1,,,/III/////////;;/,’/” / ' ll/I/I/I/I/Io/ _ ,/ 3‘ ”Ill/Il/Illllfigr III/644’» / , - 4 I ll, ' - a ‘ ‘. .§ III ’, - 5 ) \\ ““.....:::’t , 1 “‘22.:fo'4'. . Figure 2.17: The solution of the high order MIB method for the 2-D Poisson equa- tion in Example 2.5.8. —u(x, y) is plotted. 62 33:5... as 6x a: ea x Se 39:85 38 81 m? 88.83% was ea Q; can x can gasses. Se 31 a: awewsa 2% a; mas c.3888 :3 a: a: 2: x 9: 835.8. a; a: «3 seems? as 9: Be 8392.: as 3.x m3 cm x cm A3332 :3 as 5.34:.” mos a: :3 38.302; as cl as a. x a. l A8332 3x 5% $8,032.15 A3 of cm x om. Duo 820 8a Duo 885 88 :88 SEC 8o s: x 9: m3: uQUuOufififlm mug uQUuOufipuzom fizz nmfiuonficoowm .8828 8:88th $03 We :ofleaozow ofi Sm £0258.» E 8:5 DAD 23 Be. 88.5283 2: E 235:: 25. .8 0.988 :osswavm cemmmom 2: 8“ 88.538 a2 umwaoéufim Use -558 .écoomm 9: mo :OmEeQEoU “wd @358 63 Chapter 3 MIB Method: 3-D Formulation and Numerical Experiments 3.1 Second-Order 3-D MIB Scheme for Irregular Interfaces Consider a 3-D elliptic interface problem V ' (fiVUO‘D - K-(I‘MO‘) = 00‘) (3-1) with interface conditions [11] = u+—-u_, (3.2) [dug] = fl+ug—/3_ug. (3.3) Here the coefficient function 19(1) and the source function p(r) are given. The linear Poisson-Boltzmann equation to be derived in the next chapter is a special case of Eq. (3.1). As the condition (3.3) is defined in the normal direction 75, a local coordinate transformation is needed to relate this normal gradient 115 to the partial derivatives 64 with respect to the three Cartesian coordinates um, uy, uz, (11{,u,7,u<)T = A-(ux,11y,uz)T (3.4) where the superscript T denotes the transpose of the transformation matrix A which is defined as cos it cos 6 cos d1 sin 6 sin w A = — sin 6 cos 6 0 - (3'5) — sin 1/) cos 6 — sin 1/) sin 6 cos 1b The coordinate transformation represented by A is accomplished through two steps. In the first step, the x — y plane is rotated with respect to the z— axis by an angle 6 such that the normal vector would be located on the new x — 2 plane. Denoting this new coordinate by (x', y’, z’), in the second step, the x’ — 2’ plane is rotated with respect to the y’—axis by an angle 11} such that the normal vector is aligned with the x’ — axis after rotation. The coordinate after the second rotation is essentially (it n? ()1 Since cos 111 cos 6 cos d1 sin 6 sin 1b cos 1b 0 sin 111 -— sin 6 cos 6 0 = 0 1 0 — sin w cos 6 — sin w sin 6 cos 1,: — sin 1b 0 cos 1!) c056 sin6 0 X —sin6 0036 0 1 (3-6) 0 0 1 where in the right—hand side the first matrix is the transformation matrix for the second rotation with respect to the y'— axis and the second matrix represents the first rotation with respect to the z—axis. For points located on the contact surface, it is pretty easy to obtain these two angles 6,1!) as in this case 5 is actually the outer normal direction of a sphere with the given center. Details for computing these angles will be presented in the next section. W ith local transformation (3.4), interface condition (3.3) can be recasted as [13115] = ,6+(cos 1b cos6u; +coswsin6uy+ +sinwuj) — ,6“ (cos w cos 6a; + cos w sin 6n; + sin wuz— ). (3.7) Nevertheless, one can not simply use (3.7) and (3.2) as the two conditions to solve for the fictitious values fi,j,ks f,+1, 33k in the x—direction because Eq. (3.7) involves four partial derivatives in the y- and 2- directions which can not be approximated using only the grid points in the x- mesh line. To reduce the number of partial derivatives in the interface conditions, we can differentiate interface condition (3.2) with respect to the tangent and binormal coordinates respectively to provide two additional interface conditions 11:; — u; = (— sin 611},L + cos 611;) — (— sin6u; + cos 6ug), (3.8) 112' — uC— : (— sin 1,?) cos 6n: — sin w sin 611; + cos 9511:) _ (_ sin 1; cos 611; — sin wsin 611,7 + cos 1111;), (3-9) and then cancel any two derivatives which are not in the primary direction from Eqs. (3.7),(3.8) and (3.9). It turns out that there are 12 possible combinations for such cancellation, 4 for x— being the primary direction, 4 for y— and 4 for z—, which are listed as follows 1. Cancel a; and u; from Eqs. (3.7),(3.8) and (3.9) to have [13116] — [3- tanu’1[uC] + (Mun) = C311: + C1711; + (jug) + C211: (3.10) 66 where 2. Cancel u; = 58“ tan6/ cos d1 = 6+ cos w cos6 + [3— tan 1b sinw/ cos6 + 6_ cosw sin6 tan6 = —fi—/(cos1,bcos6) = (8+ —/3_)coswsin6 = (13+ — 71‘) sin ,1 and u; from Eqs. (3.7),(3.8) and (3.9) to have [13nd + C1,[u,,] — 6+ tan[uC] = Cfu: + C; u; + (7;qu + Cz-u; (3.11) where 3. Cancel u; = —,8"‘ cos it" tan 6 — 6+ tan w sin w tan 6 = 13+ cos it cos 6 + 13— cos it" sin 6 tan 6 + 13+ tan w sin 111/ cos 6 2 —(/3+ tanwsinw + 6‘ cos 1/1)/cos6 : (6+ — fi—) coswsin6 = (71+ — [1“)51111/1 and u: from Eqs. (3.7),(3.8) and (3.9) to have [13115] + Cn[u,,] — ,6“ tan Mac] 2 C’jlfuj;r + 0,711,: + C3711; + €211;F (3.12) 67 where C7) = —[3+ cos d1 tan 6 —- 1’3— tan w sin 1p tan 6 C; = 03+ cos w +13— tan If; sin w)/ cos 6 C; = —/3— cos w cos6 — 6+ cos 1.1 sin 6 tan6 — 13— tan 1b sin w/ cos6 Cy— : (6+ — 13-)coswsin6 0: = (1.1+—1r)sinw 4. Cancel u; and u: from Eqs. (3.7),(3.8) and (3.9) to have [181%] + Cn[u,,] — [3+ tan Mac] = C311: + C1711; + Cy—uy— + Cgu; (3.13) where C,, 2 —,(3+ tan 6 / cos 1!) C; = 13+/(cos1t'cos 6) C; = —,/3‘" cos 1!) cos 6 — [3+ cos d sin 6 tan 6 — 13+ tan d1 sin it'/ cos 6 Cy— : (13+ — {3')coswsin6 (I; 2 (13+ — [3“)31111/1 5. Cancel u; and u; from Eqs. (3.7),(3.8) and (3.9) to have [Bug] + C7,[u,,] — {3_ tan w[uC] 2 C31; + Cy— u; + Cgu: + (7:11.: (3.14) 68 where CD = 6‘ cot 6/ cosw C; = 6+ cos '1]; sin6 + 6‘ cosw cos6cot6 + 6‘ tanw sin w/ sin6 Cy- : —6—/(cos¢sin 6) C: = (6+ — 6-)coswcos6 0: = (13+—/i‘)sinw 6. Cancel u; and u: from Eqs. (3.7),(3.8) and (3.9) to have [Bug] + C,,[u,,] — 13+ tan 1,6[ud = Cju; + CJu; + 0:21;; + C;u; (3.15) where C,, = 16 _ cos d! cot 6 + {3+ tan :1) sin in cot 6 C; 2 6+ cos ”69 sin 6 + 13- cos u") cos 6 cot 6 + {3+ tan 112 sin w/ sin 6 Cy— : —(J6_ cos w + 6+ tan 1,0 sin d1)/ sin6 C; 2 (6+ — /3_)cosu’icos6 (I; = W — n-i siin 7. Cancel u: and u; from Eqs. (3.7),(3.8) and (3.9) to have [Hag] + 0,,[u,,] — 5' tan u’i['uC] = 0511'; + Cy—u; + cgn; + Cju: (3.16) 69 where C? 8. Cancel a? = [3+ cos w cot 6 + [3— tan 9; sin 1,0 cot 6 = (16+ c051!) + 6_ tan I/J sin w)/ sin6 = :13— cos iii sin6 — 6+ c0511) cos6 cot 6 — 6’ tan d} sin w/ sin6 = (fl+ —B-)cos¢'cos6 = (6+ — 13-) sis-up and u: from Eqs. (3.7),(3.8) and (3.9) to have [Bud + C7,[u,,] — 6+ tanfluC] 2 C511; + Cgu; + Cat—u; + Cz—uz- (3.17) where 9. Cancel u; 2 {3+ cot 6 / cos :1; 2 6+ / (cos 1,»; sin 6) = 56— cos iii sin6 — [3+ cos ([2 cos 6 cot 6 —- 6+ tan 1,6 sin w/ sin 6 2 (6+ —fi—)cos¢'cos6 2 (6+ — (3‘) sin 1,6 and u; from Eqs. (3.7),(3.8) and (3.9) to have [and + (7,,[u.,,] + Cduc] = Cjuj + Cgug + Cguj + C; a; (3.18) where C,, = 0 CC = 6‘ cotz/J C: 2 [3+ Sim/2 + ,6‘ cosw cot w Cz‘ = —,6‘/sinw C; = (/3+ — fi‘)cos¢cos6 C; = (x3+—/’3‘)cosi/;sin6 10. Cancel u; and a; from Eqs. (3.7),(3.8) and (3.9) to have [Bug] + C,,[u,,] + Cduc] = Cju: + Cz‘u; + C321: + Cgug (3.19) where C,, : —(/3+ — 113‘) cos 1;”) sin6cos6 CC 2 6+ cot 25' sin2 6 + [3‘ cot 1,0 cos2 6 C: 2 13+ Sim/2 +/3+ coswcotwsin26 +fi‘ c0811: COtt/JCOS26 C.‘ = — [1‘ sin ‘(/) — [3+ cos d) cot 7,6 sin2 6 -— ,(3‘ cos w cot 11) cos2 6 C: 2 (6+ — ,6‘)cosz/icos6 C‘ = (,8+—fi‘)coswsin6 11. Cancel a} and u; from Eqs. (3.7),(3.8) and (3.9) to have [Hug] + C77[u.,7] + CduC] = Cj'u;Jr + Cz‘u; + C;u; + Cju; (3.20) 71 where 0,, 2 (6+ — 6‘) cosy’isin6cos6 CC = 6 + cot 11' cos2 6 + 6‘ cot 2p sin2 6 C: = 6+ sin 1;, + 6+ cos 111 cotz/Jcosz6 + 6‘ cos 1/2 cotwsin26 Cz‘ = —6‘ sinzg’i — 6+ cos 1,!) cot 1,11 cos26 — 6‘ coswcotzpsin26 C; = (6+ —— 6‘) cos w cos6 C; = (/3+—[3—)COS't/)Sin0 12. Cancel a} and 11.7; from Eqs. (3.7),(3.8) and (3.9) to have [6116] + C,,[u,,] + Cduc] = Cjuj + Cgu; + Cgu; + Cy‘u; (3.21) where C,, = 0 CC = 6+ cot 1;! C 2‘ = 6 + / sin 11’) Cz‘ = —6‘ sin 6: — 6+ cos 2/2 cot V) C; = (6+—6‘)cosz/)cos6 C‘ = (6+—6‘)coswsin6 One of these twelve conditions will be chosen in addition to Eq.(3.2) for the solution of the a pair of fictitious values. The details of this solution will be illustrated in the 2:— direction with interface conditions (3.2) and (3.10). The solution with the other conditions can be accomplished similarly. Consider a situation as shown in Figure 3.1, where the interface intersects a .1:— mesh line at .s z (.130, ya, 20) between (i, j, k) and (1+ 1.j,ls:), which are then classified as irregular points. The solution of interior domain will be extended to the point 72 n2 1 k+l k 1” I / (OJ—2,10 Figure 3.1: Local topology around irregular point (2', j, k). The interface crosses the sis-mesh line at the point 5 = (1:0, yo. 20) between (1,)", k) and (z' + 1, j, k), which are two irregular points, on which two fictitious values, fi,j,k and fi+1,j,k are defined (marked with green dots). Two auxiliary lines (dashed line) are sketched passing through S, one on the :1: — :1; plane and the other on the TL‘ - 2 plane. Two auxiliary points (in empty circle) (0, j, k + 1) and (0, j, k + 2) are placed on the auxiliary line on the x — 2 plane to facilitate the'discretization of uz‘; Also, two auxiliary points (0, j, k + 1) and (o, j,k + 2) are placed on the auxiliary line on the :1: — y plane to facilitate the discretization of 119‘. (i + 1, j, k) as fi+1‘j.k whereas the exterior domain will be continuated to the point (i,j, k) as fi,j.k- With these two fictitious values the interface condition (3.2) can be approximated to be _ + . . + . . + . . (wonfidik + w0.i+1“2+1,J.k + “’0.i+2“z+2,J-k) (116-1-11% ”1‘3“. + “flint“ + l"’(-):i+1fi+l.j.k) = [u]. (3.22) In a similar way we can also get the approximations to the partial derivatives 11; and u: in Eq. (3.10). The approximation of the other two derivatives 11,7 and u; at the point (10.310. :0), however, remains a major problem because (1'0, yo, :0) is not 73 located on any y— mesh line or z— mesh line. To overcome the difficulty due to the lack of real grid points for finite difference, we add two auxiliary lines passing (1:0, ya, 20), one on the .1: — y plane and the other on the .1: — 2 plane, and set two auxiliary points on each auxiliary line to support a one-sided finite difference for uy and uz‘. This makes it possible to approximate Eq. (3.10) as - i, __ + + + + [(3115] — (3 WWW] + Cnlunl - Cm (u'1,,-fi,j,k + wl,i+1ui+l,j.k + w1,i+2ui+2,j,k) + Cs‘ (wii—iui—mn + Wham + wii+1fi+1.j.k) + ()3!— (povjluosjvk + p0,j—1“'0,j-l,k + p0,j-2“'0,j-‘2,k) + Cz_ (qo,kuo,j,k + qo,k+luo,j,k+l + 90,k+2uo,j,k+2)i (3.23) where p0,; are the finite difference weights for 11.3; at auxiliary nodes (0, j,k), (0, j — 1, k), (0, j —- 2, k) and go; are the finite difference weights for u; at auxiliary nodes (o,j,k), (o,j,k + 1) and (o,j,k + 2). Here 11+(0,j, k) is the limiting value of the solution u(r. y, 2) at the point S from the + side of the interface. Values fz-J-JC and fi+1,j,k can be solved from Eqs. (3.22) and (3.23) as the linear combinations of the solution values at the real grid points u,_11j‘k,11,,j,k,u,-’j,k,u,-,j,k and the auxiliary points uo.j,ki uo,j—l,ki Haj—2J0 uo,j,kv uo,j,k+l: uo,j,k+2 as well as the four given jumps [11.],[6115],['un],['u,<]. These values are given in the following general form fz',j,k = (l”'i—l.j,k + (?'§"li,j,k + (73"22+1,j,k + (-"2i'"i+2,j,k + (75%).}.- + Cfiuo.j+l,k + Cillo.j+2.k + 0511031 + Csiuo,j,k—1 + (/'iOILO,J-.k__2(x'{l[u] + C7{2[6u{] + Ci3[u,,] + Chluc], (3.24) fi+1,j,k = CiHUi—LM + C§+lui,j.k + C§+lui+1,j,k + 6?thsz + Cg+luo.j.k + Cé+17‘().j+l.k + (7;+17‘0.j+2.k + Céi-lluonc + Cé+1”0.j.k—l + (13111,, “.42 + off] [11.] + ('{glt-iug + (i'{§1[u,,] + (731 [11C]. (3.25) 74 In these two linear representations of the fictitious values, the solution values at auxil- iary points have not yet been defined. In MIB, we obtain these values by interpolating the available normal grid points on the same subdomain. For example, 11(0, j — 1, k) can be interpolated on u(i — 1,j — 1,k),u(i,j -— List) and u(2' + 1,j —— 1,k), and u(o,j — 2, k) can be interpolated 011 u(i-1,j— 2, k), u(z',j — 2, k) and 11(2' + 1,j- 2, k). Similarly, auxiliary values 12(0. 3', k+1) and u(o, j, k+2) are interpolated at u(i—2, j, 15+ 1), u(2' — 1,j, k +1),u(i,j,k +1) and u(z' — 1,j, k + 2),u(z',j, k + 2), u(i + l,j, k + 2), respectively. It is also noticed that the auxiliary value u+ (0, j, k) can be interpo- lated as (“E-i-[isjsk + u’fiti-Flu’H-lsjik + “’oti+2”i+2,j,k) or (ul(ii_l'll.i_1,j,k + main“); + 1116: 2' +1 f)“, j,k) + [u], see Eq.(3.22). The expansion coefficients of the fictitious values on these auxiliary points can therefore be distributed to their respective interpolation nodes. Let U(0=J} k) = (ufii_1Ui—1,j,k + ”’01“ij + w()_,i+1fzi+l,j,k) + [u], u(o,j — 1,1.) = I,-_1 . (11(2' — 1,) — 1,k),u(2l,j — 1,k),u(i+1,j — 1,1:))T, 11(2),) — 2,1.) = 1,-_2 . ("(2 — 1,) — 2,1,), ”(m — 2, 1.), 11(2' + 1,) — 2,k))T, "(viii/n+1) = It.”-(u(-i—2,.i'.k+1),u(ii-lijik+1),'u(vi,j.k+1))T. u(0,j, k + 2) = Ik+2 - (21(1' -— 1,j, k + 2), u(i,j, k + 2),u(i + Lj, k + 2))T, where I is the vector consisting of corresponding interpolation weights. Then fi,j,k = Cilia-14¢ + C2”i,j,k + C3”i+1,j,k + Ccl"i+2,j,k + ( '5’ + Ci‘§)(1Li()‘yi_lzr,_l,J-7k + wit-um), + w()—,z'+1fi+1sjsk + [12]) + egg--1 - (11(2— 1,j— 1,k),u(z',j — 1,k),u(2'+1,j— 1,k))T + C; . Ij_2 - (21(2' — 1,j - 2, k),u(z',j — 2, k), 11(2' + 1,j — 21))T + (7611.“ -(~(i? — 2,11 k + 1),“.(1 — 1,.2', k + 1)» "(17,19 k + 1))T + (11,, ~ 1H2 . (11(1 — 1.311: + 2), 11(i.j, k + 2),u(i + 1.311»: + 2))T + Call] + Ci0[6u{] + Cflll’nl + C‘ightc], (3.26) ‘1 C51 _ 2+1 n+1 n+1 n+1 f2+l,j,k — C1 “i—l,j.k+C2 “2',j,k+C'3 “2+1,j,k+('4 ui+2,j,k + (Cg+1 + (}§+1)(U’(i)_1ui—1.j,k + toil-um”); + u’(;,+1fi+l,j,k + lul) + Cg“ .1,._1 . (“(1 — 1,) — 1,k),u(2',j—1,k),1t(i+1,j— lik))T + C;+1.1j_2.(,,,(.)_1,)— 2, k),u(i,j — 2,k),u(i + 1,2 — 2sk))T + 0;,“ - 1H1 - (22(2' — 2,), k +1),u(i-1,]}k + 1), 11(2)]; k +1))T + (:{3‘1 .1“), - (22(2' — 1,), k + 2), “(77,)“, k + 2), 22(2 + 1,), k + 2))T + Cfi” [u] + of,” [ring] + 631 [11),] + 01:1 [224]. (3.27) This finished the solution of fictitious values f2”,j,k and f2+1,j,k in the :c—direction. The standard finite difference scheme for an- at the points (21, j, k) and (2' + l, j, k) will then be modified accordingly by incorporating these fictitious values (i— fi+1,j,k — 2112),): + Iii—1),}; 21.” = 2 at (i,j,k), (3.28) h 12' '. ~2uv - + -_ .- . a” = 113+ ”2’1"" ”22"” f2 1"“ at (n+1, j,k), (3.29) and replacing f2,j,k and fi+1,j,k with their respective representations in Eq.(3.26) or Eq.(3.27). Although solved in the .zz—(lirection, the fictitious value (2.3)}: or ./'21+1,j, k represents the smooth continuation of the solution on the respective subdomain to the current irregular point. Therefore it can be used for the formulation of difference schemes for any partial derivatives rather than a”. In particular, if the central difference scheme of am, at some irregular point in the subdomain 0+ involves the grid point (2', j, k), the fictitious value f,, j, k can be directly supplied. This makes it possible to handle an interface with large Gaussian curvature such that a fictitious value can only be solved in one particular direction while it will be used for the central difference scheme in other direction(s). Comments 3.1.1. It is noted that tan6, cot 6. tand) 01' cot (f) is nndrfincdfor partic- ular values of 6 or (I). This happens when the normal vector at the intevface is on a coordinate plane, i.e., on the :1? — y, :r — z 07' y — 2' plane. One can then apply the MIB 76 method formulated for 2-D problems in Chapter 2 to solve for the fictitious values in these cases. Further reduction to the 1-D problem might happen if the normal vector is aligned with a coordinate direction. The original two interface condition instead of reorganized ones are to be used to solve for the fictitious values in this case, as described in Section 1.2. w / + k+2 / Q i i r 3 k+1$ if 4 (iJ,k) i+l i+2 + x Figure 3.2: Local topology with interface intersecting at a grid point (i, j, k). Two fictitious values, fi,j,k and fi+1,j,k are defined (marked with green dots). Comments 3.1.2. In case that the interface crosses the mesh line at a grid point, the related fictitious values can also be solved in a similar manner as described in Comm- ment 2.1.2. In Figure 3.2, for example, one fictitious value fi,j,k = u(i, j, k) + [u] and the other fictitious value 6+1”; 1: can be solved with interface condition 3.13, in which 22,"; is to be directly approximated using Uiaj‘k,uni’j_1,k,u.i’j_2,k, uz‘ is to be approri- mated using Ui,j1k,Ui’j7k+l, u,,j,k+2, u; is to be approximated using 22,41,331? ui,j,ki £41,3-}; and u: is to be approximated using fi.2'.k~ u,;+1‘j_k, u,+2,j?k. Auxiliary lines and points are therefore no longer needed. 77 3.2 Numerical Experiments Example 3.2.1. We consider a 3-D Poisson equation (SHIELD + (flux/)3; + (5112)}: : f(xa ya Z) (330) defined in a cube [—1,1] x [—1, 1] x [—1,1] with a spherical interface 1 r2Ez2+y2+z2=Z The exact solution is designed to be (12 + y2 + z2)2 r g 0.5 u(:vi y, z) = (3.31) 2(132 + 312 + ZZ)2 otherwise with the diffusion coefficient 1 r < 0.5 6(23, 2)) = . (3.32) 80 otherwise Such a designated solution forces the discontinuous inhomogeneous term f (:r, y, z) to be 20(2:2 + y? + 22) r g 0.5 MI. y) = . (3.33) 3200(1:2 + y2 + 2.2) otherwise Table(3.l) lists the global errors in the infinity norm. Table 3.1: Numerical convergence test for Example(3.2.1) n), :2 ny Loo Order 20 1.56E-2 40 4.16E—3 1.91 80 8.93E-4 2.22 78 Example 3.2.2. We consider a 3-D Poisson equation (flu-TLC + (fiuyly + (fluids? = f(x7 y: Z) (334) defined in a cube [~1,1] x [—1, 1] x [—1,1] with an ellipsoidal interface :22 y2 22 0—49 + 0.16 + 0.09 The exact solution is designed to be u(x,y, z) : sin(x) cos(y) sin(z) r S 0.5 (3.35) exp(:1:2 + y2 + 22) otherwise with the difiusion coefficient 1 r S 0.5 3(1. 32) = a (3-35) 10 otherwise from which one can derive the inhomogeneous term f (:r) and the jumps [u] and [Bun] required for the solution of the Poisson equation. Table(3.2) lists the global errors in the infinity norm. Table 3.2: Numerical convergence test for example(3.2.2) n) = ny LOO Order 20 9.33E—3 40 2.42E-3 1.95 80 5.98E—4 2.02 3.3 Comparison of MIB and HM It is interesting to examine the fundamental difference (or similarities if any) between MIB and HM since both of them manage to impose continuity conditions across the 79 interface by locally modifying the finite difference schemes according to the interface conditions. MIB uses only the zeroth- and first-order interface conditions, while IIM, by contrast, derives a number of additional second-order interface conditions. In IIM, as the solution on either side of the interface is approximated by a multi-dimensional Taylor expansion with respect to the solution at the interface, it allows(actually re- quires) one to choose in the whole domain the grid points near the host irregular point. For 2-D problem, one needs 10 conditions in HM to fully determine the 10 co- efficients in two third-order Taylor polynomials of two variables. These 10 conditions are chosen to be 4 interface conditions (2 original ones plus 2 derived second-order ones) and 6 independent expansions of Taylor polynomials at 6 grid points. More grid points can be chosen to provide an under-determined linear system which essentially provides necessary space for the optimization of these polynomial coefficients. Such flexibility, however, does not exist for the MIB method as here one always chooses a primary direction to solve for the fictitious values and this is equivalent to the deter- mination of 1-D polynomials, as illustrated in Section 2.4 of interpolation formulation. As a result, one has to choose different sets of grid points for the approximation of partial derivatives with respect to x and those with respect to y, therefore MIB has more grid points involved in the finite difference schemes at irregular points than HM and the resultant linear system might be more ill-conditioned. On the other hand, as the fictitious values are always solved in one direction, the generalization of MIB to higher order is straightforward and the small linear system for fictitious values is always invertible. For the HM method, however, the generalization to higher order might be hindered by the lack of a sufficient number of interface conditions for the determination of more coefficients in higher order multi-dimensional polynomials. For example, to determine two fifth-order 2-D Taylor polynomials, a total of 28 coefficients are to be determined. It remains a problem to derive a sufficient. number of high order interface conditions, and the linear system for the coefficients of these high-order multi-dimel1sional polynomials may be not invertible. 80 Chapter 4 Poisson-Boltzmann Equation and Interface Methods In this chapter, we will discuss the Poisson-Boltzmann equation from the perspective of implicit solvent modeling, of which the fundamental theoretical underpinning is the Debye—Hiickel theory describing the distribution of the electrostatic potential ¢(x) in 3—D space. The importance of the electrostatic potential motivated many interests of scientists over the years and reviews of these efforts appear on a frequent basis. In the first section, we will discuss the major considerations and critical assumptions in the derivation of the Poisson-Boltzmann equation. Three most important mathematical issues will be highlighted in the next section. One of them, the discontinuity of the dielectric constant, will be further addressed, but this section is mainly devoted to the brief review of the analytical solution and the numerical methods for the Poisson- Boltzmann equation. A short discussion on the dielectric interface will be given in the last section as this interface is the position where discontinuity occurs and therefore is critical to the formulation of the MIB method and to the solution of the Poisson- Boltzmann equation. 81 4.1 Implicit Solvent Modeling and the Poisson- Boltzmann Equation The importance of electrostatic interactions to chemistry, physics, material science and biology has been well established. Because of the crucial roles of electrostatic properties in investigating electrostatic binding and solvation energies [23], electro— static steering of ligands toward proteins [67], protein conformational change [73, 85], folding stability [73, 95] and many other important subjects, the accurate evalua- tion of electrostatic properties has therefore always been a major concern in molec- ular/ structural biology. This task is further complicated by the consideration of the solvent environment as most biochemical reactions of biomolecules occur in salty solution or at interfaces between lipid bilayers and aqueous phases. A complete ener- getic description of biomolecular processes therefore has to include the surrounding aqueous solvent environment, either explicitly or implicitly [56]. With an explicit model the biomolecule(s) to be studied are embedded in a large number of solvent molecules. The interaction between two charges is described by Coulomb’s law and the total electrostatic energy of the system is derived as the sum of these mutual Coulomb interactions. The explicit model might provide the most accurate and de- tailed description of the solvent but it also significantly increases the computational cost as the number of solvent particles can be more than 90% for a real simulation [4]. Some quantities, such as solvation energy, converge slowly as one has to average over all degrees of freedom of all solvent particles. These critical limitations reduce the time scale, the size of the molecule and the amount of sampling that can be achieved in practical simulations. An implicit solvent model, however, describes the solvent as a structureless dielectric continuum and replaces the solvent interactions with an equivalent energetic term based on the mean field behavior. Similarly, the biomolecule itself can also be modeled as dielectric media with partial charges at atomic positions. The electrostatic potential induced by the partial charges inside the molecule and the free ions inside the solvent is therefore governed by a Poisson 82 equation: v - (sinivnmi = —— + 0(a)) (4.1) where e(x) is the spatial dependent dielectric constant and 50 is the dielectric per- mittivity of the vacuum. Here pf (x) is the fixed charge density inside the molecule and 0(x) is the charge density of the free ions among the solvent. The distribution of free ions depends on the electrostatic field. Therefore the charge density inside the solvent is a function of ¢(x). In particular, one can relate the density of the ionic particles subjected to the external field with 00(1), the ionic density in the absence of external electrostatic potential field, via 525(1) /kT 0(1) = 00(11), (4-2) Il/kT is used to describe the redistribution of where a Boltzmann distribution e‘E( the charge density in the background of potential of mean force (PMF) E (x), which represents the average effect of the entire system on a single particle. This PMF can be further related with the average electrostatic potential if one assumes that the tendency of free ions toward regions of low electrostatic potential energy is dependent of the magnitude and sign of its charge. With this crucial assumption, usually referred to as Gouy-Chapman or Debye-Hiickel theories of ion distributions, the ionic potential of the mean force is replaced by the product of average electrostatic potential and the charge of the ion. Then the Poisson equation becomes the Poisson-Boltzmann equation: V-(c(x)Vq)(x))=—4—7r pf(.r)+Zciziqe‘zqu(I)/kT/\(x) , (4.3) 60 i where 2',- and z,- are the bulk concentration and valence of ion 2, respectively. q is the unit( proton) charge, k is the Boltzmann constant, T is the absolute temperature and /\(.r) is the accessibility of the ion at position .17. There are a number of variations of 83 this primitive Poisson—Boltzmann equation under different assumptions. For instance, if the potential is weak, equation (4.3) can be linearized to be 7r -c,-z.2 2' x V - (€(I)V¢(I)) = -::—0- [291(2) - 2‘ k1; (M AW] ; (44) while for symmetric ionic solution with bulk concentration c and valency 2, equation (4.3) can be reduced to be v . (e(x)V¢(x)) = —4—" [pf (x) — 2czqsinh (“1““) Am] . (4.5) 60 kT l __ eOkT D — 8'2rczzq2 can be derived from the linear Poisson-Boltzmann equation for symmetric ionic solu- The so—called Debye length tion _ 81rcz2 2 4n V - (2 (.I:)VQ(.2:)) + -—-"—q—-¢(.'r))\(.r) = ——pf(.r)/\(x), (4.6) eOhT 60 8rrc22q2 . . 2 . where T has a dunenslon of 1 / (length) . Moreover, it can be found that 150 7 8nc32 2 v - (we) + ——_"—n = o. (4.7) CohT admits a general solution d)(:r) = Ae‘I/ZD (4.8) which describes the exponential decay of the electrostatic potential in ionic solution. Here k D = 1/lD is the Debye screening constant. There are a number of significant assumptions in the above derivation of the Poisson-Boltzmann equation. The critical examination on these assumptions can provide not only the deep understanding of the BBB but also motivations to the mod- 84 ification of Poisson-Boltzmann theory for other important or specific circumstances. Below we will summarize these assumptions and the modifications or extensions con- cerning these assumptions: 1. Fixed charge density in biomolecules. This assumption does not consider the re—arrangement of the polar and charged groups in the external electrostatic potential field. The average polarizability of the protein is measured by its dielectric constant which also macroscopically characterizes the relaxation of the protein structure in response to the charge perturbation. Therefore the neglection of polarizability of the charged groups in an implicit model would lead to considerable deviation of the dielectric constant inside the biomolecules. Shark [72] and Simonson [76] have demonstrated that a larger dielectric constant should be used in solving the Poisson-Boltzmann equation if polarizability is neglected. 2. Uniform dielectric constant in biomolecules. The issue about the dielectric constant in biomolecules is closely related to the model of fixed charge density. As mentioned above, the dielectric constant of a protein is a macroscopic quantity to characterize its relaxation in response to charged perturbations. Therefore the hydrophobic core and the charged protein side chains at the protein/solvent interface might have different response and different dielectric constant, as the latter are usually undergoing large, mutually uncorrelated and anisotropic motions. In this regard, molecular dynamics (MD) simulations are usually used for the detailed prediction of the protein dielectric properties in solution. Simonson illustrated that the Frohlich-Kirkwood theory could provide a consistent prediction with the MD simulation if the charge por- tion of the side chains are treated as part of the solvent and assign the remaining part of the protein a uniform dielectric constant of 3—4 [76]. A more complicated inhomogeneous dielectric map can be generated by assigning every amino acid a different dielectric constant based on its unique intrinsic polarizability [75, 77]. There are more physical considerations of the effects of polarizability on the di- 85 electric function besides the layered dielectric model for proteins. For example, a local dielectric function (.(r, 2") (depending on the response point r and the perturbation point r’) or 602' — r’ I) (depending on their separation only) rather than a single dielectric constant can be supplied to the Poisson-Boltzmann equa- tion to introduce the polarization induced by the external field, the fluctuations of electronic polarization within an individual molecule, and the local field in a cluster of non-overlapping molecules [34]. Non-local dielectric functions have been used to account for the effects of solvent shells surrounding charged species in work on electron transfer, the double layer at electrode-solvent interface and intermolecular interactions near interfaces. . Singular charge density in biomolecules. The physical modeling of charge distribution inside biomolecules remains an open question. It appears that there is a tendency to use more physical finite- size models such as homogeneously charged sphere, the Fermi two-parameter charge distribution or the Guassian charge distribution instead of the usual point charge model [80]. The advantages of these finite-size models to the solution of the Poisson-Boltzmann equation, however, have not been explored in depth. Particularly, because the characteristic size in these finite-size models is of the same scale as the root-mean-square (RMS) radius of the nucleus ma?) = (0.836A1/3 + 0.570) x 10‘15 meters, it is very difficult to provide a discrete description of these finite-size models on a mesh of size ranging 0.1A to 0.5A. It might be possible to investigate the consequence of these finite-size models in the framework of the boundary element method, whereby the Green function for the point charge can be simply replaced by the smooth potential corresponding to a finite-size model. . Boltzmann distribution of the ions in the solvent. When a biomolecule with charged surface is immersed in an electrolyte solution 86 (ionic solution, for example), ions of opposite charges to the surface are electro— statically attracted to the surface while ions of like charges are repelled. The counter ions are not bounded to the surface but remain dispersed and mobile in the vicinity of the surface due to their thermal motions. According to the Gouy-Chapman—Stern model, the region next to the surface is called the Stern layer where the distribution of ions is determined by the Coulomb interactions and special absorption. The region next to the stern layer is called the diffuse layer and ions in this layer can move freely (diffuse) in any direction. This spatial separation of ions is denoted as the electric double layer [37]. In addition to the local density distribution represented by the electric double layer, the packing effects caused by the finite size of the ion and solvent molecules can also be considerable and the resultant spatial variations of ionic density and strong correlation between ion-ion, ion-solvent or solvent-solvent molecules can not be modeled by a Boltzmann distribution. A number of modifications have been developed to address the deviation of real ionic distribution and the Boltzmann model. For instance, Eigen and Wicke [16] modified the Boltzmann distribution by introducing steric and confinement effects of ions, while Luo et. al. [54] conducted molecular dynamics simulations including the solvent and ionic structures to obtain the potential of mean force (PMF) of a single ion. This PMF is then supplied to the generalized Poisson-Boltzmann equation (combining Eq.(4.2) and Eq.(4.1)) 47r v . vn(s)> = 10% (s) + fan/”noon. (4.9) The solved real ionic distributions were found to agree with the x-ray reflectivity measurements. 87 4.2 Analytical Solutions and Numerical Techniques for the Poisson-Boltzmann Equation: A Brief Review A close look at the Poisson-Boltzmann equation will reveal four major features, namely, the irregular interface between the solute and the solvent, the discontinuous dielectric function on this dielectric interface, the singular partial charge distribution and the nonlinearity in the case the strong potential. We are particularly interested in the handling of the dielectric constant with discontinuity at the irregular molecular surface, not only because dielectric properties of macromolecules are crucial to their stability and activity but also because this discontinuity poses one of the most difficult challenges in seeking the numerical solution of the Poisson-Boltzmann equation and all of the current methods for the Poisson-Boltzmann equation are lacking of specific treatments of this discontinuity. In particular, we are interested in the incorporation of these two conditions ¢p = 052v: (pa; = (.212 (:3 1 (4°10) at the dielectric interface into the numerical solution of the Poisson—Boltzmann equa- tion. The first condition describes the continuity of the electrostatic potential across the interface while the second characterizes the balance of the derivative of the electro- static potential field at the interface. The subscript p denotes protein and w denotes solvent. We will refer to the second as the potential flux continuity condition, which is not considered in any finite-difference or finite-element based Poisson-Boltzmann solvers. Because of the complicated dielectric interface and charge distribution, the ana- lytical solution to the Poisson-Boltzmann equation is not available for almost all the pratical systems 1 and a variety of computational techniques and software packages 1There are very few systems having analytical solutions of linear Poisson-Boltzmann equation. To the author's best. knowledge, they are: 88 have been developed for the numerical solution fo the Poisson-Boltzmann equation. These numerical techniques cover all the major discretization methods such as finite difference, finite element, boundary element and Fourier spectral methods. One of the first such Poisson-Boltzmann solvers was developed by Orttung [62] in 1979 based on the finite element method. After that a number of other Poisson-Boltzmann solvers emerged, such as a finite difference method by Warwicker and Watson [86] and a boundary element technique by Zauhar and Morgan [93]. The successful applica- tions of these solvers, espically the method by Warwicker and Watson [86], on the biological macromolecules, motivated the establishment of the numerical technique of higher accuracy and efficiency. A list of major Poisson-Boltzmann solvers with a short description of each solver was given by Baker and is reproduced here. The construction of accurate and efficient numerical algorithms is restricted mainly by the point charges and the discontinuous dielectric constant at the molecular sur- face. Despite the slow advancement of numerical techniques for the singular charge density, a number of techniques were devised to improve the accuracy of the solution near the dielectric interface. For example, APBS manages to discretize the equation on a mesh a priori refined in the vicinity of the molecular surface. Some robust a posteriori error estimation is used to drive the adaptive mesh refinement when multi- level meshes are used in APBS [27]. Similarly, Delphi [24] attmpts to improve the accuracy of the solution, especially in the vicinity of the dielectric interface, by first solving the Poisson—Boltzmann equation in a computational domain much larger than the molecule(s) enclosed. Then, the second round of computation is carried out on the smaller domain with a refined mesh. The boundary condition of this second 1. A single sphere with single or multiple point charges. They are to be discussed in this dissertation. 2. A single sphere with a uniform distributed charge on the surface [28] 3. Multiple point-charged molecule with complete solvent penetration [28]. 4. Two centered charged spheres. It seems that an analytical solution might be obtained using the bispherical coordinate system. Despite a number of research articles [20, 61, 78] in the past decades, there is no simple solution of electrostatic potential in closed form for this case Some solutions in semi-analytical form are series of slow convergence and are too complicated to be implemented. 5. Systems with cylindrical geometries [57. 74]. 89 run comes from the last computation. Such a ‘focusing’ procedure can be hierarchi- cally implemented to acquire a sufficient resolution and accuracy near the molecular surface. All these techniques are based on the local mesh refinement and there is no special treatment of the interface conditions in the fundamental discretization scheme. A noteworthy class of numerical techniques for the Poisson-Boltzmann equation is the boundary element method. In this approach, two boundary integrals, one in the interior of the molecule and the other in solvent, are defined. The electrostatic potential and its normal gradients are determined from the matching of both inter- face conditions on the molecular surface with these two integrals. Once the surface potential is solved, the potential in the entire domain can be analytically computed. Theoretically, boundary element approach is superior to all the other techniques be— cause it allows the discretization of the equation in lower dimension, as both interface conditions are satisfied, because the singularity of the point charges are removed through boundary integrals, and because the far field boundary conditions are en- forced exactly in an integral. The very first boundary element method proposed by Zauhar and Morgan [93] assumes that the ionic strength is zero but this assumption was soon removed in the formulation of Rashin [66]. The first boundary element approach based on Green’s function for the Poisson-Boltzmann equation, however, probably is due to Juffer et. al. [66]. Compared with the wide applications of finite difference or finite element based Poisson-Boltzmann solvers, the applications of boundary element methods are found to be limited. This might be explained by a number of fundamental issues in the boundary element methods for the Poisson-Boltzmann equation. For example, the Green’s function of the Poisson-Boltzmann equation may be difficult to find if the molecules have complex internal structures. Additional volume integrals have to be calculated if the boundary element method is used to solve the nonlinear Poisson- Boltzmann equation [9]. The final linear system might be ill-conditioned if weak formulation of the linear Poisson-Boltzmann equation or the numerical quadrature in the surface integral is not minimized. In particular, the matching of the two interface 90 conditions with the surface integrals gives rise to 4N2 influence coefficients where N is the number of surface elements. The calculation and storage of this large number of quantities severely limits the efficiency of the boundary element method. It has been demonstrated that the application of a fast multipole technique could reduce the computational complexity of these quantities to an order of N log(N) [8]. It is expected that further investigation on the fast multipole or related method such as tree code, and on the fast calculation of the Maxwell stress tensor on the molecular surface [53] could provide necessary thrust to the application of the boundary element method on the numerical solution of the Poisson-Boltzmann equation. 4.3 Dielectric Interface The definition of dielectric interface is one of critical issues in the implicit solvent model. Any biomolecular surface, such as van der Waals surface, solvent accessible surface, solvent-excluded surface [12], Gaussian surface [88], multiresulution molecular surface [87], can be chosen as the dielectric interface in computation but the choice of which biomolecular surface to use as the dielectric interface is still open to debate. Most of these biomolecular surfaces consider only the constant van der Waals radii of the atoms and the radius of the probe which is the spherical model of the solvent molecule. They are defined a priori and are supplied into a particular implicit solvent model, the Poisson-Boltzmann equation for example, to determine the energetic state of the system. In contrast to these surface definitions, the biomolecular surface can be derived from the minimization of an energy functional of the solute cavity shape [15]. In this case, the biomolecular surface is the output of the implicit solvent model rather than its input. The implementation of a well-defined molecular surface in the solution of the Poisson-Boltzmann equation depends on the underlying numerical technique. In finite difference based Poisson-Boltzmann solvers such as DelPhi, MEAD and PBEQ, the molecular surface is used to define the map of the dielectric function only and the position of the molecular surface is not explicitly computed. In PBEQ, for example, 91 a set of probes centered at the van der Waals surface are defined, and any grid point enclosed in this envelope is regarded to be in the solute (otherwise it is in the solvent). In APBS with finite element discretization, nevertheless, the position of the molecular surface is explicitly calculated to provide the exact position for local mesh refinement. None of these Poisson-Boltzmann solvers calculates the normal direction of the molecular surface as they do not consider the discontinuity of the potential gradient at the molecular surface. Few molecular surface generators calculate the normal direction of the molecular surface. We notice that the normal vector and other geometrical quantities can be derived from the spherical harmonic expansion of the molecular surface [58]. The slow convergence of such expansion toward an accurate description of the molecular surface prevents us from implementing the same molecular surface as other Poisson- Boltzmann solvers. The MSMS program due to Michel Sanner et. al. provides an efficient algorithm for the molecular surface, its normal directions, its total area and its triangulation [68]. We will implement the MSMS surface in our solution of the Poisson-Boltzmann equation. However, we do prefer the dielectric interface defined as the iso—surface of some analytical function such as a Gaussian surface or the recently developed multi-resolution model of the molecular surface. The former is the iso- surface of additive approximate electron density centered at each partial charge and the latter is the iso—surface of the solution to a diffusion equation. Because the surface is a particular level set of an analytical function, the solvent accessibility and the surface geometry can be readily derived from the function [70]. The molecular surface derived from the minimization of an energy functional can also be described by an isosurface hence it can be easily implemented into our interface algorithm. 92 wmm< at? $83: on :do use 630m 02 Oh $55ng 23 3:333 coca .XHZD =< avoéigniamidg no 52 :o £32325 £3 23605 $0533-53“: 22m<=0 2038 00 was .MOm .02 mm 33555» 365588 JQZD “no: 3.3.3:?“\mauavenqaoo.Sm:%oh~on.333\\REE auto €43 :o Ema—ESQ APE 58mg... $0353-53—2 Human. $38 mmm mOm On 352258 858 £30 .85 40m ~83.xovobaEFvedowEHanu.EEQEEE 2me ”am no fiflnaaw ~23 :8...»th 80339532 onoBeE $38 mmm mOm Gm Rsnescem 23333 oomeQZD =< Buffing—Ezvo.vwo=doE§ooE\\duE macho 5m no E22350 533 Seaman emoahaaszz QmmD 850m :30 68m .XHZD :.<. £3355) EoEwB\=vo.maqum.333\\”3.2 $38 mOm Om a 55 323:0:me mmm 830m 0354 mi: -amaoa mm 353523 .3 533238 33538 28 38238 coma 40m amauw\=vo.d553oo.u05d3:a5\\”a3: .wo 602998 :o mmmdcafim £23 5.3on 203325;; mm< ED coiniowea omdxodm oudBdom 93 AQEVSEEEAU $30208 .QfiémeEeaooE 318208 .AZGVmoEefioE 83:38 .AQmeB -8588 osggosm .AMOmvcosaxfiou $>o wiwmmoosm .AOOvfiBEEw maewaaoU .AmOVEoEemwman Ammvuaefiflo 3:5 .Emvoesmcomac SEC AUE< I12 3 4 5 6 7 8 91011121314151617 Figure 5.2: 2—D diagram of the accessibility determined by chasing intersec points of mesh lines and molecular surface. Accessible grid points are marked by solid red circles and the inaccessible grid points are marked by empty black circles. For the upper mesh line (y = y4) which has 6 intersection points with the molecular surface(the right most one is counted twice as the mesh line is tangent to the molecular surface at this point), the grid points ( 1,4) through (3,4) are accessible because they are located before the first intersection point; grid points (4,4) through (5,4) are inaccessible because they are located between the first and the second intersection point; grid points (6,4) through (8,4) are accessible again because they are located between the second and the third intersection point; grid points (9,4) and (10,4) are inaccessible again because they are located between the third and the fourth intersection point; grid points (11,4) through (15,4) are accessible because they are located between the fourth and the fifth intersection point. Between the fifth and the sixth intersection points there are no grid points. After the sixth intersection point which is the last one on this mesh line, all the grid points are accessible. importantly, Table 5.1 provides a guide line for the selection of surface triangulation. In particular, in order to get a second-order convergence up to a mesh size of 0.1A, the density should be larger than 50. On the other hand, if one sticks to a mesh size of 0.5A in the calculation of electrostatic potential for the real biomolecules, a triangulation density of 10 is sufficient since for this mesh size there is no essential difference between the solutions at density of 10 and density of 40. Assured by this observation, an MSMS surface with a density of 10 will be used for numerical so- lutions of the Poisson-Boltzmann equation with a mesh size of 0.5A. When a finer mesh such as 0.2A or 0.1A is necessary for the examination of the convergence, we 99 will instead use the MSMS surface with a density of 40. Our implementation of the Table 5.1: Positions and radii of atoms of cyclohexane. 1‘ y z vdW radius -2.0270 0.9540 -0.6510 1.7 -1.6690 0.2340 0.6650 1.7 -0.4530 -0.6870 0.4410 1.7 0.7510 0.1480 -0.0400 1.7 0.3930 0.8680 -1.3560 1.7 -0.8230 1.7880 -1.I320 1.7 -2.2840 0.2080 -1.4180 1.2 -2.8880 1.6170 -0.4830 1.2 -2.5270 -0.3680 0.9970 1.2 -1.4260 0.9800 1.4350 1.2 -0.1960 -1.l890 1.3850 1.2 -0.7010 -1.4410 -0.3200 1.2 1.0070 0.8940 0.7270 1.2 1.6120 -O.5150 —0.2080 1.2 0.1490 0.1210 -2.l260 1.2 , 1.2510 1.4700 -1.6880 1.2 -1.0810 2.2910 -2.0760 1.2 -0.5750 2.5430 -0.3710 1.2 Table 5.2: Numerical convergence test for Example (5.1.1). Density of Mesh Loo triangulation size(A) error Order 1/2 103(0) 5 1/4 9.83(-1) 1/8 116(0) 1/2 1.38(-1) 10 1/4 3.63(-2) 1.93 1/8 1.65(—2) 1.14 1/2 137(1) 20 1/4 3.134(2) 1.95 1/8 9.83(-3) 1.85 1/2 1.38(-1) 40 1/4 3.54(-2) 1.96 1/8 9.07(-3) 1.96 MSMS surface for the MIB method is also applicable for the HM method, i.e., we first find the intersect points of each mesh line with the interface and then determine the accessibility of each grid point. It is noticed, however, the HM scheme at any irregular point only requires the geometry of a nearby interface point, which is not necessary I 00 to be the intersect point of the mesh line and the molecular surface. We choose one of the nearby intersect points as the required interface point since it produces minimum extra calculation. ____________________________ Figure 5.3: The grid points and interface point needed for the formulation of the HM scheme [13] at the irregular point(i, j, It). The IIM method manages to find the Taylor expansions of the solution accurate to 0(h3) at both sides of the interface, for which a total of 16 expansion coefficients are to be determined. Original IIM requires 10 grid points and 6 interface conditions at an interface point(p’ or p) to solve for the 16 polynomial coefficients. These 10 points(red) are chosen to be (i,j. k), other 6 grid points directly connecting to (i, j, k) and 3 more grid points inside the cube, preferably at the other side of the interface. IIM with optimization employs all the 27 grid points (red and green) inside the cube to produce constrained quadratic optimization problems from which 16 polynomial coefficients are solved. See [13] for the detailed formulation of 3—D IIM scheme. It is not necessary to choose the interface point to be the nearest one as 11'. For simplicity we choose p’, the intersect point of the interface with the mesh line, as the interface point for the irregular point (i, j, k). 5.2 Numerical validation of interface method for the Poisson-Boltzmann equation We consider a sphere of radius 2A, centered at the origin and with charges inside whose number and positions are to be tuned. A sphere with single or multiple charges admits an analytical solution of electrostatic potential and has been used for evaluation and 101 accuracy assessment for various Poisson-Boltzmann solvers. For noncentered charge or multiple charges one can seek the analytical solution as the expansion of spherical harmonic functions [39]. Both the exact boundary conditions and the exact values of the solvation energy can be calculated from these analytical solutions. The derivation of analytical solutions for multiple charges is given in Appendix A. The first computation is on the sphere with a single unit charge. Table 5.3 lists the convergence of two interface methods in comparison with the PBEQ and the APBS. In addition to the electrostatic solvation energy, we also calculate the maximum absolute error e1, the maximum percentage error e2 and the average percentage error e3 of the electrostatic potential near the surface of the sphere, as ~ 81 = max|¢(-'v,y,z) -¢>(:v.y,z)l, (5-4) ¢(x,y,Z) Ni" a > a > __ :r,y,z — 1:,y,z . e3 _ 100x 2| (r,y,z) |/N,,.,.. (5.6) Here (E(zr, y, z) and ¢(.r, y, z) are the computed and exact electrostatic potential, re- spectively. The maximum values and summation are taken over all the irregular grid points near the molecular surface, whose number is given by Ni". Given a suffi- ciently small grid spacing, a convergent solvation energy could be solved with all three methods. The necessary resolution required for a given accuracy of solvation energy, however, differs across the methods due. to their different convergence proper- ties. MIB and HM have a convergence rate of 2, which is indicated by about fourfold decreasing in all three errors as the grid spacing is halved. Such a uniform convergence is attributed to the modified difference scheme near the interface. In contrast, PBEQ and APBS have similar convergence rates of around 0.3. As a result, a grid spacing of 0.5A suffices for both MIB and HM to yield an almost exact solvation energy whereas for PBEQ and APBS one has to use a mesh size smaller than 0.05A to get a solvation energy of similar accuracy. The most. significant difference occurs at the accuracy of the surface potential, which can be seen from the maximum absolute error, the 102 maximum relative errors and the average relative errors in Table 5.3. For mesh size of 0.5A, the surface potential of PBEQ and APBS deviate from the exact value by more than 84% at its maximum and by more than 23% on average. The errors in the surface potential calculated with PBEQ or APBS are consistently much larger than those of the interface method despite a slight decrease with the refinement of the mesh. The maximum error is around 50% even for a grid spacing of 0.05A. With either interface method, however, the surface potential is accurate to about 8% at its maximum and only about 2% on average for a grid spacing of 0.5A. If h = 0.2A is chosen the maximum error is only about 2% with corresponding average error below 1%, which is almost unachievable by PBEQ and APBS due to their slow convergence. The second system we studied is a sphere with a noncentered charge located at (a,0,0). As the charge approaches the surface, a stronger interaction between the charge and the surface is anticipated and a larger error may occur if the two interface conditions are not handled correctly. The electrostatic solvation energies for 0 g a S 1.5 were computed and collected in Table 5.4-5.5. In addition to the growth of the solvation energy these two tables also show that the computation errors, both for PBEQ and the two interface methods, increase as the charge approaches the surface. It is interesting to see that interface methods and PBEQ have different responses to the displacement of the charge. As the charge is moving from the center to (0.2, 0,0) the relative errors, both the maximum and the average, increase about twofold for the MIB method whereas they remain in the same high level for PBEQ. A slight increase in the maximum relative error can also be seen for the HM method. A larger error in the surface potential of both interface methods occurs as the displacement of charge increases but the average relative error is well below 5% except for a = 1.5. For PBEQ, although there is a slight decrease in the accuracy of surface potential, the maximum and average relative errors are always larger than 20% for the various displacements of charge. Furthermore, the average relative error in the interface methods can be efficiently reduced to less than 1% for all displacements by refining the mesh to a grid 103 spacing of 0.2A but the large error in the traditional finite difference method PBEQ decreases only slightly from 23% to 12%. Table 5.5 also shows that a very large error could occur on interface methods as the charge is very close to surface, such as a = 10,12 or 1.5. Such a large error results from the mixing of the grid points for singular charge distribution and for the modified difference schemes near the interface. In particular, the solution on the grid points within the trilinear distribution of the singular charge is of the lowest accuracy. The inclusion of these grid points into the difference scheme at the irregular points, as in the case of a > 0.6 and h = 0.5, will degrade the accuracy of the solution near the interface. However, because the van der Walls radii of the partially charged atoms are usually larger than 1.0A such a small distance between the charge and the surface will not occur in pratical biomolecular simulations. A stronger interaction of charge and surface can be found by placing two noncen- tered charges in a single sphere. In the third system for our validation, one charge is placed at (a, 0, 0) and the other at (0, 0, a) where a varies from 0.2A to 1.5A. Except for the increase of the solvation energy, the results for both interface methods and PBEQ display a striking qualitative similarity with the case with a single noncentered charge, as presented in Table 5.6-5.7. The amplification of the charge-surface inter- action results in a larger maximum relative error in the surface potential of interface methods, which could be 34% for a = 0.6A. However, the average relative error does not grow with the increase of the number of charges, and is smaller than 3% for a < 1.0A, similar to the case with single charge. For all three systems, PBEQ provides a more negative (i.e. more favorable) free energy of solvation, and the deviation from the analytical value is around 5% for a mesh size of 0.5A. Such a systematic percentage error amounts to about 3.87kcal/mol for a centered-charged single sphere and will lead to about a deviation of 17kcal/mol for sphere with two noncentered charges (a. = 0.6). These three numerical experiments also demonstrate some noticeable differences between the two interface methods. Compared with the MIB method, the smaller 104 interface discretization stencil of HM indicates that its interface discretization will less likely be influenced by the singular charge distribution. This speculation is validated by the smaller el and eg of HM in Tables 5.3 - 5.7 at a mesh size of 0.5A if the displacement of the charge is smaller then 0.6A. In particular, the maximum absolute error in the surface potential is 2.63kcal / mol / e for MIB and is about 1.98kcal / mol / e for IIM if the displacement of charge is zero. Substantial increases of the maximum absolute error and the maximum relative error are seen in the MIB method. Such increases for IIM are nevertheless much smaller. This relative advantage of smaller discretization stencils in HM will diminish when the charge is close to the interface. If that is the case, the grid points participating in the distribution of singular charges have to be involved in the discretization at irreg- ular points and these two types of errors are always mixed together. A considerable increase in the maximum absolute error is seen in Table 5.3 and Table 5.6 when the displacement is increased from 0.6A to 0.8A or larger. On the other hand, MIB has a more accurate approximation to the interface conditions and the resultant discretization at the irregular points as illustrated in Table 2.1 and Table 2.2. This explains its smaller average relative errors for all the tested cases in the above numerical experiments, despite its larger maximum absolute errors and large relative errors at small charge displacement. Furthermore, the errors in the solvation energies computed with IIM are much larger than those computed with MIB, especially at large charge displacement. For this reason we will only report the results of MIB in all the other numerical experiments and calculations for biomolecules. For general molecules, however, the evaluation of the solutions of the variety of PB solvers remains difficult due to the lack of analytical solutions for reference. The slow convergence or stagnance of many Poisson-Boltzmann solvers, as found in practical simulations, prevents the achievement of a convergent solution through simple mesh refinement for self-checking. In this numerical experiment, we are going to investigate the convergence of PBEQ and MIB on the problems without analytical solutions. In particular, we are interested in whether the numerical results provided by the MIB method can be used as an alternative reference to measure the accuracy of popular PB solvers. Two unit positively charged spheres, of the same radius 2A, are placed on the r—axis symmetrically with respect to the origin. The molecular surface of this model molecule is composed of two spherical caps at two ends and an inward tori surface in the middle, each of which is analytically given. The solvation energy is computed for systems with different sphere-sphere distances, which vary from 2A to 5A. The electrostatic potential computed from MIB at the finest mesh of h=0.05A is saved as reference to measure the accuracy of the surface potential computed with MIB and PBEQ at three coarse meshes of h = 0.5, 0.2, 0.1. The data given in Table 5.8 indicate that the solvation energy is almost convergent for MIB with a grid spacing of 0.5A and the further refinement of the mesh can only provide minor improvement. For PBEQ, a deviation of 4-6kcal/mol in solvation energy at the coarse mesh can be reduced by mesh refinement. Hence the solution of PBEQ essentially converges to that of MIB as far as the solvation energ is concerned. The convergence rate of PBEQ’S solution to MIB’s solution matches well with the rate observed in the last three systems with analytical solutions as reference. For all the charge distance considered, the maximum and average relative errors of the surface are about more than 80% and 20%, respectively for a grid spacing of 0.5A, which are of the same magnitude as observed in the first three test systems, see Figure 5.4. The average absolute error in the surface potential at this resolution is about 3-4kcal/mol whereas the largest error could be 17-22kcal/mol. 5.3 Application of interface method for biomolecules Motivated by the superior accuracy and convergence property of MIB, we then move to the numerical experiments of MIB on biomolecules. We are going to investigate if traditional Poisson-Boltzmann solvers will converge to the solution of MIB at coarse mesh for large macroinolecules and if we can achieve the savings in the computational 106 till l b 1D[I ? l l I» o I} 1 N01 2 3 4 5 3 4 5 Figure 5.4: Maximum absolute error e1 (upper left), Maximum relative error e2( upper right), average absolute errors(lower left) and average relative error e3 (lower right) in the surface potential of two positively charged spheres. The errors of MIB and PBEQ are plotted in black and blue, respectively. The average absolute error is the total of the absolute errors at all irregular points divided by the number of irregular points. The absolute error is in kcal / mol / e and the relative error is the per- centage. The :r-axis is the distance between centers of two spheres. Circle: h = 0.5; Diamond: h = 0.2; Square: Ii = 0.1. time and the memory due to the use of coarse mesh. Twenty four structures, mostly adopted from test set 1 of [19], were used in this study to examine the accuracy and robustness of MIB in handling different protein conformations. The missing hydrogen atoms in the original PDB files were added using the HBUILD function in CHARMM. The all-atom representations were generated and partial charges were assigned, both with the CHARMM22 force field. The coordinates and partial charges of the atoms prescribed as such were also supplied to APBS. It is noticed that hydrogen atoms in some of the amino groups are completely buried inside the nitrogen atom and MSMS can not handle such cases. These hydrogen atoms then have to be removed 107 before being supplied to MSMS as they have no contribution to the molecular surface. The partial charges at these hydrogen atoms, however, have to be retained in the solution of the Poisson-Boltzmann equation as they do contribute to the electrostatic potential. The deviation of electrostatic solvation energies calculated by PBEQ with different mesh size from those calculated by MIB with a fixed mesh size of 0.5A are plotted in Figure 5.5. For all 24 proteins, PBEQ delivers a more negative solvation energy at a mesh size of 0.5A or 0.25A. At a mesh size of 0.15A, nevertheless, there are two proteins for which the PBEQ gives a slightly larger solvation energy. For protein 1ajj, the solvation energy solved from PBEQ is larger than that from MIB by 0.7kcal/mol or 0.06%. For protein luxc, this difference amounts to 1.5kcal/mol or 0.13%. The overall monotonous decreasing of the deviation from the MIB with the improvement of resolution proves the reliability of MIB at low resolution as the reference for traditional Poisson-Boltzmann solvers. The saving of computational time and memory usage are critical for simulations of large molecules, for which it is usually impossible to use a resolution of 0.2A or smaller to obtain a convergent solution. In particular, there is a strong correlation between this deviation and the total partial charge of the biological molecules. For protein 1a7m with 180 residues, for example, the total charge is about 677e and the deviation is about -103.5, —36.5 and —12.5kcal/mol/e for h = 0.5, 0.25 and 0.15, respectively. A much larger error of solvation energy can be projected for molecules of large size or highly charged, and might contribute to the over stabilization of slat-bridges as found in many simulations [55, 95]. This correlation also explains the aforementioned overshoot for the small proteins 1ajj and luxc. For larger proteins with large total absolute charges, the solvation energies computed from PBEQ at a mesh size of 0.15A is always smaller than those from MIB by more than 10kcal/mol. The surface electrostatic potential of 1a63 generated with PBEQ is shown in F ig- ure 5.6, along with the potential difference between the PBEQ solution and the MIB solution. There are scattered regions of negative potential but most of the potential surface is largely positive where the numerical error is also large. as seen in the right 108 O455755;::aa".:,;,:::r;;;;, 100 -20 5 : '1??? itiiitiiiiii 5 6 "'::E’::$tl'200 “4°” =2 515i? 300 if??? ‘80- 0: ii‘4°° —1oo~ é 'élsoo -120- -600 -140- 700 Figure 5.5: Left y-axis: the deviation of electrostatic solvation energy calculated with PBEQ at h = 0.5 (circle), h = 0.25(diamond) and h = 0.2(square) from that calculated with MIB at h = 0.5. Right y-axis: the solid line plots the total absolute charge of the proteins. x-axis from left to right: 1ajj, 1vii, 2erl, lbbl, lcbn, 2pde, Ishl, lfca, .lptq, lbor, luxc, lvjw, 1fxd, 1hpt, 1mbg, lbpi, 1r69, 451e, 1neq, 1a2s, lsvr, 1frd, 1a63 and 1a7m. chart. These errors in the surface potential, usually as large as —3kcal/mol/e to 3kcal / mol / e, may contribute to large discrepancies in calculating the binding affinity and reaction rate of protein-ligand association. More detailed work is required to elucidate and justify the role of this discrepancy. 109 was 2.8 a: a; 8.8 53 8s a; as 8s 2s as 8%. Sam. was- was- was 3.2 88 new 8.2 8.8 as.” as was 8s 8s and 3s Raw- 3%. was- was- as 3.2 as: Sc moi 33s as a: s: was and was one 3%. swam- 8%. was- as 3mm 8% 8.: 3mm 8% 8.: 9% Ba was a: as was swam- mass- 2%. as- as s s a s s s s s s a. s s 8%. 8mm 2: 92 mm? Game 35 m3 o< s 110 EMU—cocoa cacti E Eotm .BEbefi wmAw- mm $.85 nose/Em 898 BC. .0396 3:: @2828 a £5 83% c. .8“ 35:32“ 83:5 2: 5 Ste 2: 95 BE>s§ E .04 38:0 :oEcZOm osfimosoflm an...“ 2an 111 w”: mods dde mmd :d mad dud dad and dmdd- mmfid- dvsd- Hd 2.3 mdAw 84cm mwd 26 mm: mvd mwd Ed add. 3.3. 8.3- ddvsd- Nd wd dwdm wndw mmdv dmd «v.5 5.3 mwd made 3.3 dwddT dde- dmdd- md dvdfi dmds dvd Hmd wad mad wad edd 3d dbdd- dvdw- vddd- Hd no.3 mod» 5.: mwd 2d mad dvd and SA wHAd- wvdw- «odd- dmddd- md dd mmdm mmdw dwdm add «6.2 mmd and ddgm 5.» Quad. wmdw- wadd- md Rafi 38d ddd add add ddd had Ed 3d 8.3- dsdw- cadm- Hd 8.3 dads 3.2 nvwd vmd dmd wmd mod Ed dwdw- Sawm- hmdw- mamdw- ad ad 3.6m mvdw mdsm mmd ddd ddé 3d 3.3 mad mwdw- mddw- gamm- md wddfl ammo 3d mmd mvd mdd ddd mvd odd mmdw- Edd- dde- dd 3.3 wmdb 3d mwd mm; and and dug med cudd- mmdw- dddd- mow.mw- Nd md madam mwdw no.8 mod mddfi ddd mad 3.: dds vwdw- mwdw- wwdw- md no no 5 no mm 5 mm mm 5 Gama 23 mag 58am Gama 2: E2 o4 s s E52303 832% E £02m .omamao dccofisosoc Edam a 53» @8ng a 8d Banged 833m 2: 5 note 95 database 5 Cd 3.85 2232?. US$35.55 Jim 2an. 112 dad dads dem mmd Hmd dmd mud dwd de ddSwT ENS- mwdwT Md dvdfi dfldw dmdm mmA No.3 vmdfi mwd dHNm ddém nmAdH- Sad:- SAdH- vmmdwT Nd ma dfidm dddw ddfim deH mddv mmdm 5.3 Hddm Sada ode-d7 mmAdm- 31me md Hmdfi dash. vwdm mmd dmd 2A Ed wwé 3: Eme NHdNH- stmfi- Hd 3.3 wwdw dmdv dwd dad vwé 3d mdd mwd 3.de- mmde- mag-NH- ddhsg- Nd mg mmdm dmdw wddv mos mdsm 8.3 and Edd 5me 8.2:. 3.954- eddmfi- dd 5.: mad» wwdfi mmd mmA hvd :d voA wmd :dz- HmddH- dedH- Hd Sim; Edd dwdm vwd mud Hdd dvd and and dddz- Hmddfi- :ddT Haddfi- md dA mmmm dmdw Ede had. dwdm Sand mod Hmdv vmdfi mde- dvfidfi- ddddH- md no we 5 no me 5 no we 5 Ommm 2: mdz Seam momma 2: m3 oq s s .csguoa 8.35m E maotm .owaazo Cocoucoococ Swim a fig 085% a cod 35:89“ 852% 2: E She can EEEQQ E 04 xucccc :cCcZCm ctecmcbcim "min @338 113 wmdfi Hddo swofi mad odd mad ddd ood dmd Sdom- omflom- odoom- Hd no.3 wads no.5 owd Hod so; mad wvd «mm oonm- wonm- odooo- eddoom- Nd wd 26m mood oddo woo 3.3 3.3 SMN oodm ode gonn- sooom- $.wom- od Eda Sumo vmdfi mad ood Ed ddd Hod wad woodm- smdvm- wmdvm- Hd wag; 2d» 8.3 vod 2mm odA odd odd wmé omfivm- endem- hmdwm- www.mvm- Nd od deem Hood womb mod 3.: mad. mod mhgm no.3 mwdom- Simeo- vddvm- od dddfi owoo Vddfi mud dod mod and vod :d mwomm- vmdmm- 36mm- dd 3.3 moon doAm odd odd Rd odd boa SA omdmm- wmdmm- modmm- duodmm- ad ad oosvm wood mvdv mvd nod mod on; dodH No.3 owdom- 5.3%- wo.vmm- od No.2 ofivo Ed mmd ovd sdd :d mvd de doAmm- mmdmm- vodmm- Hd omdfi mud» oddH mod on; dod omd mN-A whd oddmm- hmdmm- Edmo- ooodmm- Nd md soda stvw ods? wvd and dds dog oddfl dde modem- dwdmm- oodmm- od no no 5 no mo 5 no mo 5 Gmmm 2: mdz 3er 3mm 2: m2: o4 s s 35:80am media 5 £02m .mowceso B83820: 93 £33 29on a nod Raccoon 8.35m 23 E .555 2: one 3538.4 5 04 compose :ofleZOm 0338585 "do 23mm- 114 dnd n5: Ego and odqV NNo Hnd onN Hod uNdHo- nNddv- onoB- do 3.: no.5 dodo ooé 2.3 nnoH ddd nnnn oan dHoNo- 2.591 SoNo- nSoHo- Nd oa NHAN Sdn dddo No.3 NoSn 34» no.3 non: 3.3 dnAvo- oddvo- ondvo- od 3.: Nion ddAn dnd 2N nHA an ova wna vnnfiv- dNodn- odddv- Hd anH nnAn nod“V ndA doo n56 nod nn._.V ond ooog- Nnodn- onddfi ndnddv- Nd N; Nan nodn ddoo Hon :Nn dodo 3o oodm :An nosnv- SEN? oH.oNv- od 3.: 3N5 Nn.NN nNd ndd nNd de wad nod dndnn- unson- nnshn- Hd 2.3 nndn den Ndd one an 3d Non non $.Nnn- Soon- dvsnn- Snswn- Nd da NonN NN.on Ndso $6 on.: SoN Sn dodn oddH onsdn- duomn- oddnn- od nu Nm 5 mm mm 5 mm Nm 5 Ommm 2: BE 3me Gama 2: E2 64 a a i Eggnog moflSw E 20th .mowzfiu 8.83025: 93 its 223m a so ESQBOQ oofism oi E Hobo 2: 95 ochsné E Dd kahuna :cSaZOm 3:306:85 ”Nam 23m? Table 5.8: Electrostatic solvation energy in kcal/mol of two positively charged spheres at four different mesh sizes h = 0.5, 0.2, 0.1 and 0.05. The distance between the centers of two spheres D vary from 2A to 5A. D MIB PBEQ 0.5 0.2 0.1 0.05 0.5 0.2 0.1 0.05 2.0 276.39 276.78 276.83 276.84 283.12 279.66 277.70 277.10 2.5 266.68 -267.15 267.22 267.24 272.28 269.05 268.15 267.50 3.0 257.69 258.00 258.06 258.07 266.02 260.34 259.10 258.60 3.5 249.28 249.39 249.43 249.44 252.72 250.82 250.33 249.44 4.0 241.28 241.45 241.49 241.50 247.92 243.39 242.40 241.90 4.5 234.39 234.38 234.38 234.38 238.76 235.73 235.22 234.80 5.0 228.10 228.18 228.16 228.17 235.70 230.04 229.10 228.60 ewe 5m: :53 $me 833 $8: 38: 60$ 38V 8me :5 :3 A959 Boo €2.30me 3%. 3m:- OSQ- 3%- 92mm- oéwm- 22- gm. oki- Emw- c.25- 33- ea . at am: 83 a: and 32V 28.: as 33 E: 83 an eac- Dmov Qweeomma 3mm- ofiz- £3:- admw- 33”. 3mg- «.82- 22”- 322. 3.8- 32:- 337 o4 so a: a: a: a: a: a: A: 2: so E a: $8:- omov 2233mm 0.2a- 30m:- m.wa2- 33. @82- adam- 23.2- 3mm- @82- 33- $2:- 35- ea 33 S: 88 and :3 so: Gd 3i :8 :8 an an $83. EB 2353 3%- 3E. 232- SE- H.Smm- 3.3%. 282- 3%- 33. 3mm- 33. man:- we n3 :2: mom mew 3% we; 88 we mam Sm Sn 2... .3885 25 385:2 32 ~55 was“ 8:: 3: E: 8: 22 afl .53 .a: .5: 9 mom .mofim :88 ucouotmo do Gama cad maz its meow-$5 coSwZOm 053.6658? woosafioo 23 98 $2805 ammo o0 .Aa-mEEsm ”do 2an 116 $83 £83 35 GE 8me 38V 3% $.me cam-V 3:: 8va as: sag 50V 22.335 98:- :33 c.82- ogs- 3%. 3.3- 332- 38- 337 922. Sé- 9.82- uq .. Ame-V 59 $3 33 V g V av V 32V an V at as V $3 38 A259 50V 2333mm Saa- oéfi- c.32- OQS- 3mm- 33- @32- 33- 3a:- o.o§- 3:. m2:- uq $3 :3 :mV EV E 8V EV 8V 3V a: CV 3V 38:.- DmoV 2983mm 38m- 53m- 982- 3%. m.%m- woma- fimfi- $8- 3o:- mewz- Sow- 23H- sq 32V 38V 32V $8 83 3V 23V 33 a: 32V $3 QwV $8:- DmoV £3.8sz 32m- ofimm- 32:- gms- 98w- ©.m§- 35- 28- m3:- HQE- SE- Sme- sq 8% 38 ES 22 EV as. mam cam 8w £3 NE :3 was“ V0 3952 52 as $3 23 03“ EN 3.5 a: as: :2 22 mo: a mom .mmfim :88 €98me Hm Ummm was 52 fits mo_m$:m c2338 osfimoSooVo @85928 St was wars-ca “mop Vo 5&5:an "OH-m 035- 117 Figure 5.6: The map of electrostatic potential 011 the surface of protein IRS-3. Upper: potcntiul (:omputvd with MIB; Lower: l)iffcrc11cc of pot-cutiuls computed with PBEQ and with MIB. 118 Chapter 6 Thesis Contribution and Future Research Plan 6.1 Thesis Contribution 6.1.1 The Matched Interface and Boundary Method A new approach for solving elliptic interface problems, the matched interface and boundary (MIB) method, is generalized from its primitive 1-D formulation to higher dimensions in this dissertation. The MIB method, after generalization, is capable to solve 2—D or 3-D elliptic problems with discontinuous diffusion coefficients or singular sources at irregular internal interfaces which are usually not aligned with the mesh lines. The essential ingredient of the MIB method is the matching of the lowest order interface conditions simultaneously and exactly at select points on the interface with fictitious values. T hcse fictitious values represent the extension of the solutions in its subdomain to other subdomains. For an interface that is not aligned with the mesh lines, appropriate coordinate transformations are defined such that the interface conditions in the normal direction can be handled in the Cartesian grid. These interface conditions are re-organized such that there is a primary direction in the final set of interface conditions. All the fictitious values are solved along its primary direction but can be used for the discretization of partial derivatives with 119 respect to any coordinate. The matching of interface conditions provides essential algebraic relations for the expansion of these fictitious values in terms of real solution values and given jumps. Once the fictitious values are solved, in other words, once each subdomain is extended to a sufficiently large margin, the standard central finite difference scheme can be applied on the extended smooth domain to obtain desired high convergence and high accuracy. The same set of interface conditions can be iteratively used to solve for more fictitious values to support high order difference schemes. Among a variety of numerical methods for elliptic interface problems such as immersed interface method (IIM), immersed boundary method and ghost fluid method, MIB is the first numerical method that can be systmatically generalized to arbitrary high order for a general interface. The fourth-order and a sixth-order MIB schemes are successfully formulated in this work. An interpolation formulation is also proposed in this dissertation whereby the in- terface conditions are matched via explicitly defined polynomials instead of through fictitious values. The comparison between these two implementations describes the similarities and relations among different interface methods based on Taylor expan- sion(IIM), Lagrange interpolation(MIB) and explicit polynomials(interpolation for- mulation). 6.1.2 Interface Methods for the Poisson-Boltzmann equation The application of interface methods can be found in a variety of disciplines such as fluid dynamics, material science, electromagnetics and molecular biology. However, this dissertation describes the first application of interface methods to the numerical solution of the Poisson-Boltzmann equation to conserve the continuity of the electro- static potential and its flux at the dielectric interface. It is noticed that in most finite difference Poisson-Boltzmann solvers the molecular surface is only used to determine the accessibility of grid points but the exact location of the molecular surface is not considered. Therefore it was impossible for these methods to maintain the continuity of the potential flux at the molecular interface. In contrast, the analytical molecular 120 surface generated by MSMS is incorporated into MIB for the discretization of the equation near the molecular surface. The accuracy of the solution, measured in ei- ther electrostatic potential field or the derived quantity such as electrostatic solvation energy, is significantly improved due to the consideration of the discontinuous poten- tial gradient at the molecular surface. The accurate surface potential, solved with interface methods, makes it possible to quantitatively investigate the electrostatic complementary of the molecular surfaces involved in many bimolecular interactions such as electrostatically directed ligand association. The convergent electrostatic po— tential also provides an accurate calibration for the Born radii in the Generalized Born methods. 6.2 Future research plan 6.2.1 Theoretical analysis and fast algebraic solvers for the MIB method The matched interface and boundary method discussed in this dissertation is proved to be very promising for the interface problems, as illustrated by our extensive numerical experiments. Nevertheless, there are many directions for extending this thesis work. Firstly, we expect to develop a solid analysis of the convergence of the MIB method. The approximation of the interface conditions in the current MIB method implies that the local truncation errors at the irregular points are one order lower than that at the regular points. Although it is commonly believed and actually demon— strated in this thesis that the convergence rate of the global error is not affected, we want to find out the rigorous conditions under which this claim is true. Secondly, there are several numerical issues associated with the MIB method. The coefficient matrix of the linear system for the MIB method is always not symmetric and might not be diagonally dominant. Techniques to solve such systems efficiently are important for the applications of MIB method. The convergence might become a major problem with the increasing of the percentage of numbers of irregular points. 121 It is expected that some of the state-of-the—art techniques such as the multi—grid method, GMRES, QMR, etc, could by customized for the fast solution of the linear system derived from our MIB method. Thirdly, the MIB method formulated in this thesis is based on the finite difference on a Cartesian grid. However, the essential idea of MIB, the matching of the interface conditions with fictitious values, might be extended to other formulations like varia- tional principles or finite element methods on adaptive and / or composite meshes to attain more flexibility in handling the complex boundary or internal interface. Finally, there are a large number of interface problems in which the interface is moving. It is therefore worthwhile to combine our MIB with the level set or similar marching methods for the-moving interfaces to provide an accurate solutions for this type of problems. 6.2.2 Investigation of biological significance of accurate sur- face potential In this thesis, the MIB method was used to solve for the convergent and accurate electrostatic potential from the Poisson-Boltzmann equation. The significance of the MIB method, however, is mostly substantiated by the accurate electrostatic solva— tion energies and the implication of an accurate surface potential in the analysis of the protein/ligand binding and the steering of ligands toward the protein. It is these applications that would signify the applications of the interface methods for the Poisson-Boltzmann equation. It is necessary to quantify the surface potential in analysis of the binding or association. The physical model and mathematical tools for such quantitative analysis, however, have not yet been well established. Most studies were based on the binding free energy or association entropy and molar partition coefficient[5]. Only a small number of works were based on the similarity analysis [82] on the electrostatic potential of different enzymes and the scoring functions using geometric-electrostatic correlation functions [26]. We hope that the accurate surface potential solved with the MIB method can be helpful in the accurate identification 122 and prediction of binding site and binding affinity with these new analytical tools and protocols. 1‘23 ...—. . .... Zia—1.1. a.‘......, .. f S E m D N E P P A .l. . ’2‘“ 124 "Mi. APPENDIX A Kirkwood Solution to Non-Centered Charged Sphere We consider a sphere of radius b in a solvent of dielectric constant cw. The sphere has dielectric constant 6p and M discrete point charges q1,q2, . . . ,qM inside. The sphere-solvent complex is schematically show below. Following the general solution due to Kirkwood [39], one can choose a spherical coordinate system with origin at the center of the sphere. The positions of M charges are given by (rk, 9k, (pk), k = 1, . . . , M. It should be noted that here 0 E [0, 7r] is the polar angle from z—axis and 90 E [0, 27r) is the azimuthal angle in the :c -— y plane from the x—axis. This notational convention is different from that we used in the derivation of the interface scheme. It is adopted here to maintain the consistence with the original notation used by Kirkwood and the convention normally used in physics. The electrostatic potential inside the sphere gbp satisfies the Poisson equation and is given by a spherical harmonic series .M 00 n ¢p = Z—qL—i-Z Z Bn,mrnP,T,"'(cosd)eim‘p (A.1) k=1 Cplr _ Tkl 11:20 1712—71 00 Tl E 00 n Tm , . . . i p = Z 2 E r711+1P,’,"(cosfi)ezm*°+ E Z BnmrnP,T(cos6)e’m“’(,A.2) "Zomz‘" n=0m=—n where Ir — Tkl is the distance of the charge qk from the point (136,39) where the 125 electrostatic potential is to be measured. The first spherical harmonic represents the potential induced by the point charges (II: and the second is the dielectrical contri- butions of the solvent outside of the sphere. P,’,"‘(cos 0) are the associated Legendre functions of the first kind and i is the imaginary unit. Coefficients Enm are deter- mined by the distribution of the point charges, as (n - __|__m|)! : n 1n(—ima,ok Enm (71’ + ——|——"n1|)! '2 rik (COS 6k)e (A.3) Assume the ionic concentration is zero in the solvent, the electrostatic potential outside of the sphere (bu- shall satisfy the Laplace equation and have a spherical harmonic expansion OO Tl ‘ C , (Du! = E : E : (MT; +(1nm7-")1,§" (cosfl)e "W (AA) All the coefficients, 13,-,r,,,,(7n,n and (Inm, have to be fixed by fitting the appropriate boundary conditions. As potential aw must be vanishing as r approaches to 00, Gnm = 0 follows. The remaining coefficients, Bnm and Cnm, will be determined by fitting the continuity conditions on the surface of the sphere, i.e., ”QIL— Lap ‘ 7 = ' , f y—‘— — A.5 ¢u 95]) u 67‘ (p— (91‘ ( ) which provide two equations for every harmonic component of the potential: Enm 271+] —6 + b Bnm 2 Cum, (A6) P (n + 1)Enm —- nepb2n+anm = (n + Dru-Cnm. (A.7) 12G Solution of these two equations gives the expression of Bnm and Cnm: B = _ , A.8 ”m [nep + (n + 1)eu-]b2n+1 ( l C = . A.9 nm ncp + (n + 1)éw ( ) The final real-valued solution to electrostatic potential can be obtained by con- sidering only the real part of equations (A2) and (AA), although the coefficients Enm, Cnm, Bnm determined as above are complex functions. 1‘27 References [ll [2] [3] l4] [5] [6] [7] [8] M [10] [11] [12] L. Adams and TR Chartier. New geometric immersed interface multigrid solvers. SIAM Journal on Scientific Computing, 25:1516—1533, 2004. B. H. Ameur, M. Burger, and B. Hack]. Level set methods for geometric inverse problems in linear elasticity. Inverse Problems, 20:673—696, 2004. I. Babuska. The finite element method for elliptic equations with discontinuous coefficients. Computing, 5:207—213, 1970. NA. Baker. Biomolecular applications of Poisson-Boltzmann methods. John Wiley and Sons, 2005. ' ‘ N. Ben-Tal, B. Honig, C. Miller, and SJ McLaughlin. Electrostatic binding of proteins to membranes. theoretical predictions and experimental results with charybdotoxin and phospholipid vesicles. Biophysical Journal, 73(4):1717—1727, 1997. ' ' PA. Berthelsen. A decomposed immersed interface method for variable coeffi- cient elliptic equations with non-smooth and discontinuous solutions. Journal of Computational Physics, 197:364—386, 2004. G. Biros, L.X. Ying, and D. Zorin. A fast solver for the stokes equations with distributed forces in complex geometries. Journal of Computational Physics, 193:317—348, 2004. A. H. Boschitsch, M.O. Fenley, and H.X. Zhou. Fast boundary element method for the linear poisson-boltzmann equation. Journal of physical Chemistry B, 106(10):2741—2754, 2002. A. H. Boschitsch and O. Fenley. Hybrid boundary element and finite differ- ence method for solving the nonlinear poisson-boltzmann equation. Journal of Computational Chemistry, 25(7):935—955, 2004. J. Bramble and J. King. A finite element method for interface problems in domains with smooth boundaries and interfaces. Advances in Cmnputational Mathematics, 6:109—138, 1996. Z. Chen and J. Zou. Finite element methods and their convergence for elliptic and parabolic interface problems. Numerische Mathematik, 79:175-202, 1998. ML. Connolly. Molecular surfaces: A review. 1996. 128 [13] S.Z. Deng, K. Ito, and Z.L. Li. Three-dimensional elliptic solvers for inter- face problems and applications. Journal of Computational Physics, 184:215-243, 2003. [14] MA. Dumett and JP. Keener. An immersed interface method for solving anisotropic elliptic boundary value problems in three dimensions. SIAM Journal on Scientific Computing, 25:348—367, 2003. [15] J. Dzubiella, J.M.J. Swanson, and J.A. McCammon. Coupling hydrophobic- ity, dispersion, and electrostatics in continuum solvent models. Physical review letters, 96(8):087802+, 2006. [16] M. Eigen and E. Wicke. The thermodynamics of electrolytes at higher concer- tration. Journal of Chemical Physics, 58(9):702—714, 1954. [17] EA. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof. Combined immersed- boundary finite-difference methods for three-dimensional complex flow simula- tions. Journal of Computational Physics, 161:30—60, 2000. [18] RP. Fedkiw, T. Aslam, B. Merriman, and S. S. Osher. A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method). Journal of Computational Physics, 152:457—492, 1999. [19] M. Feig and CL. Brooks III. Recent advances in the development and applica- tion of implicit solvent models in biomolecule simulations. Current Opinions in Structural Biology, 14:217—224, 2004. [20] M. Feig, A. Onufriev, MS. Lee, W. Im, D.A. Case, and CL. Brooks III. Perfor- mance comparison of generalized born and poisson methods in the calculation of electrostatic solvation energies for protein stuctures. Journal of Computational Chemistry, 25:265-284, 2004. [21] J .Q. Feng. Electrostatic interaction between two charges dielectric spheres in contact. Physical Review E, 62(2):2891—2897, 2000. [22] B. Fornberg. Calculation of weights in finite difference formulas. SIAM Review, 40:685—691, 1998. [23] F. Cibou and RP. Fedkiw. A forth order accurate discretization for the laplace and heat equations on arbitrary domains, with applications to the stefan prob- lem. Journal of Computational Physics, 202:577—601, 2005. [24] M.K. Gilson and B. Honig. Calculation of the total electrostatic energy of a macromolecular system: Solvation energies, binding energies and conformational analysis. Proteins, 427—18, 1988. [25] M.K. Gilson, K.A. Sharp, and RH. Honig. Calculating the electrostatic potential of molecules in solution: method and error assessment. Journal of Computational chemistry, 9(4):327—335, 1987. [26] SC. Gupta. The classic Stefen problem. JAI Press, 2003. [27] A. Heifetz, E. Katehalski-Katzir, and M. Eisenstein. Electrostatics in protein- protein docking. Protein Science, 11(3), 2002. 129 [28] M. Holst, N. Baker, and P. Wang. Adaptive multilevel finite element solution of the poisson-boltzmann equation i. algorithms and examples. J Ournal of compu- tational chemistry, 21(15):1319—1342, 2000. [29] M.J. Holst. The Poisson-Boltzmann Equation: Analysis and Multilevel Numeri- cal Solution. PhD thesis, University of Illinois, Urbana-Champaign, 1993. [30] S. Hou and X.-D. Liu. A numerical method for solving variable coefficient elliptic equation with interfaces. Journal of Computational Physics, 202:411—445, 2005. [31] T.Y. Hou, Z.L. Li, S. Osher, and H. Zhao. A hybrid method for moving inter- face problems with application to the hele-shaw flow. Journal of Computational Physics, 134:236—252, 1997. [32] H. Huang and Z.L. Li. Convergence analysis of the immersed interface method. IMA Journal of Numerical Analysis, 19:583—608, 1999. [33] J.K. Hunter, Z.L. Li, and H. Zhao. Reactive autophobic spreading of drops. Journal of Computational Physics, 183:335—366, 2002. [34] G. Iaccarino and R. Verzicco. Immersed boundary technique for turbulent flow simulations. Applied Mechnics Reviews, pages 331—347, 2003. [35] CC. Jenkins and K.L.C. Hunt. Nonlocal dielectric functions on the nanoscale: Screened forces from unscreened potentialsa. Journal of Chemical Physics, 119(16):8250—8256, 2003. [36] S. Jin and X.L. Wang. Robust numerical simulation of porosity evolution in chemical vapor infiltration ii. two-dimensional anisotropic fronts. Journal of Computational Physics, 179:557—577, 2002. [37] H. Johansen and P. Colella. A cartesian grid embedding boundary method for poisson’s equation on irregular domains. Journal of Computational Physics, 147:60—85, 1998. [38] L. Joly, C. Ybert, E. Trizac, and L. Bocquet. Hydrodynamics within the electric double layer on slipping surfaces. Physical review letters, 932257805, 2004. [39] JD. Kandilarov. Immersed interface method for a reaction-diffusion equation with a moving own concentrated source, volume 2542, pages 506-513. Springer, 2003. [40] J. G. Kirkwood. Theory of solution of molecules containing widely separated charges with special application to zwitterions. Journal of Chemical Physics, 7:351—361, 1934. [41] L. L. Adams and Z.L. Li. The immersed interface/multigrid methods for interface problems. SIAM Journal on Scientific Computing, 24:463—479, 2002. [42] L. Lee and R.J. LeVeque. An immersed interface method for incompressible navier-stokes equations. SIAM Journal on Scientific Computing, 25:832—856, 2003. 130 [43] R.J. LeVeque and Z.L. Li. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM Journal on Numerical Analysis, 31:1019—1044, 1994. [44] Z.L. Li. A fast iterative algorithm for elliptic interface problems. SIAM Journal on Numerical Analysis, pages 230—254, 1998. [45] Z.L. Li. The immersed interface method using a finite element formulation. Applied Numerical Mathematics, 27(3):253—267, 1998. [46] Z.L. Li. An overview of the immersed interface method and its applications. Taiwanese Journal of Mathematics, 7:1—49, 2003. [47] Z.L. Li and K. Ito. Maximum principle preserving schemes for interface problems with discontinuous coefficients. SIAM Journal on Scientific Computing, 23:339— 361, 2001. [48] Z.L. Li and K. Ito. The Immersed Interface Method — Numerical Solutions of PDEs Involving Interfaces and Irregular Domains. SIAM, 2006. [49] Z.L. Li and SR. Lubkin. Numerical analysis of interfacial two-dimensional stokes flow with discontinuous viscosity and variable surface tension. International Journal for Numerical Methods in Fluids, 37:525—540, 2001. [50] Z.L. Li, W.—C. Wang, I.-L. Chem, and M.-C. Lai. New formulations for interface problems in polar coordinates. SIAM Journal on Scientific Computing, 25:224— 245, 2003. [51] S. Lim and CS. Peskin. Simulations of the whirling instability by the immersed boundary method. SIAM Journal on Scientific Computing, 25:2066—2083, 2004. [52] X.D. Liu, Fedkiw RP, and M. Kang. A boundary condition capturing method for poisson’s equation on irregular domains. Journal of Computational Physics, 160:151—178, 2000. [03] B. Lombard and J. J. Piraux. How to incorporate the spring-mass conditions in finite-difference schemes. SIAM Journal on Scientific Computing, 24:1379—1407, 2003. [54] B.Z. Lu, D.Q. Zhang, and J .A. McCannnon. Computation of electrostatic forces between solvated molecules determined by the poisson-boltzmann equation using a boundary element method. Journal of Chemical Physics, 122(21):214102+, 2005. [55] G. Luo, S. Malkova, J. Yoon, C. Schultz, B. Lin, M. Meron, L. Benjamin, P. Vanysek, and M. L. Schlossman. Ion distributions near a liquid-liquid in- terface. Science, 311(5758):216—218, 2006. [56] R. Luo, H. David, Land Hung, J. Devaney, and MK. Gilson. Strength of solvent- exposed salt-bridges. Journal of Physical Chemistry: B, 103:727—736, 1999. [57] GS. Manning. Limiting laws and counterion condensation in polyelectrolyte solutions v. further development of the chemical model. Biophysical Chemistry, 9(1y65—60,1978. 131 [58] XL. Max and E. D. Getzoff. Spherical harmonic molecular surfaces. IEEE Computer Crphics, 8(4):42—50, 1988. [59] A. Mayo. The fast solution of poisson’s and the biharmonic equations on irregular regions. SIAM Journal on Numerical Analysis, 21:285—299, 1984. [60] A. Mckenney, L. Greengard, , and A Mayo. A fast poisson solver for complex geometries. Journal of Computational Physics, 118:348—355, 1995. [61] H. Ohshima, E. Mishonova, and E. Alexov. Electrostatic interaction between two charged spherical molecules. Biophysical Chemistry, 57(2):189—203, 1996. [62] W.H. Orttung. Direct solution of the poisson equation for biomolecules of arbi- trary shape, polarizability density and charge distribution. Annals of the New York Academy of Sciences, 303:22—37, 1977. [63] C .S. Peskin. Numerical analysis of blood flow in heart. Journal of Computational Physics, 25:220—252, 1977. [64] CS. Peskin. Lectures on mathematical aspects of physiology. Lectures in Applied Mathematics, 19:69—107, 1981. [65] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and HP. F lannery. Numerical recipes in Fortran: the art of scientific computing. Cambridge University Press, 1999. [66] AA. Rashin. Hydration phenomena, classical electrostatics, and the boundary element method. Journal of physical chemistry, 94(5):1725—1733, 1990. [67] VA. Rokerts, H.C Freeman, A.J. Olson, J.A. Tainer, and ED. Getzoff. Elec- trostatic orientiation of the electron-transfer complex between plastocyanin and cycochrome c*. Journal of Biological Chemistry, 266:13431—13441, 199L. [68] M.F. Sanner, A.J. Olson, and J.C. Spelmer. Fast and robust computation of molecular surface. In Proceedings of the 11th symposium on Computational ge- ometry, pages C6—C7, 1995. [69] M. Schulz and G. Steinebach. Two-dimensional modelling of the river rhine. Journal of Computational and Applied Mathematics, 145:11—20, 2002. [70] J.A. Sethian. Level Set Methods and Fast Marching Methods : Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science. Cambridge University Press, 1999. [71] J .A. Sethian and A. Wiegmann. Structural boundary design via level set and immersed interface methods. Journal of Computational Physics, 163, 2000. [72] K. Sharp, A. Jean-Charles, and B. Honig. Journla of Physical Chemistry, 96:3822—3828, 1992. [73] EB. Sheinerman, R. Norel, and B. Honig. Electrostatic aspects of proteinprotein interactions. Current Opinion in Structural Biology, 10:153—159, 2000. 132 [74] [76] [77] [78] l79l [80] [81l [82] [83l [84] [85] [86] LA. Shkel, O.V. Tsodikov, and M.T. Jr. Record. Asymptotic solution of the cylindrical nonlinear poisson-boltzmann equation at low salt concentration: an- alytical expressions for surface potential and preferential interaction coefficient. Proceedings of the National Academy of Sciences of the United States of America, 99(5):2597—2602, 2002. T. Simonson. Dielectric constant of cytochrome c from simulation in a water droplet including all electrostatic interactions. Journal of American Chemisty Society, 120:4875—4876, 1998. T. Simonson and D. Perahia. Proceedings of the National Academy of Sciences of the United States of America, 92:1082—1086, 1995. X. Song. An inhomogeneous model of protein dielectric properties: intrinsic polarizabilities of amino acids. Journal of Chemical Physics, 116(21):9359—9363, 2002. N.V. Sushkin and G.D.L. Phillies. Charged dielectric spheres in electrolyte so- lutions: Induced dipole and counterion exclusion effects. Journal of Chemical Physics, 103(11):4600—4612, 1995. AK. Tornberg and B. Engquist. Numerical approximations of singular source terms in differential equations. Journal of Computational Physics, 200:462—488, 2004. L. Visscher and K. G. Dyall. Dirac-fock atomic electronic structure calculations using different nuclear charge distributions. Atmoic data and nuclear data tables, 67:207—224, 1997. J .V. Voorde, J. Vierendeels, and E. Dick. Flow simulations in rotary volumetric pumps and compressors with the fictitious domain method. Journal of Compu- tational and Applied Mathematics, 168:491—499, 2004. RC. Wade, R.R. Gabdoulline, S.K. Liidemann, and V. Lounnas. Eletrostatic streering and ionic techering in enzyme-ligand binding: insights from simulations. Proceedings of the National Academy of Sciences of the United States of America, 95(11):5942—5949, 1998. J. H. Walther and G. Morgenthal. An immersed interface method for the vortex- in-cell algorithm. Journal of Turbulence, 3(39):1—9, 2002. W.-C. Wang. A jump condition capturing finite difference scheme for elliptic interface problems. SIAM Journal on Scientific Computing, 25:1479—1496, 2004. A. Warshel and A. Papazyan. Electrostatic effects in macromolecules: funda- mental concepts and pratical modeling. Corrent Opinion in structural biology, 8:211—217, 1998. J. \Varwicker and H.C. Watson. Calculation of the electric potential in the active site cleft due to a-helix dipoles. Journal of Molecular Biology, 157(4):671~679, 1982. 133 [87] G.W’. Wei, Y.H. Sun, Y.C. Zhou, and M. Feig. Molecular multiresolution sur- faces. 2005. arXivzmath-ph/0511001. [88] J. Weiser, P.S. Shenkin, and WC. Still. Optimization of gaussian surface calcula- tions and extension to solvent-accessible surface area. Journal of Computational chemistry, 20(7):688-703, 1999. [89] A. Wiegmann and KP. Bube. The immersed interface method for nonlinear differential equations with discontinuous coefficients and singular sources. SIAM Journal on Numerical Analysis, 35:177—200, 1998. [90] A. Wiegmann and KP. Bube. The explicit-jump immersed interface method: finite difference methods for pdes with piecewise smooth solutions. SIAM Journal on Numerical Analysis, 37:827-862, 2000. [91] J .C. Xu. Error estimates of the finite element method for the 2nd order elliptic equation with discontinuous coefficients. Journal of Xiangtan University, 1, 1982. [92] Y.-H. Tseng Y.-H. and J .H. Ferziger. A ghost-cell immersed boundary method for flow in complex geometry. Journal of Computational Physics, 192:593—623, 2003. [93] R.J. Zauhar and RS. Morgan. A new method for computing the macromolecular electric potential. Journal of Molecular Biology, 186(4):815—820, 1985. [94] S. Zhao, , and G.W. Wei. High order fdtd methods via derivative matching for maxwell’s equations with material interfaces. Journal of Computational Physics, 200:60—103, 2004. [95] R. Zhou and B.J. Berne. Can a continuum solvent model reproduce the free energy landscape of a B-hairpin folding in water? Proceedings of the National Academy of Sciences of the United States of America, 99:12777—12782, 2002. [96] Y.C. Zhou, S. Zhao, M. Feig, and G.W. Wei. High order matched interface and boundary (mib) schemes for elliptic equations with discontinuous coefficients and singular sources. Journal of Computational Physics, 213:1-30, 2006. 134