Additive Schwarz method

Last updated

In mathematics, the additive Schwarz method, named after Hermann Schwarz, solves a boundary value problem for a partial differential equation approximately by splitting it into boundary value problems on smaller domains and adding the results.

Contents

Overview

Partial differential equations (PDEs) are used in all sciences to model phenomena. For the purpose of exposition, we give an example physical problem and the accompanying boundary value problem (BVP). Even if the reader is unfamiliar with the notation, the purpose is merely to show what a BVP looks like when written down.

(Model problem) The heat distribution in a square metal plate such that the left edge is kept at 1 degree, and the other edges are kept at 0 degree, after letting it sit for a long period of time satisfies the following boundary value problem:
fxx(x,y) + fyy(x,y) = 0
f(0,y) = 1; f(x,0) = f(x,1) = f(1,y) = 0
where f is the unknown function, fxx and fyy denote the second partial derivatives with respect to x and y, respectively.

Here, the domain is the square [0,1] × [0,1].

This particular problem can be solved exactly on paper, so there is no need for a computer. However, this is an exceptional case, and most BVPs cannot be solved exactly. The only possibility is to use a computer to find an approximate solution.

Solving on a computer

A typical way of doing this is to samplef at regular intervals in the square [0,1] × [0,1]. For instance, we could take 8 samples in the x direction at x = 0.1, 0.2, ..., 0.8 and 0.9, and 8 samples in the y direction at similar coordinates. We would then have 64 samples of the square, at places like (0.2,0.8) and (0.6,0.6). The goal of the computer program would be to calculate the value of f at those 64 points, which seems easier than finding an abstract function of the square.

There are some difficulties, for instance it is not possible to calculate fxx(0.5,0.5) knowing f at only 64 points in the square. To overcome this, one uses some sort of numerical approximation of the derivatives, see for instance the finite element method or finite differences. We ignore these difficulties and concentrate on another aspect of the problem.

Solving linear problems

Whichever method we choose to solve this problem, we will need to solve a large linear system of equations. The reader may recall linear systems of equations from high school, they look like this:

2a + 5b = 12 (*)
6a − 3b = −3

This is a system of 2 equations in 2 unknowns (a and b). If we solve the BVP above in the manner suggested, we will need to solve a system of 64 equations in 64 unknowns. This is not a hard problem for modern computers, but if we use a larger number of samples, even modern computers cannot solve the BVP very efficiently.

Domain decomposition

Which brings us to domain decomposition methods. If we split the domain [0,1] × [0,1] into two subdomains [0,0.5] × [0,1] and [0.5,1] × [0,1], each has only half of the sample points. So we can try to solve a version of our model problem on each subdomain, but this time each subdomain has only 32 sample points. Finally, given the solutions on each subdomain, we can attempt to reconcile them to obtain a solution of the original problem on [0,1] × [0,1].

Size of the problems

In terms of the linear systems, we're trying to split the system of 64 equations in 64 unknowns into two systems of 32 equations in 32 unknowns. This would be a clear gain, for the following reason. Looking back at system (*), we see that there are 6 important pieces of information. They are the coefficients of a and b (2,5 on the first line and 6,−3 on the second line), and the right hand side (which we write as 12,−3). On the other hand, if we take two "systems" of 1 equation in 1 unknown, it might look like this:

System 1: 2a = 12
System 2: -3b = −3

We see that this system has only 4 important pieces of information. This means that a computer program will have an easier time solving two 1×1 systems than solving a single 2×2 system, because the pair of 1×1 systems are simpler than the single 2×2 system. While the 64×64 and 32×32 systems are too large to illustrate here, we could say by analogy that the 64×64 system has 4160 pieces of information, while the 32×32 systems each have 1056, or roughly a quarter of the 64×64 system.

Domain decomposition algorithm

Unfortunately, for technical reasons it is usually not possible to split our grid of 64 points (a 64×64 system of linear equations) into two grids of 32 points (two 32×32 systems of linear equations) and obtain an answer to the 64×64 system. Instead, the following algorithm is what actually happens:

1) Begin with an approximate solution of the 64×64 system.
2) From the 64×64 system, create two 32×32 systems to improve the approximate solution.
3) Solve the two 32×32 systems.
4) Put the two 32×32 solutions "together" to improve the approximate solution to the 64×64 system.
5) If the solution isn't very good yet, repeat from 2.

There are two ways in which this can be better than solving the base 64×64 system. First, if the number of repetitions of the algorithm is small, solving two 32×32 systems may be more efficient than solving a 64×64 system. Second, the two 32×32 systems need not be solved on the same computer, so this algorithm can be run in parallel to use the power of multiple computers.

In fact, solving two 32×32 systems instead of a 64×64 system on a single computer (without using parallelism) is unlikely to be efficient. However, if we use more than two subdomains, the picture can change. For instance, we could use four 16×16 problems, and there's a chance that solving these will be better than solving a single 64×64 problem even if the domain decomposition algorithm needs to iterate a few times.

