{"id":308,"date":"2020-08-09T00:30:57","date_gmt":"2020-08-09T00:30:57","guid":{"rendered":"http:\/\/www.ethanepperly.com\/?p=308"},"modified":"2021-03-28T18:13:27","modified_gmt":"2021-03-28T18:13:27","slug":"big-ideas-in-applied-math-galerkin-approximation","status":"publish","type":"post","link":"https:\/\/www.ethanepperly.com\/index.php\/2020\/08\/09\/big-ideas-in-applied-math-galerkin-approximation\/","title":{"rendered":"Big Ideas in Applied Math: Galerkin Approximation"},"content":{"rendered":"\n

<\/p>\n\n\n\n

My first experience with the numerical solution of partial differential equations (PDEs)<\/a> was with finite difference methods<\/a>. I found finite difference methods to be somewhat fiddly: it is quite an exercise in patience to, for example, work out the appropriate fifth-order finite difference approximation to a second order differential operator on an irregularly spaced grid and even more of a pain to prove that the scheme is convergent. I found that I liked the finite element method<\/a> a lot better1<\/a><\/sup>Finite element methods certainly have their own fiddly-nesses (as anyone who has worked with a serious finite element code can no doubt attest to).<\/span> as there was a unifying underlying functional analytic theory, Galerkin approximation<\/a>, which showed how, in a sense, the finite element method computed the best possible approximate solution to the PDE among a family of potential solutions. However, I came to feel later that Galerkin approximation was, in a sense, the more fundamental concept, with the finite element method being one particular instantiation (with spectral methods<\/a>, boundary element methods<\/a>, and the conjugate gradient method<\/a> being others). In this post, I hope to give a general introduction to Galerkin approximation as computing the best possible approximate solution to a problem within a certain finite-dimensional space of possibilities.<\/em><\/p>\n\n\n\n

Systems of Linear Equations<\/h2>\n\n\n\n

Let us begin with a linear algebraic example, which is unburdened by some of the technicalities of partial differential equations. Suppose we want to solve a very<\/em> large system of linear equations \"Ax = b\", where the matrix \"A\" is symmetric and positive definite (SPD)<\/a>. Suppose that \"A\" is \"N\times N\" where \"N\" is so large that we don’t even want to store all \"N\" components of the solution \"x\" on our computer. What can we possibly do?<\/p>\n\n\n\n

One solution is to consider only solutions \"x\" lying in a subspace<\/a> \"\mathcal{X}\" of the set of all possible solutions \"\mathbb{R}^N\". If this subspace has a basis \"x_1,x_2,\ldots,x_M \in \mathcal{X}\", then the solution \"x \in \mathcal{X}\" can be represented as \"x = a_1x_1 + \cdots + a_Mx_M\" and one only has to store the \"M < N\" numbers \"a_1,\ldots,a_M\". In general, \"x\" will not belong to the subspace \"\mathcal{X}\" and we must settle for an approximate solution \"\hat{x} \in \mathcal{X}\".<\/p>\n\n\n\n

The next step is to convert the system of linear equations \"Ax = b\" into a form which is more amenable to approximate solution on a subspace \"\mathcal{X}\". Note that the equation \"Ax = b\" encodes \"n\" different linear equations \"a_i^\top x = b_i\" where \"a_i^\top\" is the \"i\"th row of \"A\" and \"b_i\" is the \"i\"th element of \"b\". Note that the \"i\"th equation is equivalent to the condition \"e_i^\top A x = e_i^\top b\", where \"e_i\" is the vector with zeros in all entries except for the \"i\"th entry which is a one. More generally, by multiplying the equation \"Ax = b\" by an arbitrary test row vector \"y^\top\", we get \"y^\top Ax = y^\top b\" for all \"y \in \mathbb{R}^N\". We refer to this as a variational formulation of the linear system of equations \"Ax = b\". In fact, one can easily show that the variational problem is equivalent to the system of linear equations:<\/p>\n\n\n

(1) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

Since we are seeking an approximate solution from the subspace \"\mathcal{X}\", it is only natural that we also restrict our test vectors \"y\" to lie in the subspace \"\mathcal{X}\". Thus, we seek an approximate solution \"\hat{x}\" to the system of equations \"Ax = b\" as the solution of the variational problem<\/p>\n\n\n

(2) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

One can relatively easily show this problem possesses a unique solution \"\hat{x}\".2<\/a><\/sup>Here is a linear algebraic proof. As we shall see below, the same conclusion will also follow from the general Lax-Milgram theorem<\/a>. Let \"P\" be a matrix whose columns form a basis for \"\mathcal{X}\". Then every \"y\in \mathcal{X}\" can be written as \"y = Pu\" for some \"u \in \mathbb{R}^M\". Thus, writing \"\hat{x} = Pw\", we have that \"y^\top A\hat{x} = u^\top P^\top A P w = u^\top P^\top b\" for every \"u\in \mathbb{R}^M\". But this is just a variational formulation of the equation \"(P^\top A P)w = P^\top b\". The matrix \"P^\top A P\" is SPD since \"v^\top(P^\top A P)v = (Pv)^\top A (Pv) > 0\" for \"v\ne 0\" since \"A\" is SPD. Thus \"(P^\top A P)w = P^\top b\" has a unique solution \"w = (P^\top A P)^{-1}P^\top b\". Thus \"\hat{x} = Pw = P(P^\top AP)^{-1}P^\top b\" is the unique solution to the variational problem Eq. (2).<\/span> In what sense is \"\hat{x}\" a good approximate solution for \"Ax = b\"? To answer this question, we need to introduce a special way of measuring the error to an approximate solution to \"Ax = b\". We define the \"A\"-inner product of a vector \"x\" and \"y\" to be \"\langle x, y \rangle_A := y^\top Ax\" and the associated \"A\"-norm \"\|x\|_A = \sqrt{\langle x, x\rangle_A} = \sqrt{x^\top A x}\".3<\/a><\/sup>Note that \"A\"-norm can be seen as a weighted Euclidean norm, where the components of the vector \"x\" in the direction of the eigenvectors of \"A\" are scaled by their corresponding eigenvector. Concretely, if \"x = a_1 q_1 + \cdots + a_N q_N\" where \"q_j\" is an eigenvector of \"A\" with eigenvalue \"\lambda_j\" (\"Aq_j = \lambda_j q_j\"), then we have \"\|x\|_A^2 = \lambda_1 a_1^2 + \cdots + \lambda_N a_N^2\".<\/span> All of the properties satisfied by the familiar Euclidean inner product and norm<\/a> carry over to the new \"A\"-inner product and norm (e.g., the Pythagorean theorem<\/a>). Indeed, for those familiar, one can show \"\langle \cdot, \cdot \rangle_A\" satisfies all the axioms for an inner product space<\/a>.<\/p>\n\n\n\n

We shall now show that the error \"x - \hat{x}\" between \"x\" and its Galerkin approximation \"\hat{x}\" is \"A\"-orthogonal to the space \"\mathcal{X}\" in the sense that \"\langle y, x - \hat{x}\rangle_A = 0\" for all \"y \in \mathcal{X}\". This follows from the straightforward calculation, for \"y \in \mathcal{X}\",<\/p>\n\n\n

(3) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n

where \"y^\top A x = y^\top b\" since \"x\" solves the variational problem Eq. (1) and \"y^\top A \hat{x} = y^\top b\" since \"x\" solves the variational problem Eq. (2).<\/p>\n\n\n

The fact that the error \"x - \hat{x}\" is \"A\"-orthogonal to \"\mathcal{X}\" can be used to show that \"\hat{x}\" is, in a sense, the best approximate solution to \"Ax = b\" in the subspace \"\mathcal{X}\". First note that, for any approximate solution \"z \in \mathcal{X}\" to \"Ax = b\", the vector \"\hat{x} - z \in \mathcal{X}\" is \"A\"-orthogonal to \"x - \hat{x}\". Thus, by the Pythagorean theorem,<\/p>\n\n\n

(4) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

Thus, the Galerkin approximation \"\hat{x}\" is the best approximate solution to \"Ax = b\" in the subspace \"\mathcal{X}\" with respect to the \"A\"-norm<\/em>, \"\|x - z\|_A \ge \|x - \hat{x} \|_A\" for every \"z \in \mathcal{X}\". Thus, if one picks a subspace \"\mathcal{X}\" for which the solution \"x\" almost lies in \"\mathcal{X},\"4<\/a><\/sup>In the sense that \"\inf_{z \in \mathcal{X}} \|x - z\|\" is small<\/span> then \"\hat{x}\" will be a good approximate solution to \"Ax = b\", irrespective of the size of the subspace \"\mathcal{X}\".<\/p>\n\n\n\n

Variational Formulations of Differential Equations<\/h2>\n\n\n\n

As I hope I’ve conveyed in the previous section, Galerkin approximation is not a technique that only works for finite element methods or even just PDEs. However, differential and integral equations are one of the most important applications of Galerkin approximation since the space of all possible solution to a differential or integral equation is infinite-dimensional: approximation in a finite-dimensional space is absolutely critical. In this section, I want to give a brief introduction to how one can develop variational formulations of differential equations amenable to Galerkin approximation. For simplicity of presentation, I shall focus on a one-dimensional problem which is described by an ordinary differential equation (ODE)<\/a> boundary value problem<\/a>. All of this generalized wholesale to partial differential equations in multiple dimensions, though there are some additional technical and notational difficulties (some of which I will address in footnotes). Variational formulation of differential equations is a topic with important technical subtleties which I will end up brushing past. Rigorous references are Chapters 5 and 6 from Evans’ Partial Differential Equations<\/em><\/a> or Chapters 0-2 from Brenner and Scott’s The Mathematical Theory of Finite Element Methods<\/em><\/a>.<\/p>\n\n\n\n

As our model problem for which we seek a variational formulation, we will focus on the one-dimensional Poisson equation<\/a>, which appears in the study of electrostatics, gravitation, diffusion, heat flow, and fluid mechanics. The unknown \"u\" is a real-valued function on an interval which take to be \"[0,1]\".5<\/a><\/sup>In higher dimensions, one can consider an arbitrary domain \"\Omega \subseteq \mathbb{R}^d\" with, for example, a Lipschitz boundary<\/a>. <\/span> We assume Dirichlet boundary conditions<\/a> that \"u\" is equal to zero on the boundary \"u(0) = u(1) = 0\".6<\/a><\/sup>In higher dimensions, one has \"u = 0\" on the boundary \"\partial \Omega\" of the region \"\Omega\".<\/span> Poisson’s equations then reads7<\/a><\/sup>\"-\Delta u = f\" on \"\Omega\" and \"u = 0\" for higher dimensions, where \"\Delta\" is the Laplacian operator<\/a>.<\/span>\n\n\n

(5) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

We wish to develop a variational formulation of this differential equation, similar to how we develop a variational formulation of the linear system of equations in the previous section. To develop our variational formulation, we take inspiration from physics. If \"u(x)\" represents, say, the temperature at a point \"x\", we are never able to measure \"u(x)\" exactly. Rather, we can measure the temperature in a region around \"x\" with a thermometer. No matter how carefully we engineer our thermometer, our thermometer tip will have some volume occupying a region \"R\" in space. The temperature \"u_{\rm meas}\" measured by our thermometer will be the average temperature in the region \"R\" or, more generally, a weighted average \"u_{\rm meas} = \int_0^1 u(x) v(x) \, dx\" where \"v(\cdot)\" is a weighting function which is zero outside the region \"R\". Now let’s use our thermometer to “measure” our differential equation:<\/p>\n\n\n

(6) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

This integral expression is some kind of variational formulation of our differential equation, as it is an equation involving the solution to our differential equation \"u\" which must hold for every<\/em> averaging function \"v\". (The precise meaning of every<\/em> will be forthcoming.) It will benefit us greatly to make this expression more “symmetric” with respect to \"u\" and \"v\". To do this, we shall integrate by parts<\/a>:8<\/a><\/sup>Integrating by parts is harder in higher dimensions. My personal advice for integrating by parts in higher dimensions<\/a> is to remember that integration by parts is ultimately a result of the product rule. As such, to integrate by parts, we first write an expression involving our integrand using the product rule of some differential operator and then integrate by both sides. In this case, notice that \"\nabla \cdot (v \nabla u) = v\Delta u + \nabla v \cdot \nabla u\". Rearranging and integrating, we see that \"\int_\Omega -v\Delta u \, dx = \int_\Omega \nabla v \cdot \nabla u \, dx - \int_\Omega \nabla \cdot (v\nabla u) \, dx\". We then apply the divergence theorem<\/a> to the last term to get \"\int_\Omega \nabla \cdot (v\nabla u) \, dx = \int_{\partial \Omega} v \nabla u \cdot n dS(x)\", where \"n\" represents an outward facing unit normal to the boundary \"\partial \Omega\" and \"dS(x)\" represents integration on the surface \"\partial \Omega\". If \"v\" is zero on \"\partial \Omega\", we conclude \"\int_\Omega \nabla v \cdot \nabla u \, dx = \int_\Omega v f \, dx\" for all nice<\/em> functions \"v\" on \"\Omega\" satisfying \"v = 0\" on \"\partial \Omega\".<\/span> <\/p>\n\n\n

(7) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

In particular, if \"v\" is zero on the boundary \"v(0) = v(1) = 0\", then the second two terms vanish and we’re left with the variational equation<\/p>\n\n\n

(8) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

Compare the variational formulation of the Poisson equation Eq. (8) to the variational formulation of the system of linear equations \"Ax = b\" in Eq. (1). The solution vector \"x\" in the differential equation context is a function \"u\" satisfying the boundary condition of \"u\" being zero on the boundary \"u(0) = u(1) = 0\". The right-hand side \"b\" is replaced by a function \"f\" on the interval \"[0,1]\". The test vector \"y\" is replaced by a test function \"v\" on the interval \"[0,1]\". The matrix product expression \"y^\top A x\" is replaced by the integral \"\int_0^1 v'(x)u'(x) \, dx\". The product \"y^\top b\" is replaced by the integral \"\int_0^1 v(x) f(x) \, dx\". As we shall soon see, there is a unifying theory which treats both of these contexts simultaneously.<\/p>\n\n\n\n

Before this unifying theory, we must address the question of which functions \"v\" we will consider in our variational formulation. One can show that all of the calculations we did in this section hold if \"v\" is a continuously differentiable function<\/a> on \"[0,1]\" which is zero away from the endpoints \"0\" and \"1\" and \"u\" is a twice continuously differentiable function on \"[0,1]\". Because of technical functional analytic considerations, we shall actually want to expand the class of functions in our variational formulation to even more functions \"v\". Specifically, we shall consider all functions \"v\" which are (A) square-integrable<\/a> (\"\int_0^1|v(x)|^2 \,dx\" is finite), (B) possess a square integrable derivative9<\/a><\/sup>More specifically, we only insist that \"v\" possess a square-integrable weak derivative<\/a>.<\/span> \"v'\" (\"\int_0^1|v'(x)|^2 \,dx\" is finite), and (C) are zero on the boundary. We refer to this class of functions as the Sobolev space<\/a> \"H_0^1((0,1))\".10<\/a><\/sup>The class of functions satisfying (A) and (B) but not necessarily (C) is the Sobolev space \"H^1((0,1))\". For an arbitrary function in \"H^1((0,1))\", the existence of a well-defined restriction to the boundary \"0\" and \"1\" is actually nontrivial to show, requiring showing the existence of a trace operator<\/a>. Chapter 5 of Evan’s Partial Differential Equations<\/em><\/a> is a good introduction to Sobolev spaces. The Sobolev spaces \"H_0^1((0,1))\" and \"H^1((0,1))\" naturally extend to spaces \"H_0^1(\Omega)\" and \"H^1(\Omega)\" for an arbitrary domain \"\Omega \subseteq \mathbb{R}^d\" with a nice boundary.<\/span> <\/p>\n\n\n\n

Now this is where things get really strange. Note that it is possible for a function \"u\" to satisfy the variational formulation Eq. (8) but for \"u\" not to satisfy the Poisson equation Eq. (5). A simple example is when \"f\" possesses a discontinuity (say, for example, a step discontinuity where \"f\" is \"0\" and then jumps to \"1\"). Then no continuously differentiable \"u\" will satisfy Eq. (5) at every point in \"\Omega\" and yet a solution \"u\" to the variational problem Eq. (8) exists! The variational formulation actually allows us to give a reasonable definition of “solving the differential equation” when a classical<\/em> solution to \"-u'' = f\" does not exist. Our only requirement for the variational problem is that \"u\", itself, belongs to the space \"H_0^1((0,1))\". A solution to the variational problem Eq. (8) is called a weak solution<\/a> to the differential equation Eq. (5) because, as we have argued, a weak solution to Eq. (8) need not always solve Eq. (5).11<\/a><\/sup>One can show that any classical solution to Eq. (5) solves Eq. (8). Given certain conditions on \"f\", one can go the other way, showing that weak solutions are indeed bonafide classical solutions. This is the subject of regularity theory<\/a>.<\/span>\n\n\n\n

The Lax-Milgram Theorem<\/h2>\n\n\n\n

Let us now build up an abstract language which allows us to use Galerkin approximation both for linear systems of equations and PDEs (as well as other contexts). If one compares the expressions \"y^\top A x\" from the linear systems context and \"\int_0^1 v'(x)u'(x) \, dx\" from the differential equation context, one recognizes that both these expressions are so-called bilinear forms<\/a>: they depend on two arguments (\"x\" and \"y\" or \"u\" and \"v\") and are a linear transformation<\/a> in each argument independently if the other one is fixed. For example, if one defines \"a(x,y) = y^\top A x\" one has \"a(x,\alpha_1 y_1 + \alpha_2 y_2) = \alpha_1 a(x,y_1) + \alpha_2 a(x,y_2)\". Similarly, if one defines \"a(u,v) = \int_0^1 v'(x)u'(x) \, dx\", then \"a(\alpha_1 u_1 + \alpha_2 u_2,v) = \alpha_1 a(u_1,v) + \alpha_2 a(u_2, v)\".<\/p>\n\n\n\n

Implicitly swimming in the background is some space of vectors or function which this bilinear form \"a(\cdot,\cdot)\" is defined upon. In the linear system of equations context, this space \"\mathbb{R}^N\" of \"N\"-dimensional vectors and in the differential context, this space is \"H_0^1((0,1))\" as defined in the previous section.12<\/a><\/sup>The connection between vectors and functions is even more natural if one considers a function \"u : [0,1] \to \mathbb{R}\" as a vector of infinite length, with one entry for each real number \"x \in [0,1]\".<\/span> Call this space \"\mathcal{V}\". We shall assume that \"\mathcal{V}\" is a special type of linear space called a Hilbert space<\/a>, an inner product space (with inner product \"\langle \cdot, \cdot \rangle_\mathcal{V}\") where every Cauchy sequence<\/a> converges to an element in \"\mathcal{V}\" (in the inner product-induced norm).13<\/a><\/sup>Note that every inner product space has a unique completion<\/a> to a Hilbert space. For example, if one considers the space \"C_{0}^\infty((0,1))\" of \"C^\infty\" smooth<\/a> functions which are zero away from the boundary of \"(0,1)\" with the inner product \"\langle u, v \rangle_{H^1_0((0,1))} = \int_0^1 u(x)v(x) + u'(x)v'(x) \, dx\", the completion is \"H_0^1((0,1))\". A natural extension to higher dimensions hold.<\/span> The Cauchy sequence convergence property, also known as metric completeness<\/a>, is important because we shall often deal with a sequence of entries \"u_1,u_2,\ldots \in \mathcal{V}\" which we will need to establish convergence to a vector \"u \in \mathcal{V}\". (Think of \"u_1,u_2,\ldots\" as a sequence of Galerkin approximations to a solution \"u\".)<\/p>\n\n\n\n

With these formalities, an abstract variational problem takes the form<\/p>\n\n\n

(9) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

where \"a(\cdot,\cdot)\" is a bilinear form on \"\mathcal{V}\" and \"\ell(\cdot)\" is a linear form<\/a> on \"\mathcal{V}\" (a linear map \"\ell: \mathcal{V} \to \mathbb{R}\"). There is a beautiful and general theorem called the Lax-Milgram theorem<\/a> which establishes existence and uniqueness of solutions to a problem like Eq. (9).<\/p>\n\n\n\n

Theorem (Lax-Milgram):<\/strong> Let \"a(\cdot,\cdot)\" and \"f\" satisfy the following properties:<\/p>\n\n\n\n

  1. (Boundedness of \"a\") There exists a constant \"C \ge 0\" such that every \"u,v \in \mathcal{V}\", \"|a(u,v)| \le C \|u\|_{\mathcal{V}}\|v\|_{\mathcal{V}}\".<\/li>
  2. (Coercivity) There exists a positive constant \"c > 0\" such that \"a(u,u) \ge c \|u\|_{\mathcal{V}}^2\" for every \"u \in \mathcal{V}\".<\/li>
  3. (Boundedness of \"\ell\") There exists a constant \"K\" such that \"|\ell(v)| \le K \|v\|_{\mathcal{V}}\" for every \"v\in \mathcal{V}\".<\/li><\/ol>\n\n\n\n

    Then the variational problem Eq. (9) possesses a unique solution<\/em>.<\/p>\n\n\n\n

    For our cases, \"a\" will also be symmetric \"a(u,v) = a(v,u)\" for all \"u,v \in \mathcal{V}\". While the Lax-Milgram theorem holds without symmetry, let us continue our discussion with this additional symmetry assumption. Note that, taken together, properties (1-2) say that the \"a\"-inner product<\/em>, defined as \"\langle u, v \rangle_a = a(u,v)\", is no more than so much bigger or smaller than the standard inner product \"\langle u, v\rangle_{\mathcal{V}}\" of \"u\" and \"v\".14<\/a><\/sup>That is, one has that the \"a\"-norm and the \"\mathcal{V}\"-norm are equivalent<\/a> in the sense that \"\sqrt{c}\|v\|_{\mathcal{V}} \le \|v\|_a \le \sqrt{C} \|v\|_{\mathcal{V}}\". so the norms \"\|\cdot\|_\mathcal{V}\" and \"\|\cdot\|_a\" define the same topology.<\/span>\n\n\n\n

    Let us now see how the Lax-Milgram theorem can apply to our two examples. For a reader who wants a more “big picture” perspective, they can comfortably skip to the next section. For those who want to see Lax-Milgram in action, see the discussion below.<\/p>\n\n\n

    <\/span>Applying the Lax-Milgram Theorem<\/div>
    Begin with the linear system of equations with \"\mathcal{V} = \mathbb{R}^N\" with inner product \"\langle x, y \rangle_{\mathbb{R}^N} = y^\top x\", \"a(x,y) = y^\top Ax\", and \"\ell(y) = y^\top b\". Note that we have the inequality \"\lambda_{\rm min} x^\top x \le x^\top A x \le \lambda_{\rm max} x^\top x\".15<\/a><\/sup>This is an easy consequence of the Courant-Fischer theorem<\/a>. More generally, note that, since \"A\" is symmetric, \"A\" has an orthonormal basis of eigenvectors \"q_1,\ldots,q_N\" with eigenvalues \"\lambda_{\rm max} = \lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_N = \lambda_{\rm min}\". Then \"x = x_1q_1 + \cdots + x_Nq_N\" and \"x^\top Ax = \lambda_{\rm max} x_1^2 + \cdots + \lambda_{\rm min} x_N^2\". The inequalities follow from noting the Parseval relation \"\|x\|^2 = x^\top x = x_1^2 + \cdots + x_N^2\" and noting that \"x^\top Ax / x^\top x\" is a convex combination<\/a> of the eigenvalues of \"A\".<\/span> In particular, we have that \"\|x\|_A = \sqrt{x^\top A x} \le \sqrt{\lambda_{\rm max}} \|x\|\". Property (1) then follows from the Cauchy-Schwarz inequality<\/a> applied to the \"A\"-inner product: \"|a(x,y)| = |\langle x, y\rangle_A| \le \|x\|_A \|y\|_A \le \lambda_{\rm max} \|x\|\|y\|\". Property (2) is simply the established inequality \"a(x,x) = x^\top Ax \ge \lambda_{\min} x^\top x\". Property (3) also follows from the Cauchy-Schwarz inequality: \"|\ell(y)| = |y^\top b| \le \|b\|\|y\|\". Thus, by Lax-Milgram, the variational problem \"y^\top Ax = y^\top b\" for \"y \in \mathbb{R}^N\" has a unique solution \"x\". Note that the linear systems example shows why the coercivity property (2) is necessary. If \"A\" is positive semi-definite but not positive-definite, then there exists an eigenvector \"v\" of \"A\" with eigenvalue \"0\". Then \"v^\top Av = 0 \not\ge c v^\top v\" for any positive constant \"c\" and \"A\" is singular, so the variational formulation of \"Ax = b\" has no solution for some choices of the vector \"b\".<\/p>\n

    Applying the Lax-Milgram theorem to differential equations can require powerful inequalities. In this case, the \"\mathcal{V} = H_0^1((0,1))\"-inner product is given by \"\langle u, v \rangle_{H_0^1((0,1))} = \int_0^1 v(x)u(x) + v'(x)u'(x) \, dx\", \"a(u,v) = \int_0^1 v'(x)u'(x) \, dx\", and \"\ell(v) = \int_0^1 v(x)f(x) \, dx\". Condition (1) is follows from a application of the Cauchy-Schwarz inequality for integrals<\/a>:16<\/a><\/sup> In higher dimensions, we need even more Cauchy-Schwarz! First, we note that the absolute value of integral is less than the integral of absolute value \"|a(u,v)| \le \int_\Omega |\nabla v \cdot \nabla u| \, dx\". Second, we apply the Cauchy-Schwarz inequality for the vectors \"\nabla u\" and \"\nabla v\" to get \"|\nabla v \cdot \nabla u| \le |\nabla v| |\nabla u|\" where, e.g., \"\|\nabla v\|\" is the Euclidean norm of the vector \"\nabla v\". This gives \"|a(u,v)| \le \int_\Omega \|\nabla v\| \|\nabla u\| \, dx\". Next, we apply the Cauchy-Schwarz inequality for integrals<\/a>\u00a0to get \"|a(u,v)| \le \le \left(\int_\Omega \|\nabla v\|^2 \, dx\right)^{1/2} \left(\int_\Omega \|\nabla u\|^2 \, dx\right)^{1/2}\". Finally, we note that \"\int_\Omega \|\nabla v\|^2 \, dx \le \int_\Omega v^2 + |\nabla v|^2 \, dx = \langle v, v\rangle_{H_0^1(\Omega)} = \|v\|_{H_0^1(\Omega)}^2\" and thus obtain \"|a(u,v)| \le \|v\|_{H_0^1(\Omega)}\|u\|_{H_0^1(\Omega)}\". This is the desired inequality with constant one.<\/span>\n

    (10) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n

    Let’s go line-by-line. First, we note that the absolute value of integral is less than the integral of absolute value. Next, we apply the Cauchy-Schwarz inequality for integrals<\/a>. Finally, we note that \"\int_0^1|v'(x)|^2 \, dx \le \int_0^1|v(x)|^2 + |v'(x)|^2 \, dx = \langle v, v\rangle_{H_0^1((0,1))} = \|v\|_{H_0^1((0,1))}^2\". This establishes Property (1) with constant \"C = 1\". As we already see one third of the way into verifying the hypotheses of Lax-Milgram, establishing these inequalities can require several steps. Ultimately, however, strong knowledge of just a core few inequalities (e.g. Cauchy-Schwarz) may be all that’s needed.<\/p>\n

    Proving coercivity (Property (2)) actually requires a very special inequality, Poincar\u00e9’s inequality<\/a>.17<\/a><\/sup>In higher dimensions, Poincar\u00e9’s inequality takes the form \"\int_\Omega \|\nabla u\|^2 \, dx \ge k\int_\Omega u^2 \, dx\" for a constant \"k\" depending only on the domain \"\Omega\".<\/span> In it’s simplest incarnation, the inequality states that there exists a constant \"k\" such that, for all functions \"u \in H_0^1((0,1))\",18<\/a><\/sup>A simple proof of Poincar\u00e9’s inequality<\/a>\u00a0for continuously differentiable \"u\" with \"u(0) = u(1) = 0\" goes as follows. Note that, by the fundamental theorem of calculus, \"u(x) = u(0) + \int_0^x u'(y) \, dy = \int_0^x u'(y) \, dy\". Applying the Cauchy-Schwarz inequality for integrals gives \"|u(x)|\le \left(\int_0^x 1^2 \, dy\right)^{1/2}\left(\int_0^x |u'(y)|^2 \, dy\right)^{1/2} \le \sqrt{x} \left(\int_0^x |u'(y)|^2 \, dy\right)^{1/2} \le \left(\int_0^1 |u'(y)|^2 \, dy\right)^{1/2}\" since \"0\lex\le 1\". Thus \"|u(x)|^2 \le \int_0^1 |u'(y)|^2 \, dy\" for all \"x \in [0,1]\" integrating over \"x\" gives \"\int_0^1\int_0^1|u'(y)|^2 \, dy\, dx = \int_0^1|u'(x)|^2 \, dx \ge \int_0^1|u(x)|^2 \, dx\". This proves Poincar\u00e9’s inequality with the constant \"k = 1\".<\/span>\n

    (11) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n

    With this inequality in tow, property (2) follows after another lengthy string of inequalities:19<\/a><\/sup>The same estimate holds in higher dimensions, with the appropriate generalization of Poincar\u00e9’s inequality.<\/span>\n

    (12) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n

    For Property (3) to hold, the function \"f\" must be square-integrable<\/a>. With this hypothesis, Property (3) is much easier than Properties (1-2) and we leave it as an exercise for the interested reader (or to a footnote20<\/a><\/sup>The proof is similar in one dimension or higher dimensions, so we state it for arbitrary domain \"\Omega\" for brevity. By Cauchy-Schwarz, we have that \"\ell(v) = \int_\Omega vf \, dx \le \left( \int_{\Omega} f^2 \, dx \right){1/2}\left( \int_{\Omega} v^2 \, dx \right){1/2} \le \left( \int_{\Omega} f^2 \, dx \right) \|v\|_{H_0^1(\Omega)}\".<\/span> for the uninterested reader).<\/p>\n

    This may seem like a lot of work, but the result we have achieved is stunning. We have proven (modulo a lot of omitted details) that the Poisson equation \"-u'' = f\" has a unique weak solution as long as \"f\" is square-integrable!21<\/a><\/sup>And in the footnotes, we have upgraded this proof to existence of a unique weak solution to the Poisson equation \"-\Delta u = f\" on a domain \"\Omega\".<\/span> What is remarkable about this proof is that it uses the Lax-Milgram theorem and some inequalities alone: no specialized knowledge about the physics underlying the Poisson equation were necessary. Going through the details of Lax-Milgram has been a somewhat lengthy affair for an introductory post, but hopefully this discussion has illuminated the power of functional analytic tools (like Lax-Milgram) in studying differential equations. Now, with a healthy dose of theory in hand, let us return to Galerkin approximation. <\/div><\/div>\n\n\n

    General Galerkin Approximation<\/h2>\n\n\n\n

    With our general theory set up, Galerkin approximation for general variational problem is the same as it was for a system of linear equations. First, we pick an approximation space \"\mathcal{X}\" which is a subspace of \"\mathcal{V}\". We then have the Galerkin variational problem<\/p>\n\n\n

    (13) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    Provided \"a\" and \"\ell\" satisfy the conditions of the Lax-Milgram theorem, there is a unique solution \"\hat{u}\" to the problem Eq. (13). Moreover, the special property of Galerkin approximation holds: the error \"u-\hat{u}\" is \"a\"-orthogonal to the subspace \"\mathcal{X}\". Consequently, \"\hat{u}\" is te best approximate solution to the variational problem Eq. (9) in the \"a\"-norm<\/em>. To see the \"a\"-orthogonality, we have that, for any \"v \in \mathcal{X}\", <\/p>\n\n\n

    (14) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    where we use the variational equation Eq. (9) for \"a(u,v) = \ell(v)\" and Eq. (13) for \"a(\hat{u},v) = \ell(v)\". Note the similarities with Eq. (3). Thus, using the Pythagorean theorem for the \"a\"-norm, for any other approximation solution \"w \in \mathcal{X}\", we have22<\/a><\/sup>Compare Eq. (4).<\/span>\n\n\n

    (15) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    Put simply, \"\hat{u}\" is the best approximation to \"u\" in the \"a\"-norm.23<\/a><\/sup>Using the fact the norms \"\|\cdot\|_{\mathcal{V}}\" and \"\|\cdot\|_a\" are equivalent in virtue of Properties (1-2), one can also show that \"\hat{u}\" is within a constant factor \"C/c\" of the best approximation in the norm \"\|\cdot\|_{\mathcal{V}\". This is known as C\u00e9a’s Lemma<\/a>.<\/span>\n\n\n\n

    Galerkin approximation is powerful because it allows us to approximate an infinite-dimensional problem by a finite-dimensional one. If we let \"\phi_1,\ldots,\phi_M\" be a basis<\/a> for the space \"\mathcal{X}\", then the approximate solution \"\hat{u}\" can be represented as \"\hat{u} = x_1 \phi_1 + \cdots + x_M \phi_M\". Since \"\phi_1,\ldots,\phi_N\" form a basis of \"\mathcal{X}\", to check that the Galerkin variational problem Eq. (13) holds for all \"v \in \mathcal{X}\" it is sufficient to check that it holds for \"v = \phi_1, v=\phi_2,\ldots,v=\phi_M\".24<\/a><\/sup>For an arbitrary \"v\in \mathcal{X}\" can be written as \"v = y_1 \phi_1 + \cdots y_M\phi_M\", so \"a(\hat{u},v) = a\left(\hat{u}, \sum_{j=1}^M y_j \phi_j \right) = \sum_{j=1}^M y_j a(\hat{u},\phi_j) = \sum_{j=1}^M y_j \ell(\phi_j) = \ell\left(\sum_{j=1}^M y_j \phi_j\right) = \ell(v)\".<\/span> Thus, plugging in \"\hat{u} = \sum_{j=1}^M x_j \phi_j\" and \"v = \phi_i\" into Eq. (13), we get (using bilinearity of \"a\")<\/p>\n\n\n

    (16) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    If we define \"a_{ij} = a(\phi_j,\phi_i)\" and \"b_i = \ell(\phi_i)\", then this gives us a matrix equation \"Ax = b\" for the unknowns \"x_1,\ldots,x_M\" parametrizing \"\hat{u}\". Thus, we can compute our Galerkin approximation by solving a linear system of equations.<\/p>\n\n\n\n

    We’ve covered a lot of ground so let’s summarize. Galerkin approximation is a technique which allows us to approximately solve a large- or infinite-dimensional problem by searching for an approximate solution in a smaller finite-dimensional space \"\mathcal{X}\" of our choosing. This Galerkin approximation is the best approximate solution to our original problem in the \"a\"-norm. By choosing a basis \"\phi_1,\ldots,\phi_M\" for our approximation space \"\mathcal{X}\", we reduce the problem of computing a Galerkin approximation to a linear system of equations.<\/p>\n\n\n\n

    Design of a Galerkin approximation scheme for a variational problem thus boils down to choosing the approximation space \"\mathcal{X}\" and a basis \"\phi_1,\ldots,\phi_M\". Picking \"\mathcal{X}\" to be a space of piecewise polynomial functions (splines) gives the finite element method<\/a>. Picking \"\mathcal{X}\" to be a space spanned by a collection of trigonometric functions gives a Fourier spectral method<\/a>. One can use a space spanned by wavelets<\/a> as well. The Galerkin framework is extremely general: give it a subspace \"\mathcal{X}\" and it will give you a linear system of equations to solve to give you the best approximate solution in \"\mathcal{X}\".<\/p>\n\n\n\n

    Two design considerations factor into the choice of space \"\mathcal{X}\" and basis \"\phi_1,\ldots,\phi_M\". First, one wants to pick a space \"\mathcal{X}\", where the solution \"u\" almost lies in. This is the rationale behind spectral methods. Smooth functions are very well-approximated by short truncated Fourier expansions<\/a>, so, if the solution \"u\" is smooth, spectral methods will converge very quickly. Finite element methods, which often use low-order piecewise polynomial functions, converge much more slowly to a smooth \"u\". The second design consideration one wants to consider is the ease of solving the system \"Ax = b\" resulting from the Galerkin approximation. If the basis function \"\phi_1,\ldots,\phi_M\" are local in the sense that most<\/em> pairs of basis functions \"\phi_i\" and \"\phi_j\" aren’t nonzero at the same point \"x\" (more formally, \"\phi_i\" and \"\phi_j\" have disjoint<\/a> supports<\/a> for most \"i\" and \"j\"), the system \"Ax = b\" will be sparse and thus usually much easier to solve<\/a>. Traditional spectral methods usually result in a harder-to-solve dense linear systems of equations.25<\/a><\/sup>There are clever ways<\/a> of making spectral methods which lead to sparse matrices. Conversely, if one uses high-order piecewise polynomials in a finite element approximation, one can get convergence properties similar to a spectral method. These are called spectral element methods<\/a>.<\/span> It should be noted that both spectral and finite element methods lead to ill-conditioned matrices<\/a> \"A\", making integral equation<\/a>-based approaches preferable if one needs high-accuracy.26<\/a><\/sup>For example, only one researcher using a finite-element method was able to meet Trefethen’s challenge to solve the Poisson equation to eight digits of accuracy on an L-shaped domain (see Section 6 of this paper<\/a>). Getting that solution required using a finite-element method of order 15!<\/span> Integral equations, themselves, are often solved using Galerkin approximation, leading to so-called boundary element methods<\/a>.<\/p>\n\n\n\n

    Upshot:<\/strong> Galerkin approximation is a powerful and extremely flexible methodology for approximately solving large- or infinite-dimensional problems by finding the best approximate solution in a smaller finite-dimensional subspace. To use a Galerkin approximation, one must convert their problem to a variational formulation and pick a basis for the approximation space. After doing this, computing the Galerkin approximation reduces down to solving a system of linear equations with dimension equal to the dimension of the approximation space.<\/p>\n","protected":false},"excerpt":{"rendered":"

    My first experience with the numerical solution of partial differential equations (PDEs) was with finite difference methods. I found finite difference methods to be somewhat fiddly: it is quite an exercise in patience to, for example, work out the appropriate fifth-order finite difference approximation to a second order differential operator on an irregularly spaced gridRead more<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[5],"tags":[],"class_list":["post-308","post","type-post","status-publish","format-standard","hentry","category-big-ideas-in-applied-math"],"_links":{"self":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/308","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/comments?post=308"}],"version-history":[{"count":115,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/308\/revisions"}],"predecessor-version":[{"id":669,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/308\/revisions\/669"}],"wp:attachment":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/media?parent=308"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/categories?post=308"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/tags?post=308"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}