[Math] How to apply Lagrange Multipliers to BCs of Time Dependent problems using finite elements

ap.analysis-of-pdesna.numerical-analysis

I am trying to implement a finite element scheme using the method of lines (finite difference in time and finite element in space) and enforcing boundary conditions using Lagrange Multipliers. This might be necessary for essential boundary conditions on non-Lagrangrian finite elements, such as Argyris finite elements.

For a Stationary problem I can formulate the FEM using the following linear system

$\begin{bmatrix} A & \Lambda^T \\ \Lambda & \varepsilon \end{bmatrix} \begin{bmatrix} u \\ \lambda\end{bmatrix} = \begin{bmatrix} \ell \\ \mathbf{0}\end{bmatrix}$

where $u$ is the variable we are trying to solve for, $\lambda$ is the Lagrange multipliers, $A$ might be the stiffness matrix, $\Lambda$ is the matrix of constraints, $\ell$ is the load vector, and $\varepsilon$ is a perturbation matrix to ensure the system is not singular. Thus, I can solve the system

$Au = \ell$
given the constraint

$\Lambda u = 0$,

where the constraint might be a Neumann, Dirichlet, or even Robin boundary condition. My question is how to extend this to a time dependent problem? As an easy example (read silly) how would one apply this to the problem

$\dfrac{\partial u}{\partial t} – \Delta u = f$ on $\Omega$

$u = 0$ on $\partial \Omega$

