Galerkin method in abstract settings

 
This section deals with analyzing errors in finite element method. To determine when Galerkin method will produce a good approximation, abstraction is introduced.
Let $B:V \times V \rightarrow R$ be bounded bilinear form and let $F:V\rightarrow R$ be bounded linear form.
It is assumed that the problem to be solved can be stated as, find $u \in V$ such that
$$B(u,v)=F(v),v\in V$$
  Example
To make sense of above abstraction, it is best to see how this is derived using a weak form of PDE as shown below.
Consider Poisson's equation given by:
$$-\mathrm{\Delta }u=f\text{in }\mathrm{\Omega },u=0\text{on }\partial \mathrm{\Omega }$$
where $\mathrm{\Omega }$ is a bounded domain, $f$ is a given function, and $u$ is the function to be determined.


The weak formulation involves multiplying the differential equation by a test function $v$ from the space $H_0^1(\mathrm{\Omega })$, integrating over the domain $\mathrm{\Omega }$, and applying integration by parts:
$${\int }_{\mathrm{\Omega }}\mathrm{\nabla }u\cdot \mathrm{\nabla }vdx={\int }_{\mathrm{\Omega }}fvdx$$
Here, $B(u,v)={\int }_{\mathrm{\Omega }}\mathrm{\nabla }u\cdot \mathrm{\nabla }vdx$ and $F(v)={\int }_{\mathrm{\Omega }}fvdx$ define the bounded bilinear and linear forms respectively.

Well posed problem
The abstract setting problem is considered well-posed, if for each $F \in V^*$ , there exists a unique solution $u \in V$ and the mapping $F->u$ is bounded. $V^*$ here means dual space of $V$. Alternately, if $L:V \rightarrow V^*$ given by $<Lu,v>=B(u,v)$ is an isomorphism.

Example
For Dirichelet problem of Poisson's equation, we have
$$\begin{array}{r}V=\circ H^1(\mathrm{\Sigma })\\ B(u,v)={\int }_{\mathrm{\Sigma }}(gradu(x))(gradv(x))dx\\ F(v)={\int }_{\mathrm{\Sigma }}f(x)v(x)dx\end{array}$$
 A generalized Garlekin method for the abstract problem begins with a finite dimensional normed vector space $V_h$, a bileaner form $B_h:V_h \times V_h \rightarrow \mathcal{R}$ and a linear form $F_h:V_h \rightarrow \mathcal{R}$ and defines $u_h \in V_h$ by
 $$B_h(u_h,v)=F_h(v),v\in V_h$$
 Above equation can be written in the form $L_h u_h = F_h$ where $L_h:V_h\rightarrow V_h^*$ given by $<L_hu,v>=B_h(u,v)$
 If the finite dimensional problem is non-singular, then the norm of discrete solution operator known as ``stability constant'' is defined as
 $$\parallel L_h^{-1}\parallel$$

 In this approximation of the original problem determined by $V,B,F$ by $V_h,B_h,F_h$ intention is that $V_h$ in some sense approximates $V$ and $B_h,F_h$ approximate $B,V$. This is idea behind ``Consistency''.
 The goal here is to approximate $u_h$ to $u$. This is known as ``Convergence''.

 To this end, assume there is a restriction operator $\pi_h:V \rightarrow V_h$ so that $\pi_hu$ is close to $u$.
 Using the equation $L_h u_h = F_h$, we can compute the ``consistency error'' as
 $$L_h\pi u-F_h$$
 Error we wish to control is
 $${\pi }_{h}u-u$$
 Easy to see the relation between error and consistency error.
 $${\pi }_{h}u-u=L_h^{-1}(L_h\pi u-F_h)$$
 Recall thet norm of $L_h$ is stability constant. If we take norms both sides to above equation, we see that the norm of error is bounded by product of stability constant and norm of consistency error.
 $$\parallel {\pi }_{h}u-u\parallel \le \parallel L_h^{-1}\parallel \parallel (L_h\pi u-F_h)\parallel$$
