One Hundred Years of the Galerkin Method

1 Galerkin Method

The mathematical idea inherent in the Galerkin method turned out to be extremely fruitful for creating the foundations of both qualitative and quantitative analysis of partial differential equations (PDEs), which began to be intensively studied since the middle of the 19th century. They proved to be much more difficult than ordinary differential equations that were studied earlier. Mathematical models based on PDEs arose from the problems of mechanics and physics. Of course, the question

How to find an approximate solution?

was one of the first.

Already in the late 19th century, it became clear that approximations of the problem ℒ⁢u=f (associated with a differential operator ℒ acting from a Banach space X to a Banach space Y) should be sought in the form

where N is a natural number and the coordinate functions wi∈XN⊂X form a linearly independent system. It is clear that the weights αi should be chosen from the condition that the residual ℒ⁢u-f is minimal (in some sense). At that time, the equations were studied in the framework of classical analysis and, therefore, the space X was considered as the space of continuous functions having sufficient amount of classical derivatives. It is not surprising that one of the first methods that appeared was the collocation method based on the exact satisfaction of the equation at some selected points.

However, in this way it was difficult to create a unified and mathematically justified approach, which required answers to the principal questions:

How to choose the coordinate and test functions and what form of the residual should be used?

Without having the right answers to these questions, it is impossible to prove that uN converges to the exact solution u.

Apparently the very first step in the right direction was made by Walter Ritz (1878–1909) [30] who suggested to define αi by minimization of the variational (energy) functional. The next step was taken by Boris Galerkin (1871–1945) who proposed a much more general approach suitable for partial differential equations of all types. Briefly, the essence of his approach can be formulated as follows:

The approximation is determined from the condition that the residual is orthogonal to a system of linearly independent test functions and the orthogonality relation is understood in the integral sense.

The first mathematical formulation of this method appeared in the paper [14] devoted to a biharmonic type equation (originally, it was published in Russian at the very end of 1915).

Figure 1

Boris Galerkin and the page in [14] (English translation) containing the main integral relation.

Boris Galerkin was born in Plotsk (nowadays this town is a part of Belarus). Later he moved to St. Petersburg and graduated from St. Petersburg Technological Institute. In 1909, he began to work in higher educational institutions of St. Petersburg (later Leningrad) among which we first of all mention Peter the Great Polytechnic Institute. Yet, as extraordinary as mathematical significance of Galerkin’s ideas, he mainly devoted his creative capacities to mechanics. The invention of the method was motivated by the analysis of mechanical constructions. Galerkin used an earlier publication of Ivan Bubnov (marine officer and engineer), where the method was applied to ship building problems.

In modern notation, the idea of the Galerkin method consists of finding αi by the orthogonality relation

where uN is defined by (1), the set YM⊂Y′ contains linearly independent test functions in the space conjugate to Y, and 〈⋅,⋅〉 denotes the duality pairing of Y and Y′. In general, the sets YM and XN are different (but of course they must be selected such that the system (2) is solvable).

This approach is very flexible. It includes the absolute majority of known methods. If X and Y are Hilbert spaces and the finite dimensional spaces XN and YM coincide, then the method is called Bubnov–Galerkin. This name was used by S. Mikhlin [26] who devoted a lot of time to studying the method and was the first to prove its convergence. If (2) originates from the Euler equation generated by a quadratic functional, then we arrive at the Ritz–Galerkin method. More general schemes with different spaces for coordinate and test functions were studied by G. Petrov who also extended the method to eigenvalue problems. The least squares method and the Fourier method can be also considered as particular forms of (2). The reader can find a systematic consideration of the Galerkin method in application to elliptic, parabolic, and hyperbolic equations in [35].

Figure 2

W. Ritz, I. Bubnov, and G. Petrov.

2 Galerkin Method and Correctness of Boundary Value Problems

The development of the theory of differential equations confirmed the correctness of the path outlined in the works of Galerkin and his successors, which rendered breakthrough results in qualitative and quantitative analysis of PDEs.

At the end of the 19th century it was unclear how to guarantee existence of solutions to boundary value problems for partial differential equations. This question was open even for the simplest boundary value problems. In particular, solvability of the problem

