Finite Element Method (FEM)
I've taken a little break from neural networks lately to work on some code that solves Poisson's equation in two dimensions using finite element method. This is something I've been wanting to do for a long time now, but I've spent far more time lusting than doing. FEM is cool because one can deal with irregular domains. I've messed around a little with finite difference method (FDM), but this only works on regular domains. If you wanted use FDM on a non rectangular (or circular in polar coordinates) shape, even just an 'L'-shape, you'd have to do some nasty boundary condition matching, or just give up. If you do FEM on a rectangular domain, the matrix you solve ends up looking exactly like the one you use for FDM. In this post, I'm just going to form the weak form of the Poisson equation with zero boundary Dirichlet conditions.
Instead of starting directly with the Poisson equation, I want to develop FEM from a more general form of a second order PDE. When we get to the "weak" form of the problem (we'll define that in a hot sec) we'll just introduce the Poisson equation as a special case of the more general form of the problem. Below is a general form of a second order PDE. Note that this is not the most general form of a second order PDE -- we could have first derivative terms and mixed partial derivatives, etc.
$$
-\nabla \cdot (p(x,y)\nabla u(x,y)) + q(x,y)u(x,y) = f(x,y)
$$
Here $p$ $q$ and $f$ are random scalar functions, and $\nabla$ is our old friend, the del operator. The idea is to solve for $u$ on some domain $\Omega$. In a little bit, I'll use the notation $\partial \Omega$ to indicate the boundary of the domain $\Omega$. Now the goal is to form the 'weak' form of the equation. We do that by multiplying each side by some function $v$ and integrating over the whole domain. Note that I won't use a double integral sign to indicate 2D integrals. The differentials $da$ and $ds$ will tell you whether I mean an area or line integral, respectively.
$$
\int_\Omega (-\nabla \cdot (p(x,y)\nabla u(x,y)))v\:da = \int_\Omega fv\:da \\
-\int_\Omega \nabla \cdot (p \nabla u)v \: da + \int_\Omega quv \: da = \int_\Omega fv \: da
$$
Note that I dropped the $(x,y)$ just to make everything look a little cleaner. What happens next is pretty slick.
$$
\int_\Omega p \nabla u \cdot \nabla v \: da + \int_\Omega quv\:da - \int_{\partial \Omega} p \nabla u v \cdot\:d\vec{s} = \int_\Omega fv\: da \\
\int_\Omega (p \nabla u \cdot \nabla v + quv)\:da - \int_{\partial \Omega} p \nabla u v \cdot d\vec{s} = \int_\Omega fv \: da
$$
In order to figure this out, let's imagine that $p \nabla u$ is actually just some vector field $ \vec{\beta}$. What if I were to take the divergence of the product $v\vec{\beta}$?
$$
\nabla \cdot (v\vec{\beta}) = \nabla v \cdot \vec{\beta} + v(\nabla \cdot \vec{\beta})
$$
Now let's say $\vec{F} = v\vec{\beta}$. I know from the divergence theorem that
$$
\int_\Omega \nabla \cdot \vec{F}\:da = \int_{\partial \Omega} \vec{F} \cdot d\vec{s} \\
\int_\Omega \nabla \cdot (v\vec{\beta})\:da = \int_{\partial \Omega} v\vec{\beta} \cdot d\vec{s}
$$
We're almost there:
$$
\int_{\partial \Omega} v\vec{\beta} \cdot d\vec{s} = \int_\Omega \nabla v \cdot \vec{\beta} + v(\nabla \cdot \vec{\beta}) \\
\int_\Omega (\nabla \cdot \vec{\beta})v\:da = \int_{\partial \Omega} v\vec{\beta} \cdot d\vec{s} - \int_\Omega \nabla v \cdot \vec{\beta}
$$
Substituting $p \nabla u$ in for $\vec{\beta}$:
$$
\int_\Omega (\nabla \cdot p \nabla u)v\:da = \int_{\partial \Omega} vp \nabla u \cdot d\vec{s} - \int_\Omega \nabla v \cdot p \nabla u \: da
$$
This is exactly what we had written above (with an extra minus sign, and a rearranged p). Now we can deal with the more specific case of the Poisson equation. $\nabla^2 u = f$. Note that $p=1$ and $q=0$, which simplifies things considerably.
$$
\int_\Omega (\nabla u \cdot \nabla v) \: da - \int_{\partial \Omega} \nabla u v \cdot d\vec{s} = \int_\Omega fv\:da
$$
In the problem I'm tackling, I'm considering zero Dirichlet boundary conditions. In other words
$$
u(x,y)\vert_{\partial \Omega} = 0
$$
This simplifies the problem further, because we lose the line integral part of the weak form. In the end, the weak form for the Poisson equation given zero Dirichlet boundary conditions looks as follows.
$$
\int_\Omega (\nabla u \cdot \nabla v) \: da = \int_\Omega fv\:da
$$
Cool! Next time, I'll talk about discretizing the domain, and forming basis functions from which to construct an approximate solution to this continuously defined problem.
Instead of starting directly with the Poisson equation, I want to develop FEM from a more general form of a second order PDE. When we get to the "weak" form of the problem (we'll define that in a hot sec) we'll just introduce the Poisson equation as a special case of the more general form of the problem. Below is a general form of a second order PDE. Note that this is not the most general form of a second order PDE -- we could have first derivative terms and mixed partial derivatives, etc.
$$
-\nabla \cdot (p(x,y)\nabla u(x,y)) + q(x,y)u(x,y) = f(x,y)
$$
Here $p$ $q$ and $f$ are random scalar functions, and $\nabla$ is our old friend, the del operator. The idea is to solve for $u$ on some domain $\Omega$. In a little bit, I'll use the notation $\partial \Omega$ to indicate the boundary of the domain $\Omega$. Now the goal is to form the 'weak' form of the equation. We do that by multiplying each side by some function $v$ and integrating over the whole domain. Note that I won't use a double integral sign to indicate 2D integrals. The differentials $da$ and $ds$ will tell you whether I mean an area or line integral, respectively.
$$
\int_\Omega (-\nabla \cdot (p(x,y)\nabla u(x,y)))v\:da = \int_\Omega fv\:da \\
-\int_\Omega \nabla \cdot (p \nabla u)v \: da + \int_\Omega quv \: da = \int_\Omega fv \: da
$$
Note that I dropped the $(x,y)$ just to make everything look a little cleaner. What happens next is pretty slick.
$$
\int_\Omega p \nabla u \cdot \nabla v \: da + \int_\Omega quv\:da - \int_{\partial \Omega} p \nabla u v \cdot\:d\vec{s} = \int_\Omega fv\: da \\
\int_\Omega (p \nabla u \cdot \nabla v + quv)\:da - \int_{\partial \Omega} p \nabla u v \cdot d\vec{s} = \int_\Omega fv \: da
$$
In order to figure this out, let's imagine that $p \nabla u$ is actually just some vector field $ \vec{\beta}$. What if I were to take the divergence of the product $v\vec{\beta}$?
$$
\nabla \cdot (v\vec{\beta}) = \nabla v \cdot \vec{\beta} + v(\nabla \cdot \vec{\beta})
$$
Now let's say $\vec{F} = v\vec{\beta}$. I know from the divergence theorem that
$$
\int_\Omega \nabla \cdot \vec{F}\:da = \int_{\partial \Omega} \vec{F} \cdot d\vec{s} \\
\int_\Omega \nabla \cdot (v\vec{\beta})\:da = \int_{\partial \Omega} v\vec{\beta} \cdot d\vec{s}
$$
We're almost there:
$$
\int_{\partial \Omega} v\vec{\beta} \cdot d\vec{s} = \int_\Omega \nabla v \cdot \vec{\beta} + v(\nabla \cdot \vec{\beta}) \\
\int_\Omega (\nabla \cdot \vec{\beta})v\:da = \int_{\partial \Omega} v\vec{\beta} \cdot d\vec{s} - \int_\Omega \nabla v \cdot \vec{\beta}
$$
Substituting $p \nabla u$ in for $\vec{\beta}$:
$$
\int_\Omega (\nabla \cdot p \nabla u)v\:da = \int_{\partial \Omega} vp \nabla u \cdot d\vec{s} - \int_\Omega \nabla v \cdot p \nabla u \: da
$$
This is exactly what we had written above (with an extra minus sign, and a rearranged p). Now we can deal with the more specific case of the Poisson equation. $\nabla^2 u = f$. Note that $p=1$ and $q=0$, which simplifies things considerably.
$$
\int_\Omega (\nabla u \cdot \nabla v) \: da - \int_{\partial \Omega} \nabla u v \cdot d\vec{s} = \int_\Omega fv\:da
$$
In the problem I'm tackling, I'm considering zero Dirichlet boundary conditions. In other words
$$
u(x,y)\vert_{\partial \Omega} = 0
$$
This simplifies the problem further, because we lose the line integral part of the weak form. In the end, the weak form for the Poisson equation given zero Dirichlet boundary conditions looks as follows.
$$
\int_\Omega (\nabla u \cdot \nabla v) \: da = \int_\Omega fv\:da
$$
Cool! Next time, I'll talk about discretizing the domain, and forming basis functions from which to construct an approximate solution to this continuously defined problem.