VariationalPrinciple_and_WeakForm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Variational Principle and Weak Form\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The main difference between FEM and FDM is that \n",
    "  1. FEM solves PDEs in their **weak form**\n",
    "  2. while the Taylor expansion was the central mathematical tool for deriving FDM.\n",
    "\n",
    "- As we'll see, the weak form of a PDE involves **integration** as a part of its definition.\n",
    "  - It is not a coincidence that we need fairly heavy mathematical machinery for integration and integrable functions in a general theory for FEM.\n",
    "  - For instance, modern theories for FEM assumes *Lebesgue integral*, which is more general than the Riemann integration we are familiar with.\n",
    "\n",
    "- For this reason, it's natural to assume that approximate solutions we acquire with FEM are integrable functions. \n",
    "\n",
    "- A theory of FEM often starts from defining spaces of integrable functions like $L_{2}$, **the space of square-integrable functions**,\n",
    "  - a scalar product of functions\n",
    "  $$\n",
    "    (u,v) = \\int_{\\Omega} u(x)\\,v(x)dx,\n",
    "  $$\n",
    "  - and the induced norm of a function like the famous $L_{2}$ norm: \n",
    "$$\n",
    "||u||_{L_{2}} = (u,u) = \\int_{\\Omega} u(x)\\,u(x)dx.\n",
    "$$\n",
    "\n",
    "- Since we are solving PDEs, **Sobolev space**, $H^{m}(\\Omega)$, plays a crucial role in the theory, which is a space of $L_{2}$ functions with square-integrable partial derivatives up to order $m$.\n",
    "\n",
    "\n",
    "- We are not going to go deep into these mathematical theories. Instead, we'll try to build intuitive understanding through examples.\n",
    "  \n",
    "- However, it is still worth knowing some of abstract formulation to understand the existing literature, keep up with the advances of the method and to enjoy unified treatment of many different problems.\n",
    "  \n",
    "- Those who are interested in mathematical foundation of FEM are referred to introductory books like the one by Johnson\n",
    "  - *Numerical Solution of Partial Differential Equations by the Finite Element Method*, Claes Johnson, Cambridge University Press, 1988, p 280.\n",
    "  - As of 2009, a Dover edition available! \n",
    "\n",
    "- Other than the weak form, a working FEM further needs the following components: \n",
    "  1. Choice of basis functions (e.g., orthogonal functions, hat functions, etc), \n",
    "  2. numerical integrations in space (quadrature) and time and \n",
    "  3. efficient solvers for a large sparse linear system.\n",
    "  \n",
    "- We have already dipped our toes into all of these numerical techniques. So, once we understand the weak formulation, it would be straightforward to solve some simple problems using FEM.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic FEM: Weak form of an elliptic PDE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation*}\n",
    "\\text{(D)} \\quad -u'' = f\n",
    "\\end{equation*}\n",
    "\n",
    "- We denote this formulation of the problem \"(D)\", meaning **differential**. Sometimes also called a **strong form** to distinguish it from the **weak form**.\n",
    "  \n",
    "- Now, we show that the solution $u$ of (D) is also the solution of **minimization problem** (M) and of **variational problem** (V).\n",
    "\n",
    "- Some notations and definitions:\n",
    "\n",
    "  - Scalar product of real-valued piecewise continuous bounded functions:\n",
    "\n",
    "  \\begin{equation*}\n",
    "    (v,w) = \\int_{0}^{1} v(x)\\, w(x)\\, dx.\n",
    "  \\end{equation*}\n",
    "\n",
    "  - Linear space $V$: \n",
    "  $$\n",
    "    V= \\{ v\\ |\\ v \\text{ is a continuous function on } [0,1],\\ v^{\\prime} \\text{ is piecewise continuous and bounded on } [0,1], \\text{ and } v(0)=v(1)=0 \\}.\n",
    "  $$\n",
    "\n",
    "  - A linear functional $F: V \\rightarrow R$ given by\n",
    "\n",
    "  \\begin{equation*}\n",
    "    F(v) = \\frac{1}{2} (v^{\\prime}, v^{\\prime})-(f,v).\n",
    "  \\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The problems (M) and (V) are defined as follows:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "& \\text{(M)} \\qquad \\text{Find } u \\in V \\text{ such that } F(u) \\le F(v),\\ \\forall v \\in V, \\\\\n",
    "& \\text{(V)} \\qquad \\text{Find } u \\in V \\text{ such that } (u^{\\prime}, v^{\\prime}) = (f,v),\\ \\forall v \\in V.\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "- In the solid mechanical context, $F(v)$ corresponds to the *total potential energy* associated with the displacement $v \\in V$. \n",
    "  - $\\frac{1}{2}(v^{\\prime}, v^{\\prime})$ represents the internal elastic energy and $(f,v)$ the work done by external force $f$.\n",
    "\n",
    "- For this reason, (M) corresponds to the fundamental **Principle of minimum potential energy** in mechanics. \n",
    "\n",
    "- (V), on the other hand, corresponds to the **Principle of virtual work** when $v$ is considered to be an arbitrary *virtual displacement*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Let's verify that the solution $u$ of (D) is a solution of (V).\n",
    "\n",
    "  - Multiply the equation $-u^{\\prime \\prime}=f$ by an *arbitrary* function $v \\in V$, a so-called **test function**. And integrate over the interval $(0,1)$ which yields\n",
    "\n",
    "  \\begin{equation*}\n",
    "    -(u^{\\prime\\prime},v)=(f,v).\n",
    "  \\end{equation*}\n",
    "  \n",
    "  - Integrate the l.h.s. by parts using the fact that $v(0)=v(1)=0$ to get\n",
    "\n",
    "  \\begin{equation*}\n",
    "    -(u^{\\prime\\prime},v)=-u^{\\prime}(1)v(1) + u^{\\prime}(0)v(0) + (u^{\\prime},v^{\\prime}) = (u^{\\prime},v^{\\prime}),\n",
    "  \\end{equation*}\n",
    "  \n",
    "  - $\\therefore$ $(u^{\\prime}, v^{\\prime}) = (f,v) \\quad \\forall v \\in V$, and $u$ is a solution of (V).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Showing (V) $\\Rightarrow$ (D) requires additional condition on $u$. \n",
    "\n",
    "  - Let's assume $u \\in V$ satisfies \n",
    "\n",
    "\\begin{equation*}\n",
    "  \\int_{0}^{1} u^{\\prime} v^{\\prime} dx - \\int_{0}^{1} f\\,v\\,dx = 0 \\quad \\forall v \\in V.\n",
    "\\end{equation*}\n",
    "\n",
    "  - If $u^{\\prime \\prime}$ exists and is **continuous**, then by integration by parts and the fact that $v(0)=v(1)=0$,\n",
    "\n",
    "\\begin{equation*}\n",
    "  -\\int_{0}^{1} \\left( u^{\\prime \\prime} +f \\right) v\\, dx =0 \\quad \\forall v \\in V.\n",
    "\\end{equation*}\n",
    "\n",
    "  - With the assumption that $(u^{\\prime\\prime}+f)$ is *continuous*, this relation can only hold if $(u^{\\prime\\prime}+f)=0 \\quad 0<x<1$. Thus $u$ is the solution of (D).\n",
    "\n",
    "- As just shown, a *regularity assumption* that $u^{\\prime \\prime}$ is continuous is required for the equivalence between (D) and (V). \n",
    "\n",
    "- However, finding (V) associated with (D) is (almost?) always possible.\n",
    "\n",
    "- The equivalence between (V) and (M) can be shown (cf Sec. 2.2 of Braess) without introducing extra restriction. However, the physical meaning of $F(u,v)$ is not always clear.\n",
    "\n",
    "- The formulation (V) is said to be a **weak formulation** of (D) and the solution of (V) is said to be a **weak solution** of (D).\n",
    "\n",
    "- Existence of a solution to (V) is easy to prove but that of (D) is relatively difficult. \n",
    "  - So, one often starts with the weak solution of (D) and shows that this solution is sufficiently regular to be also a solution to (D).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic FEM: Weak form of an elliptic PDE in 2 or 3D\n",
    "\n",
    "- Let's find the weak formulation for the following boundary value problem for the Poisson equation in a 2D domain:\n",
    "\n",
    "\\begin{align*}\n",
    "-\\Delta u &= f \\text{ in } \\Omega, \\\\\n",
    "        u &= 0 \\text{ on } \\partial \\Omega\n",
    "\\end{align*}\n",
    "where $\\Omega$ is a bounded open domain in the plane $R^{2}$ with boundary $\\partial \\Omega$, $f$ is a given continuous function and $\\Delta u = \\partial^{2} u/\\partial x_{1}^{2} + \\partial^{2} u/\\partial x_{2}^{2}$.\n",
    "\n",
    "- Multiply $v$, an arbitrary function that is continuous on $\\Omega$, has piecewise continuous $\\partial v/\\partial x_{1}$ and $\\partial v/\\partial x_{2}$ and 0 on $\\partial \\Omega$, to the above PDE and integrate over $\\Omega$ to get\n",
    "\n",
    "\\begin{equation}\n",
    "\\int_{\\Omega} -\\Delta u\\,v\\,d\\Omega = \\int_{\\Omega} f\\,v\\,d\\Omega.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Note that\n",
    "\n",
    "\\begin{equation}\n",
    "\\int_{\\Omega} \\nabla \\cdot (\\nabla u\\, v)\\, d\\Omega = \\int_{\\Omega} \\Delta u\\, v\\, d\\Omega + \\int_{\\Omega} \\nabla u \\cdot \\nabla v\\, d\\Omega.\n",
    "\\end{equation}\n",
    "\n",
    "- Applying **Green's formula**\n",
    "\n",
    "\\begin{equation}\n",
    "  \\int_{\\Omega} \\nabla \\cdot A \\, dx = \\int_{\\partial \\Omega} A \\cdot \\mathbf{n}\\, ds,\n",
    "\\end{equation}\n",
    "\n",
    "to the l.h.s. of the above equation, we get\n",
    "\n",
    "\\begin{equation}\n",
    "  \\int_{\\partial \\Omega} v\\,(\\nabla u \\cdot \\mathbf{n}) \\,dS = 0\n",
    "\\end{equation}\n",
    "\n",
    "since $v=0$ on $\\partial \\Omega$.\n",
    "\n",
    "- We get\n",
    "\n",
    "\\begin{equation}\n",
    "  -\\int_{\\Omega} \\Delta u\\, v\\, d\\Omega = \\int_{\\Omega} \\nabla u \\cdot \\nabla v\\, d\\Omega.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The BVP is given the following weak formulation:\n",
    "\n",
    "\\begin{equation}\n",
    "  (\\nabla u,\\ \\nabla v) = (f,\\,v) \\quad \\forall v \\in V,\n",
    "\\end{equation}\n",
    "where $V$ is the set of $v$'s with properties listed above.\n",
    "\n",
    "- We can try to find a solution a Sobolev space, which is bigger than $V$ and defined as follows: Let $\\Omega$ be a bounded domain in $R^{d}$, $d=2$ or $3$ and define\n",
    "\n",
    "\\begin{equation*}\n",
    "  L_{2}(\\Omega) = \\left\\{ v\\ |\\ v: \\Omega \\rightarrow R \\text{ and } \\int_{\\Omega} v^{2}\\,d\\Omega < \\infty \\right\\},\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "  H^{1}(\\Omega) = \\left\\{ v \\in L_{2}(\\Omega) \\ |\\ \\frac{\\partial v}{\\partial x_{i}} \\in L_{2}(\\Omega),\\ i=1, \\ldots, d \\right\\}.\n",
    "\\end{equation*}\n",
    "\n",
    "- Also introduce the corresponding scalar products and norms:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "& (v,\\ w)=\\int_{\\Omega} v\\,w\\,d\\Omega,\\ ||v||_{L_{2}(\\Omega)}=\\left( \\int_{\\Omega} v^{2}\\,d\\Omega \\right)^{1/2}, \\\\\n",
    "& (v,\\ w)_{H^{1}(\\Omega)}=\\int_{\\Omega} \\left[ v\\,w+\\nabla v \\nabla w \\right] \\,d\\Omega, \\\\\n",
    "& ||v||_{H^{1}(\\Omega)}=\\left( \\int_{\\Omega} \\left[ v^{2} + |\\nabla v|^{2} \\right] d\\Omega \\right)^{1/2}.\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "- Likewise, we can define \n",
    "\n",
    "\\begin{equation*}\n",
    "  H^{1}_{0}(\\Omega) = \\left\\{ v \\in H^{1}(\\Omega) \\ |\\ v=0\\ \\text{on}\\ \\partial \\Omega \\right\\}.\n",
    "\\end{equation*}  \n",
    "with the same scalar product and norm as $H^{1}(\\Omega)$.\n",
    "\n",
    "- Finally, the boundary value problem\n",
    "\n",
    "\\begin{equation*}\n",
    "\\text{(D)} \\qquad \\left.\\begin{aligned}\n",
    "&-\\Delta u = f \\quad \\text{ in } \\Omega,\\\\\n",
    "& u=0  \\quad \\text{ on } \\partial \\Omega.\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "\\end{equation*}\n",
    "\n",
    "can be given this **weak** or **variational** formulation:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\text{(V)} \\qquad \\text{Find } u \\in H^{1}_{0}(\\Omega) \\text{ such that } a(u,v)=(f,v) \\quad \\forall v \\in H^{1}_{0}(\\Omega),\n",
    "\\end{equation*}\n",
    "where \n",
    "\\begin{equation*}\n",
    "  a(u,v) = \\int_{\\Omega} \\nabla u \\cdot \\nabla v\\,dx, \\quad (f,\\ v)=\\int_{\\Omega} f\\,v\\,dx.\n",
    "\\end{equation*}\n",
    "\n",
    "- The B.C. given in this problem is called **Dirichlet** or **essential** condition. The same PDE can be given a different type of boundary condition\n",
    "\n",
    "\\begin{equation}\n",
    "  \\nabla\\,u \\cdot \\mathbf{n} = g \\quad \\text{ on } \\partial \\Omega,\n",
    "\\end{equation}\n",
    "which is called **Neumann** or **natural** condition\n",
    "\n",
    "- Incorporating the Neumann BC into the weak formulation is easy:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "& \\text{(V)} \\qquad \\text{Find } u \\in H^{1}(\\Omega) \\text{ such that } a(u,v)=(f,v) + <g,\\,v> \\\\\n",
    "& \\forall v \\in H^{1}(\\Omega),\n",
    "\\end{split}\n",
    "\\end{equation*}\n",
    "where \n",
    "\n",
    "\\begin{equation}\n",
    "  <g,\\,v>=\\int_{\\partial \\Omega} g\\,v\\, dS.\n",
    "\\end{equation}\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}