was an open problem intensively discussed in the literature (where it was typically understood in the classical sense, i.e., u was considered in C2⁢(Ω)∩C⁢(Ω¯)). This question was answered after many years of studies that reconstructed the theory of partial differential equations. A new conception of generalized or weak solutions was created by D. Hilbert, H. Poincaré, S. Sobolev, R. Courant, O. A. Ladyzhenskaya, and many other outstanding mathematicians.

In fact, the Galerkin method served as a turning point in the study of this problem. This is easy to see with the paradigm of (3) where the orthogonality relation (2) yields

and XN contains functions vanishing on the boundary ∂⁡Ω. Let us tend N to +∞. If we await that uN will tend (in some sense) to the exact solution, then it is natural to consider the limit form of (4), where XN is replaced by the whole space X. This way leads to the generalized statement of (2): find u∈X satisfying the boundary conditions with square summable ∇⁡u such that

where X is a proper closure of {XN}.

Integral type definitions of solutions to PDEs appeared in works of many mathematicians. A complete and systematic consideration of this approach is presented in [23, 24, 31]. It is worth noticing that integral type definitions of weak solutions are rather diverse. For example, in 1926, N. Wiener [34] suggested that “admissible” solutions of a linear equation ℒ⁢u=f should be considered among the functions satisfying

for any sufficiently smooth function w with a compact support in Ω, where ℒ* denotes the conjugate operator. This form of the integral identity is close to integral relations used in the so-called weak Galerkin method.

It is easy to show that a sequence of Galerkin solutions converges to a function which satisfies (4) provided that the subspaces {Xk} are limit dense in a reflexive Banach space X (i.e., for any w∈X, there exists a sequence wk∈Xk such that wk→w in X). Indeed, from (4) it follows that ∥∇⁡uN∥Ω≤C⁢∥f∥Ω with a constant independent of N. Hence, there exists a subsequence of uN weakly converging to a function u∈X. In view of the limit density property, for any w∈X there exists a sequence wk∈X that converges to w in V. Passing to the limit in (4), we conclude that u satisfies (5). We see that Galerkin approximations suggest a simple and efficient method for proving existence of weak solutions to various boundary value problems. E. Hopf used this idea in his famous paper [18] in order to prove existence of weak solutions to nonstationary Navier–Stokes equations. A similar approach was often used by O. Ladyzhenskaya [23] for various nonlinear problems in the theory of viscous fluids and other partial differential equations [24].

It should be noted that roots of the Galerkin method and weak solutions can be found in much earlier studies related to the Virtual Work Principle in mechanics that dates back to J. D’Alembert (who used it in order to derive equations of motion for a mechanical system of rigid bodies) and J.-L. Lagrange (who suggested a more general principle of virtual displacements for continuum media problems). Of course the development of mathematics at that time was insufficient to correctly determine what should be considered as “the set of virtual displacements”.

3 Further Development of the Method

Modern numerical methods are based on the development of principle ideas encompassed in the Galerkin method. The amount of related publications is vast, and in this short note it is impossible to present a systematic overview of them. Therefore, we will briefly mention only a few areas that seem to be the most important.

In 1943, R. Courant [12] suggested a version of the method based on locally supported test functions. This was an important step, which generated geometrically flexible numerical schemes with dispersed resolving matrixes. Soon this version of the method gained wide popularity and was named the Finite Element Method (FEM). Mathematical analysis of FEM is based upon two fundamental relations: Galerkin orthogonality and projection estimate or Cea lemma (see, e.g., [7, 10]), which states that the distance between u and the Galerkin solution uN is controlled by the distance between u and the respective finite dimensional space XN.

The projection inequality forms the basis of a priori analysis of the accuracy of Galerkin approximations. The development of this direction was closely connected with the solution of purely mathematical problems related to functional analysis, theory of interpolation, and regularity of weak solutions. It was established that ultimate accuracy of finite-dimensional approximations is determined by the quantity called Kolmogorov’s N width of a compact K⊂X (see [21]):

This quantity is extremely important in modern approximation theory, which investigates (6) for various functional and approximation spaces.