Expressing this in terms of bilinear forms the relation becomes,
$$\parallel L_h{\pi }_{h}u-V_h\parallel =su{p}_{0\ne v\in V_h}\frac{B_h({\pi }_{h}u,v)-F_h(v))}{\parallel v\parallel }$$
Above equation, especially RHS needs some explanation.
The consistency error measures how close $B_h(\pi_hu,v)$ to $F_h(v)$ for all test functions $v$ in the subspace $V_h$. Here $F_h$ represents the discrete analog of forcing function. The ratio including sup means worst or max consistency error to the norm of test function $v$ over all possible non-zero test functions in $V_h$.
Finite dimensional problem is non-singular iff
$${\gamma }_h=in{f}_{0\ne u\in V_h}su{p}_{0\ne v\in V_h}\frac{B_h(u,v)}{\parallel u\parallel \parallel v\parallel }$$
And the stability constant is given by $\gamma^{-1}_h$
Notes:

Above $\gamma_h$ measures smallest ration of bilinear form $B_h(u,v)$ to the product of norms of $u$ and $v$ overall non-zero functions $u,v$ in the subspace $V_h$. Above condition shows that the bilinear form $B_h(u,v)$ is bounded from below and has a positive lower bound. This is also called ``coercive'' over subspace $V_h$. Since $B_h(u,v)$ is a matrix, above condition shows that this matrix is non-singular. For numerical methods, this means the problem is well-posed and having a continuous solution.

PDE- Galerkin method and example

 
The Galerkin method is a numerical technique that transforms boundary value problems, particularly partial differential equations (PDEs), into a system of linear algebraic equations by projecting the problem onto a finite-dimensional subspace.
Steps of the Galerkin Method


 Problem Setup: 

Start with the weak form of the differential equation integrated against a test function.

Choice of Basis Functions: 

Select basis functions $\{{\varphi }_i{\}}_{i=1}^n$ that satisfy the boundary conditions.
 Approximation of the Solution: 

Assume $u(x)\approx u_n(x)=\sum _{i=1}^n{c}_i{\varphi }_i(x)$.
 Galerkin Projection: 

Ensure the residual is orthogonal to the span of the basis functions:
  $$\int \text{Residual}\cdot {\varphi }_{j}dx=0\text{for all }j.$$

Example: One-Dimensional Poisson Equation

Consider the problem
$$-u^{″}(x)=f(x),u(0)=u(1)=0.$$
on the interval $[0,1]$.
Step 1: Weak Form

Multiply by a test function $v(x)$ and integrate by parts to get:
$${\int }_0^1{u}^{\prime }(x)v^{\prime }(x)dx={\int }_0^{1}f(x)v(x)dx.$$
Step 2: Discretization

Choose linear basis functions and assume:
$$u(x)\approx u_n(x)=\sum _{i=1}^n{c}_i{\varphi }_i(x).$$
Step 3: Galerkin Projection

Insert the approximation into the weak form using basis functions as test functions:
$$\sum _{i=1}^n{c}_i{\int }_0^1{\varphi }_i^{\prime }(x){\varphi }_j^{\prime }(x)dx={\int }_0^{1}f(x){\varphi }_j(x)dx\text{for all }j.$$
Step 4: Matrix Formulation

Define the stiffness matrix $\mathbf{A}$ and load vector $\mathbf{b}$ as follows:
$$A_{ij}={\int }_0^1{\varphi }_i^{\prime }(x){\varphi }_j^{\prime }(x)dx,b_j={\int }_0^{1}f(x){\varphi }_j(x)dx.$$
This leads to the linear system:
$$\mathbf{Ac}=\mathbf{b},$$
where $\mathbf{c}$ is the vector of coefficients.
\subsection*{Example: Very simple and concrete}
Take a differential equation
$$\frac{d^{2}y}{d{x}^2}+x+y=0,0\le x\le 1$$
with boundary condition $y(0)=y(1)=0$

Example:

Step 1
Take a trail solution $y(x)=a_0+a_{1}x+a_2{x}^2$. Always take the number of constants one greater than the highest differential power. Here we have second order differential. So, we take $2+1=3$ constants $a_0,a_1,a_2$.
Apply boundary conditions $x=0,y=0$. This leads to $a_0=0$. Then apply boundary condition $x=1,y=0$ leads to $a_1=-a_2$. Then the trail function becomes $y(x)=a_2(x^2-x)$. This is one parameter solution - namely $a_2$.