In the case that we set $\dfrac{\partial u}{\partial t} = 0$ the method amounts to $A = (\nabla u_i, \nabla v_j)$, $\ell = (f, v_j)$, and $\Lambda$ ` is a sparse matrix consisting of ones in the locations corresponding to the basis functions on the boundary, but once I have the time derivative I don't know exactly what to do, since there are no constraints on the time derivative. Usually, to solve the time dependent problem we might formulate the problem as

$\dfrac{d u}{d t} = F(t; u)$

and then apply something like a Runge-Kutta method. Unfortunately, I am unable to see how to come up with this formulation when using Lagrange Multipliers for the BCs. So, what I would like to ask is how do I apply the above scheme to a time dependent problem? References would be greatly appreciated.

Best Answer

In this answer I assume that by finite difference you mean only Forward Euler. Backward Euler and the RK are generalizations of this scheme that I will address at the end of the answer.

Like you said in your question, our goal is to come up with something like $$ \frac{\mathrm{d}}{\mathrm{d}t}u = F(t; u), $$

To account for the Langrange multipliers and the fact that problem has no time dependent boundary conditions, we modify this to $$ \frac{\mathrm{d}}{\mathrm{d}t}z = G(z) $$

where $z = (u,\lambda)$. With this goal in mind, we proceed.

Now, after discretizing the the PDE, we will get something like

$0=\dot{u}+ Au+\Lambda^\mathrm{T}\lambda$ and $0= \Lambda u+\epsilon \lambda$. Now, all we have to do to get these equations to the correct form is take the time derivative of the second equation and rearrange. Doing so, we get

$$ \left[ \begin{array}{c} \dot{u}\\ \dot{\lambda} \end{array} \right] = \left[ \begin{array}{cc} -A& -\Lambda^{\mathrm{T}}\\ \frac{1}{\epsilon}\Lambda A& \frac{1}{\epsilon} \Lambda\Lambda^\mathrm{T}\end{array}\right]\left[\begin{array}{c} u\\ \lambda\end{array}\right]. $$

And this what we wanted. This is the equation that we need to do Forward Euler, Backward Euler and most of the RK methods.

Related Solutions

[Math] How to apply Neuman boundary condition to Finite-Element-Method problems

You need to modify right-hand vector b of an equation Kx = b, where K is your stiffness matrix.

Here's how to do it, depending on which edge is on von Neumann boundary:

edge 1-2 (i.e. connecting local nodes 1 and 2),

l = sqrt(x21*x21 + y21*y21);
b[node1] += a * l / 2.0f;
b[node2] += a * l / 2.0f;
b[node3] += 0;

edge 1-3,

l = sqrt(x13*x13 + y13*y13);
b[node1] += a * l / 2.0f;
b[node2] += 0;
b[node3] += a * l / 2.0f;

edge 2-3,

l = sqrt(x32*x32 + y32*y32);
b[node1] += 0;
b[node2] += a * l / 2.0f;
b[node3] += a * l / 2.0f;

[Math] Tangential boundary conditions for magnetostatic FEM problem

$\newcommand{\v}[1]{\boldsymbol{#1}}$

First about the appropriate boundary condition we would like to impose: The cross product(scalar-$\mathrm{curl}$) in $\mathbb{R}^2$ is done by we embed a vector into $\mathbb{R}^3$: $\v{v} = \langle v_1,v_2,0\rangle$, then apply the cross product, what we get would be

$$\nabla \times \v{v} = \left\langle 0,0,\left|\begin{matrix}\partial_x&\partial_y\newline v_1&v_2 \end{matrix}\right|\right\rangle $$ Similarly the vector-$\mathbf{curl}$ is we embed a scalar function into the third dimension of $\mathbb{R}^3$: $u = \langle0,0,u\rangle$, and what we get is the rotation of the $\nabla$-operator normally denoted by $\nabla^{\perp}$: $$ \nabla^{\perp} u = \left\langle\frac{\partial u}{\partial y},-\frac{\partial u}{\partial x},0 \right\rangle $$ I read this way of defining 2D-curl operators in a famous computational fluid dynamics book Finite element methods for Navier-Stokes equations: theory and algorithms.

Then in $\mathbb{R}^2$ your problem can be written as: \begin{equation} \nabla^{\perp}(\mu^{-1}\nabla \times \v{A}) = \mu \v{J} \tag{1} \end{equation}

Suppose we multiply a test vector $\v{v}$ vanishing on boundary on both side or the equation and integration by parts:

$$ \int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v} - \int_{\partial \Omega} \mu^{-1}(\nabla \times \v{A})(\v{v} \times \v{n}) = \int_{\Omega} \mu\v{J}\cdot \v{v} $$ where recall in $\mathbb{R}^2$, $B = \mu^{-1}(\nabla \times \v{A})$ is really a scalar instead of vector, hence for the tangential Neumann boundary condition, doing the $\mathbb{R}^3$-embedding we have: $$ \v{n}\times B = \langle n_2 B, -n_1 B,0 \rangle $$ this is the value of the scalar $B$ along the tangential direction of the domain $\v{t} = \langle n_2, -n_1,0 \rangle$, dot product with $\v{t}$ itself we could get $B$($\v{t}$ is the rotation of a unit outward normal vector), hence the Neumann type boundary condition you really wanna impose here for 2D is: $$ \mu^{-1}(\nabla \times \v{A}) = g $$ where $g$ is a scalar function, and your variational form becomes: $$ \int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v} = \int_{\Omega} \mu\v{J}\cdot \v{v} + \int_{\partial \Omega} g(\v{v} \cdot \v{t}) \tag{2} $$ The finite element method is we triangulating the domain, construct a finite dimensional Hilbert space $\v{V}_h$ (usually a subspace of the solution space $\v{V}$ to above variational problem), the boundary integral is treated as:

$$ \int_{\partial \Omega} g(\v{v} \cdot \v{t}) = \sum_{e\subset \partial \Omega} \int_{e} g (\v{v}\cdot \v{t}_e) $$

where $e$ is a boundary edge, and here is the place quadrature kicks in, that you wanna evaluate a line integral on each edge, it depends on what your finite element test function space $\v{V}_h$, nowadays a consensus is to choose Nédélec elements, in lowest order case it is $\v{v}_e = \lambda_A\nabla \lambda_B - \lambda_B\nabla \lambda_A $ for $e = \overrightarrow{AB}$.

Secondly, what is really problematic with your statement here is: for a smooth vector field $\v{A} \in \v{C}^{\infty}(\mathbb{R}^2)$ or $\v{C}^{\infty}_c(\Omega)$, the equation (1) doesn't simply reduce to the scalar Laplace equation as you mentioned in OP, a direct observation would be that homogeneous Neumann problem for $-\Delta u = f$ has a constant as null space because a solution plus a constant is still a solution, while $\v{curl}$-problem has all the gradient vector fields as null space!

Similar thing applies for variational formulation (2), the bilinear form is not coercive, simply plugging a gradient field $\v{A} = \v{v} = \nabla \varphi$ into it. What we numerical analyst usually do here is to find a solution in the quotient space $\v{H}(\mathrm{curl})/\nabla H^1$ by adding a Lagrange multiplier to let the problem fall into mixed finite element(google MFEM if you would like to know more) framework:

$$ \begin{cases} \int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v} + \int_{\Omega} \v{v}\cdot \nabla \varphi= \int_{\Omega} \mu\v{J}\cdot \v{v} + \int_{\partial \Omega} g(\v{v} \cdot \v{t}) & \text{ for } \forall \v{v}\in \v{H}(\mathrm{curl})\newline \int_{\Omega} \v{A} \cdot \nabla \psi = 0& \text{ for } \forall \psi\in H^1/\mathbb{R} \end{cases} $$ to make sure problem (2) is solvable using finite element, not only we have to use Lagrange multiplier to circumvent the "null space" problem, we also have to guarantee that the right side is in the range of the curl-operator, here plug $\v{v} = \nabla \varphi$ into (2): $$ 0=\int_{\Omega} \mu\v{J}\cdot \nabla \varphi + \int_{\partial \Omega} g(\nabla \varphi \cdot \v{t}) = -\int_{\Omega} \nabla \cdot(\mu\v{J}) \varphi + \int_{\partial \Omega} (\mu\v{J} \cdot \v{n}) \varphi + \int_{\partial \Omega} g(\nabla \varphi \cdot \v{t}) $$ in your case, $g=0$, also by the arbitrariness of $\varphi$, for eg we choose $\nabla \varphi \cdot \v{t} = 0$ for $\varphi\in H^1_0$, we would get so called "compatibility condition" in finite element method for Neumann problem: $$ \nabla \cdot(\mu\v{J}) = 0\text{ in }\Omega\quad \text{ and }\quad \mu\v{J}\cdot \v{n} = 0\text{ on }\partial\Omega $$

Lastly about the implementation of finite element method for this equation, if you have a pure Neumann problem with homogeneous boundary condition(derivative on boundary is zero), then you do not have to do anything during assembling of the stiffness matrix or right hand side integral.

Best Answer

Related Solutions

[Math] How to apply Neuman boundary condition to Finite-Element-Method problems

[Math] Tangential boundary conditions for magnetostatic FEM problem

Related Question