The Galerkin method has stimulated other important mathematical studies, especially in numerical linear algebra. It led to the creation of domain decomposition and multigrid methods. In the last decade, methods based on tensor type representations have been of exceptional interest.

In modern computational technologies, we everywhere find traces of the basic Galerkin idea. In mixed and hybrid finite element methods (see, e.g., [8]), it is applied to minimax formulations tested by pairs of functions in the respective primal–dual functional spaces. The discontinuous Galerkin method (DG) [2, 3, 5, 13] uses discontinuous test functions. Many commonly used methods (spectral methods, the finite volume method [15, 25], the weak Galerkin method, isogeometric analysis [11]) can be viewed as advanced versions of the Galerkin method.

The last but not the least direction that is inextricably linked with the Galerkin method is adaptivity. Here, the question

How to reform a subspace XN1 and get a better subspace XN2?

is in the focus of consideration. It is closely related to the theory of mesh adaptation (see, e.g., [29]) and a posteriori error indicators, which are often based on special properties of Galerkin approximations (see, e.g., [4, 9, 27, 32]).

4 International Conference CMAM

The above discussed directions and other close topics in applied mathematics and numerical analysis were widely represented at the conference Computational Methods in Applied Mathematics (CMAM-7), which was held in Jyväskylä (Finland) between July 31 and August 6, 2016, and dedicated to the jubilee of the Galerkin method. The conference was organized by the Department of Mathematical Information Technology of the University of Jyväskylä, together with the St. Petersburg Department of the V. A. Steklov Institute of Mathematics of the Russian Academy of Sciences (co-chairs Prof. P. Neittaanmäki and Prof. S. Repin). The conference was supported by the editorial board of the journal Computational Methods in Applied Mathematics and the Radon Institute of Computational and Applied Mathematics of the Austrian Academy of Sciences.

Traditionally, the CMAM conferences focus on various aspects of mathematical modelling and numerical methods for the approximate solution of problems arising in science and engineering and are supported by our journal having the same name. Previous CMAM conferences took place in Minsk (2003, 2007), Trakai (2005), Bedlewo (2010), Berlin (2012), and Strobl (2014).

Figure 3

Participants of CMAM-7.

Invited keynote lecturers of CMAM-7 were delivered by outstanding scientists in applied mathematics and numerical analysis, namely:

1. Carsten Carstensen (Humboldt University, Berlin, Germany);

2. Roland Glowinski (University of Houston, USA);

3. Boris Khoromskij (Max Planck Institute for Mathematics in the Sciences, Leipzig, Germany);

4. Yuri Kuznetsov (Institute of Numerical Mathematics, Moscow, Russia; University of Houston, USA);

5. Raytcho Lazarov (Texas A&M University, USA);

6. Markus Melenk (Technische Universität Wien, Austria);

7. Roderick Melnik (Wilfrid Laurier University, Ontario, Canada);

8. Stefan Sauter (University of Zurich, Switzerland);

9. Rolf Stenberg (Aalto University, Finland);

10. Ragnar Winther (University of Oslo, Norway);

11. Jun Zou (The Chinese University of Hong Kong, Hong Kong).

The conference also included seven minisimposia and many excellent contributed talks. All the information on the conference can be found on www.mit.jyu.fi/scoma/cmam2016/.

The majority of papers of this volume present carefully selected contributions of CMAM participants. The papers by Aghili–Di Pietro–Ruffini [1], Bærland–Lee–Mardal–Winther [6], Guo–Pan–He–Glowinski [16], Gustafsson–Stenberg–Videman [17], and Khoromskij–Repin [20] are related to the development of finite element methods and advanced computational technologies based on ideas of the Galerkin method. In the papers by Gustafsson–Stenberg–Videman [17], Khoromskij–Repin [20], Kruse–Wu [22], and Weymuth–Sauter–Repin [33], the reader will find new results related to error analysis. The paper by Picard [28] is devoted to mathematical foundations of numerical methods for piesoelectric evolutionary models. Tensor type methods for the Hartree–Fock equation are studied in the paper by Khoromskaia–Khoromskij [19]. Kruse–Wu [22] present new approaches to analysis of problems with randomized coefficients.