Step 2
Compute Weighting function
$$W(x)=\frac{\partial y}{a_2}=x^2-x$$

Step 3
Compute domain residuals. Simply substitute the trail function into original differential equation.
$$\begin{array}{rl}{R}_d& =\frac{d^{2}y}{d{x}^2}+x+y\\ R_d(x)& =\frac{d{a}_2(x^2-x)}{d{x}^2}+x+a_2(x^2-x)\end{array}$$
 Step 4
Minimization of domain residual:
$${\int }_0^{1}W(x)R_d(x)dx=0$$
Compute minimization. You will get $a_2=-\frac{5}{98}$.
Then the solution is $y(x)=-\frac{5}{98}(x-x^2)$.
Other formulations
For the differential equation $-u{}^{″}=f$, we can write this as first order by setting
$$\sigma =-u^{\prime },{\sigma }^{\prime }=f$$
The pair $\sigma ,f$ can be characterized variationally as unique critical point of the functional
$$I(\sigma ,u)={\int }_{-1}^1(\frac{1}{2}{\sigma }^2-u{\sigma }^{\prime })dx+{\int }_{-1}^{1}fudx$$
over $H^1(-1,1)\times L^2(-1,1)$
Note if we take the original $J(x)$ the Energy functional and substitute, we get above.
$$\begin{array}{rl}J(u)& =\frac{1}{2}{\int }_{-1}^1|u^{\prime }(x){|}^{2}dx-{\int }_{-1}^{1}f(x)u(x)dx\\ & =\frac{1}{2}{\int }_{-1}^1{\sigma }^2-{\int }_{-1}^{1}f(x)u(x)dx\\ & =\frac{1}{2}{\int }_{-1}^1{\sigma }^2+f(x){\int }_{-1}^{1}u(x)dx-{\int }_{-1}^1{\int }_{-1}^{1}u(x)dx{f}^{\prime }(x)dx\\ & =\frac{1}{2}{\int }_{-1}^1{\sigma }^2+0\\ & =\frac{1}{2}{\int }_{-1}^1{\sigma }^2+{\int }_{-1}^{1}f(x)u(x)dx-{\int }_{-1}^{1}f(x)u(x)dx\\ & ={\int }_{-1}^1\frac{1}{2}{\sigma }^2-f{\sigma }^{\prime }+{\int }_{-1}^{1}f(x)u(x)dx\end{array}$$

Hilbert space properties for PDEs

 In the context of PDEs the following properties are important for the Hilbert spaces.

Notation: Let $V$ be Hilbert space and let $a(.,.):V\times V \rightarrow \mathcal{R}$ be a bilinear form.
Property 1:
For PDE solutions, we need this bilinear form to be bounded.
$$|a(u,v)|<M||u||||v||\text{ for some }M>0\in \mathcal{R}$$
Property 2:
This is called $V$-Ellipticity'. The concept of $V$-ellipticity is crucial in establishing the well-posedness (existence, uniqueness, and stability of solutions) of boundary value problems formulated in a variational framework. It guarantees the uniqueness and stability of solutions to the corresponding variational problems. In essence, if the bilinear form derived from a PDE is V-elliptic, then the solution to the variational problem (and hence to the PDE) depends continuously on the data (such as boundary conditions and external forces), ensuring that small changes in input lead to small changes in the output.
Mathematically, it provides a lower bound for the bilinear form.
$$a(u,v)\ge \alpha ||v|{|}^2\alpha \text{ is a constant. }v\in V$$

Simple example
Simplest example of Hilbert space is Euclidean space with euclidean metric. It is not too difficult to verify that the euclidean metric is bilinear, it is bounded in the property $1$ sense and it has $V$-Ellipticity (Property $2$).

Weak formulation of boundary value PDE and its meaning

Energy functional
An energy functional is a mapping from a function space (often a Sobolev space) to the real numbers, which assigns a "total energy" value to each function in the space. The energy assigned typically depends on the function and its derivatives, reflecting physical or geometrical properties like potential energy, kinetic energy, or strain energy in various physical contexts.

$$J(u)=\frac{1}{2}{\int }_{-1}^1|u^{\prime }(x){|}^{2}dx-{\int }_{-1}^{1}f(x)u(x)dx$$
over $\mathring{H^1}(-1,1)$ where this Sobolev space is for functions that go to $0$ at boundary points for weak formulation.