A technical example

Here we assume that the reader is familiar with partial differential equations.

We will be solving the partial differential equation

uxx + uyy = f (**)

We impose boundedness at infinity.

We decompose the domain R² into two overlapping subdomains H1 = (− ∞,1] × R and H2 = [0,+ ∞) × R. In each subdomain, we will be solving a BVP of the form:

u( j )xx + u( j )yy = f in Hj
u( j )(xj,y) = g(y)

where x1 = 1 and x2 = 0 and taking boundedness at infinity as the other boundary condition. We denote the solution u( j ) of the above problem by S(f,g). Note that S is bilinear.

The Schwarz algorithm proceeds as follows:

  1. Start with approximate solutions u( 1 )0 and u( 2 )0 of the PDE in subdomains H1 and H2 respectively. Initialize k to 1.
  2. Calculate u( j )k + 1 = S(f,u(3 − j)k(xj)) with j = 1,2.
  3. Increase k by one and repeat 2 until sufficient precision is achieved.

See also

Related Research Articles

<span class="mw-page-title-main">Numerical analysis</span> Methods for numerical approximations

Numerical analysis is the study of algorithms that use numerical approximation for the problems of mathematical analysis. It is the study of numerical methods that attempt to find approximate solutions of problems rather than the exact ones. Numerical analysis finds application in all fields of engineering and the physical sciences, and in the 21st century also the life and social sciences like economics, medicine, business and even the arts. Current growth in computing power has enabled the use of more complex numerical analysis, providing detailed and realistic mathematical models in science and engineering. Examples of numerical analysis include: ordinary differential equations as found in celestial mechanics, numerical linear algebra in data analysis, and stochastic differential equations and Markov chains for simulating living cells in medicine and biology.

<span class="mw-page-title-main">Partial differential equation</span> Type of differential equation

In mathematics, a partial differential equation (PDE) is an equation which computes a function between various partial derivatives of a multivariable function.

Spectral methods are a class of techniques used in applied mathematics and scientific computing to numerically solve certain differential equations. The idea is to write the solution of the differential equation as a sum of certain "basis functions" and then to choose the coefficients in the sum in order to satisfy the differential equation as well as possible.

<span class="mw-page-title-main">Numerical methods for ordinary differential equations</span> Methods used to find numerical solutions of ordinary differential equations

Numerical methods for ordinary differential equations are methods used to find numerical approximations to the solutions of ordinary differential equations (ODEs). Their use is also known as "numerical integration", although this term can also refer to the computation of integrals.

In mathematics, a differential equation is an equation that relates one or more unknown functions and their derivatives. In applications, the functions generally represent physical quantities, the derivatives represent their rates of change, and the differential equation defines a relationship between the two. Such relations are common; therefore, differential equations play a prominent role in many disciplines including engineering, physics, economics, and biology.

Numerical methods for partial differential equations is the branch of numerical analysis that studies the numerical solution of partial differential equations (PDEs).

In numerical analysis, a multigrid method is an algorithm for solving differential equations using a hierarchy of discretizations. They are an example of a class of techniques called multiresolution methods, very useful in problems exhibiting multiple scales of behavior. For example, many basic relaxation methods exhibit different rates of convergence for short- and long-wavelength components, suggesting these different scales be treated differently, as in a Fourier analysis approach to multigrid. MG methods can be used as solvers as well as preconditioners.

An eikonal equation is a non-linear first-order partial differential equation that is encountered in problems of wave propagation.

Numerical linear algebra, sometimes called applied linear algebra, is the study of how matrix operations can be used to create computer algorithms which efficiently and accurately provide approximate answers to questions in continuous mathematics. It is a subfield of numerical analysis, and a type of linear algebra. Computers use floating-point arithmetic and cannot exactly represent irrational data, so when a computer algorithm is applied to a matrix of data, it can sometimes increase the difference between a number stored in the computer and the true number that it is an approximation of. Numerical linear algebra uses properties of vectors and matrices to develop computer algorithms that minimize the error introduced by the computer, and is also concerned with ensuring that the algorithm is as efficient as possible.

<span class="mw-page-title-main">Elliptic boundary value problem</span>

In mathematics, an elliptic boundary value problem is a special kind of boundary value problem which can be thought of as the steady state of an evolution problem. For example, the Dirichlet problem for the Laplacian gives the eventual distribution of heat in a room several hours after the heating is turned on.

<span class="mw-page-title-main">Domain decomposition methods</span>

