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.