Find $u\in \mathring{H^1}(-1,1)$such that
$${\int }_{-1}^1{u}^{\prime }(x)v^{\prime }(x)={\int }_{-1}^{1}f(x)v(x),v(x)\in \mathring{H^1}(-1,1)$$

Weak Formulation of a Boundary Value Problem

We consider a boundary value problem where we seek to find a function $u$ that satisfies the differential equation
$$-u^{″}(x)=f(x)\text{on}(-1,1),$$
with boundary conditions
$$u(-1)=u(1)=0.$$
Multiplying by a Test Function
To derive the weak form, multiply the differential equation by a test function $v(x)$, which is smooth and vanishes at the boundaries $(-1,1)$, hence $v(-1)=v(1)=0$. This test function $v(x)$ belongs to the space ${\mathring{H}}^1(-1,1)$, a subspace of $H^1(-1,1)$. We obtain:
$$-u^{″}(x)v(x)=f(x)v(x).$$
Integration Over the Domain
Integrate both sides over the interval $(-1,1)$:
$$-{\int }_{-1}^1{u}^{″}(x)v(x)dx={\int }_{-1}^{1}f(x)v(x)dx.$$
Integration by Parts
Use integration by parts on the left-hand side:
$$\begin{array}{rl}-{\int }_{-1}^1{u}^{″}(x)v(x)dx& ={[-u^{\prime }(x)v(x)]}_{-1}^1+{\int }_{-1}^1{u}^{\prime }(x)v^{\prime }(x)dx\\ & =0+{\int }_{-1}^1{u}^{\prime }(x)v^{\prime }(x)dx,\end{array}$$
                                                                                                        where the boundary terms vanish because $v(-1)=v(1)=0$.
Weak Formulation
Thus, we have the weak formulation of the boundary value problem:
$${\int }_{-1}^1{u}^{\prime }(x)v^{\prime }((x)dx={\int }_{-1}^{1}f(x)v(x)dx.$$
This equation must hold for all test functions $v(x)\in {\mathring{H}}^1(-1,1)$.
Meaning of the Weak Formulation
In this form, the differential equation $-u^{″}=f$ is translated into an integral equation that does not require the function $u$ to be twice differentiable. Instead, $u$ needs only to have its first derivative in $L^2(-1,1)$. This allows the inclusion of functions with less smoothness, accommodating more general solutions such as those exhibiting weak derivatives.

Simplest two-point boundary value problem

  

 Previous blog discussed Sobolev space definition and example

Here we focus on Simplest two-point boundary value problem.

Simplest two-point boundary value problem
$$-u^{{}^{″}}(x)=f(x),-1<x<1,u(-1)=u(1)=0$$
Energy functional
An energy functional is a mapping from a function space (often a Sobolev space) to the real numbers, which assigns a "total energy" value to each function in the space. The energy assigned typically depends on the function and its derivatives, reflecting physical or geometrical properties like potential energy, kinetic energy, or strain energy in various physical contexts.

$$J(u)=\frac{1}{2}{\int }_{-1}^1|u^{\prime }(x){|}^{2}dx-{\int }_{-1}^{1}f(x)u(x)dx$$
over $\mathring{H^1}(-1,1)$ where this Sobolev space is for functions that go to $0$ at boundary points for weak formulation.

Find $u\in \mathring{H^1}(-1,1)$such that
$${\int }_{-1}^1{u}^{\prime }(x)v^{\prime }(x)={\int }_{-1}^{1}f(x)v(x),v(x)\in \mathring{H^1}(-1,1)$$

 Example

Imagine a fluid flowing between two long, parallel plates separated by a distance, typically normalized for simplicity. The flow is assumed to be steady (unchanging with time) and laminar (smooth with no turbulence).
The fluid is driven by a pressure gradient or external force such as gravity.
 
 Navier-Stokes Equations: For an incompressible, steady flow with no change in flow direction, the Navier-Stokes equations can simplify significantly. For a flow driven by a pressure gradient in the x-direction, or an external force such as gravity, and assuming no dependence on the y coordinate (2D flow), the equation simplifies to:
 $$-\mu \frac{d^{2}u}{d{x}^2}=\frac{dp}{dx}+f$$
 where $\mu )$ is viscosity of the fluid, $\frac{dp}{dx}$ pressure gradient, $f$ additional force such as gravity.

 Assuming constant pressure gradient and assuming $f$ is constant, we can rewrite above equation as
 $$-u^{{}^{″}}(x)=g(x)$$
 where $g(x)=\frac{1}{\mu }(\frac{dp}{dx}+f)$.
 Assume walls are at $x=-1,x=1$ and the flow is zero here, that is $u(-1)=u(0)=0$

 