In mathematics, numerical analysis, and numerical partial differential equations, domain decomposition methods solve a boundary value problem by splitting it into smaller boundary value problems on subdomains and iterating to coordinate the solution between adjacent subdomains. A coarse problem with one or few unknowns per subdomain is used to further coordinate the solution between the subdomains globally. The problems on the subdomains are independent, which makes domain decomposition methods suitable for parallel computing. Domain decomposition methods are typically used as preconditioners for Krylov space iterative methods, such as the conjugate gradient method, GMRES, and LOBPCG.

In numerical analysis, the balancing domain decomposition method (BDD) is an iterative method to find the solution of a symmetric positive definite system of linear algebraic equations arising from the finite element method. In each iteration, it combines the solution of local problems on non-overlapping subdomains with a coarse problem created from the subdomain nullspaces. BDD requires only solution of subdomain problems rather than access to the matrices of those problems, so it is applicable to situations where only the solution operators are available, such as in oil reservoir simulation by mixed finite elements. In its original formulation, BDD performs well only for 2nd order problems, such elasticity in 2D and 3D. For 4th order problems, such as plate bending, it needs to be modified by adding to the coarse problem special basis functions that enforce continuity of the solution at subdomain corners, which makes it however more expensive. The BDDC method uses the same corner basis functions as, but in an additive rather than multiplicative fashion. The dual counterpart to BDD is FETI, which enforces the equality of the solution between the subdomain by Lagrange multipliers. The base versions of BDD and FETI are not mathematically equivalent, though a special version of FETI designed to be robust for hard problems has the same eigenvalues and thus essentially the same performance as BDD.

In mathematics, Neumann–Neumann methods are domain decomposition preconditioners named so because they solve a Neumann problem on each subdomain on both sides of the interface between the subdomains. Just like all domain decomposition methods, so that the number of iterations does not grow with the number of subdomains, Neumann–Neumann methods require the solution of a coarse problem to provide global communication. The balancing domain decomposition is a Neumann–Neumann method with a special kind of coarse problem.

In numerical analysis, the Schur complement method, named after Issai Schur, is the basic and the earliest version of non-overlapping domain decomposition method, also called iterative substructuring. A finite element problem is split into non-overlapping subdomains, and the unknowns in the interiors of the subdomains are eliminated. The remaining Schur complement system on the unknowns associated with subdomain interfaces is solved by the conjugate gradient method.

<span class="mw-page-title-main">Finite element method</span> Numerical method for solving physical or engineering problems

The finite element method (FEM) is a popular method for numerically solving differential equations arising in engineering and mathematical modeling. Typical problem areas of interest include the traditional fields of structural analysis, heat transfer, fluid flow, mass transport, and electromagnetic potential. Computers are usually used to perform the calculations required. With high-speed supercomputers, better solutions can be achieved, and are often required to solve the largest and most complex problems.

Fluid motion is governed by the Navier–Stokes equations, a set of coupled and nonlinear partial differential equations derived from the basic laws of conservation of mass, momentum and energy. The unknowns are usually the flow velocity, the pressure and density and temperature. The analytical solution of this equation is impossible hence scientists resort to laboratory experiments in such situations. The answers delivered are, however, usually qualitatively different since dynamical and geometric similitude are difficult to enforce simultaneously between the lab experiment and the prototype. Furthermore, the design and construction of these experiments can be difficult, particularly for stratified rotating flows. Computational fluid dynamics (CFD) is an additional tool in the arsenal of scientists. In its early days CFD was often controversial, as it involved additional approximation to the governing equations and raised additional (legitimate) issues. Nowadays CFD is an established discipline alongside theoretical and experimental methods. This position is in large part due to the exponential growth of computer power which has allowed us to tackle ever larger and more complex problems.

In mathematics, the walk-on-spheres method (WoS) is a numerical probabilistic algorithm, or Monte-Carlo method, used mainly in order to approximate the solutions of some specific boundary value problem for partial differential equations (PDEs). The WoS method was first introduced by Mervin E. Muller in 1956 to solve Laplace's equation, and was since then generalized to other problems.

The proper generalized decomposition (PGD) is an iterative numerical method for solving boundary value problems (BVPs), that is, partial differential equations constrained by a set of boundary conditions, such as the Poisson's equation or the Laplace's equation.

<span class="mw-page-title-main">Physics-informed neural networks</span> Technique to solve partial differential equations

Physics-informed neural networks (PINNs), also referred to as Theory-Trained Neural Networks (TTNs), are a type of universal function approximators that can embed the knowledge of any physical laws that govern a given data-set in the learning process, and can be described by partial differential equations (PDEs). They overcome the low data availability of some biological and engineering systems that makes most state-of-the-art machine learning techniques lack robustness, rendering them ineffective in these scenarios. The prior knowledge of general physical laws acts in the training of neural networks (NNs) as a regularization agent that limits the space of admissible solutions, increasing the correctness of the function approximation. This way, embedding this prior information into a neural network results in enhancing the information content of the available data, facilitating the learning algorithm to capture the right solution and to generalize well even with a low amount of training examples.

References