Sobolev Space definition with example

 Sobolev spaces are essential in functional analysis and partial differential equations.

 The space $H^1(-1,1)$ is a specific Sobolev space defined on the interval $(-1,1)$ on the real line.

Definition of $H^1(-1,1)$

The Sobolev space $H^1(-1,1)$ includes functions $u$ that satisfy both the function itself and its first derivative being in $L^2(-1,1)$. Specifically, the space is defined as follows:

•    A function $u$ belongs to $L^2(-1,1)$ if:

  $${\int }_{-1}^1|u(x){|}^{2}dx<\mathrm{\infty }$$

•   The first weak derivative $u^{\prime }$ must also be in $L^2(-1,1)$, which means:

  $${\int }_{-1}^1|u^{\prime }(x){|}^{2}dx<\mathrm{\infty }$$


Weak derivatives are used to accommodate functions that might not be differentiable in the classical sense everywhere on the interval, including functions that are continuous and differentiable almost everywhere but may have points of derivative discontinuity.

Norm in $H^1(-1,1)$}

The norm in the Sobolev space $H^1(-1,1)$ is defined to incorporate both the function and its derivative, given by:
$$\Vert u{\Vert }_{H^1}={({\int }_{-1}^1|u(x){|}^2+|u^{\prime }(x){|}^{2}dx)}^{1/2}$$

Importance and Applications

Sobolev spaces such as $H^1(-1,1)$ are crucial in the study of partial differential equations (PDEs), as they frequently arise as the solution spaces for various PDEs. Understanding Sobolev spaces is foundational for advanced studies in fields such as mathematical physics, engineering, and applied mathematics.

Example Function for $H^1(-1,1)$

Working with a concrete example clarifies above definitions.

Consider the function $u(x)=x^2$ as an example for a function in the Sobolev space $H^1(-1,1)$. We will verify that both $u$ and its first derivative $u^{\prime }$ belong to $L^2(-1,1)$.

\Verification of $u\in L^2(-1,1)$
First, we check if the function $u(x)=x^2$ is square integrable over $(-1,1)$:

$${\int }_{-1}^1|u(x){|}^{2}dx={\int }_{-1}^1{x}^{4}dx$$

Calculate the integral:

$${\int }_{-1}^1{x}^{4}dx={\left[\frac{x^5}{5}\right]}_{-1}^1=\frac{1}{5}-(-\frac{1}{5})=\frac{2}{5}$$

This result confirms that $u(x)=x^2$ is in $L^2(-1,1)$ as the integral is finite ($<\mathrm{\infty }$).

Verification of $u^{\prime }\in L^2(-1,1)$
Next, compute the derivative of $u(x)$ and check its square integrability:

$$u^{\prime }(x)=2x$$

Check if $u^{\prime }(x)=2x$ is square integrable:

$${\int }_{-1}^1|u^{\prime }(x){|}^{2}dx={\int }_{-1}^1(2x{)}^{2}dx=4{\int }_{-1}^1{x}^{2}dx$$

Computing the integral:

$${\int }_{-1}^1{x}^{2}dx={\left[\frac{x^3}{3}\right]}_{-1}^1=\frac{1}{3}-(-\frac{1}{3})=\frac{2}{3}$$

Thus:

$$4{\int }_{-1}^1{x}^{2}dx=4\cdot \frac{2}{3}=\frac{8}{3}$$

Since this integral is also finite, $u^{\prime }(x)=2x$ is in $L^2(-1,1)$.

Thus, the function $u(x)=x^2$ and its derivative $u^{\prime }(x)=2x$ are both in $L^2(-1,1)$, confirming that $u(x)=x^2$ is a valid example of a function in the Sobolev space $H^1(-1,1)$.