I. The Bubnov–Galerkin functional

Let us write once more the functional for the Euler equation (5.4$\displaystyle I(\upsilon)$ $\displaystyle = (L[\upsilon], \upsilon)$ $\displaystyle - 2 \cdot (f, \upsilon).$)

The linear operator $L[\upsilon]$ for the parabolic equation we are interested in has the form

The scalar product $(L[\upsilon], \upsilon)$ has the form

The scalar product $(f, \upsilon)$ has the form

The Ostrogradsky formula for the Laplace operator is known

where $S$ is the boundary of the submanifold $M$, and $dS$ is its vector element directed along the outward normal. Let us integrate the second integral by parts

Thus, the original functional can be written as

Let us give the increment $\widehat{\upsilon}$ $+ \epsilon \cdot \upsilon$, where $\widehat{\upsilon}$ is the exact solution, and $\epsilon$ is some small number. The full increment of the functional can be written as

Neglecting the second-order terms proportional to $\epsilon^2$, we obtain

If we take $\epsilon$ out of the brackets and take into account that $\epsilon$ can be negative, while the increment of the functional is always positive, since $\widehat{\upsilon}$ is the exact solution, we obtain as a mandatory condition

Taking into account that $\widehat{\upsilon}$ $\approx \upsilon$, and cancelling the common factor, we return to the original equation. This means that the numerical solution of the boundary value problem reduces to the minimization of the functional

The factor $2$ in front of the linear term $\int_M f \cdot \upsilon dM$ is kept deliberately: it is precisely what ensures that the minimum condition $\delta I$ $= 0$ returns the original equation rather than an equation with an extra factor of $1/2$ (when differentiating, the quadratic terms produce a factor of $2$, and the linear term must have the same one).

Let us proceed to the separation into stationary and nonstationary problems. It is fundamentally important to perform this separation and the subsequent time discretization at the level of the equation, not of the functional: the time derivative is not a self-adjoint operator, so direct substitution of $\upsilon$ $= \psi(M) \cdot \phi(t)$ into the functional would lead to incorrect coefficients. Let us write the weak form of the equation $L[u]$ $= f$, obtained above by integration by parts

which must hold for an arbitrary trial function $\upsilon$.

The stationary equation with Dirichlet and/or Neumann boundary conditions

If $\frac{\displaystyle \partial u}{\displaystyle \partial t}$ $= 0$, the time term disappears, and due to the symmetry of the spatial operator the weak form (I.4$\displaystyle \int_M \frac{\displaystyle \partial u}{\displaystyle \partial t} \cdot \upsilon \,dM$ $\displaystyle + a^2 \cdot \int_M \nabla u \cdot \nabla \upsilon \,dM$ $\displaystyle - a^2 \cdot \int_S \upsilon \cdot \nabla u \,dS$ $\displaystyle = \int_M f \cdot \upsilon \,dM,$) is equivalent to the minimization of the functional

The nonstationary equation with Dirichlet and/or Neumann boundary conditions

We discretize the time derivative by the implicit Euler scheme with step $\Delta t$, applying it to the weak form of the equation: $\frac{\displaystyle \partial u}{\displaystyle \partial t}$ $\approx \frac{\displaystyle u_n - u_{n-1}}{\displaystyle \Delta t}$, where $u_n$ is the solution on the current time layer, and $u_{n-1}$ on the previous one. Substituting this into (I.4$\displaystyle \int_M \frac{\displaystyle \partial u}{\displaystyle \partial t} \cdot \upsilon \,dM$ $\displaystyle + a^2 \cdot \int_M \nabla u \cdot \nabla \upsilon \,dM$ $\displaystyle - a^2 \cdot \int_S \upsilon \cdot \nabla u \,dS$ $\displaystyle = \int_M f \cdot \upsilon \,dM,$) and multiplying by $\Delta t$, we obtain the equation for the time step

This equation, in turn, is the minimum condition of the functional

Thus, the nonstationary equation has been reduced to successively solving stationary problems: at each step, from the known $\upsilon_{n-1}$ is found $\upsilon_n$.

Matrix form

We discretize the trial function over the mesh nodes: $\upsilon$ $= \sum_i q_i \cdot \phi_i(M)$, where $\phi_i(M)$ are the basis functions, and $\vec{q}$ $= (q_0, \ldots, q_n)^T$ is the vector of nodal values. The substitution turns each integral of the functional into a quadratic or linear form in $\vec{q}$, and their coefficients are assembled into matrices:

where $K$ is the stiffness matrix (the integral of the product of the gradients), $D$ is the damping matrix (the integral of the product of the trial functions), $\vec{F}$ is the load vector (the integral of the source); the matrices $K$ and $D$ are symmetric. The boundary integral $-a^2 \int_S \upsilon \cdot \nabla \upsilon \,dS$ does not collapse into a matrix: it is a surface term that, in discrete form, enters only the equations of the boundary nodes (for an interior node the trial function vanishes on the boundary) and is determined by the boundary condition itself.

The discrete system is obtained by the Bubnov–Galerkin method: we substitute $\upsilon$ $= \sum_i q_i \phi_i$ into the weak form (I.4$\displaystyle \int_M \frac{\displaystyle \partial u}{\displaystyle \partial t} \cdot \upsilon \,dM$ $\displaystyle + a^2 \cdot \int_M \nabla u \cdot \nabla \upsilon \,dM$ $\displaystyle - a^2 \cdot \int_S \upsilon \cdot \nabla u \,dS$ $\displaystyle = \int_M f \cdot \upsilon \,dM,$) and successively take the trial functions $\phi_j$. For the stationary equation this gives the system

where $\vec{s}$ is the vector of nodal boundary fluxes. The flux itself is not specified here — what to do with it is determined by the type of boundary condition. Accounting for boundary conditions is a nontrivial task, and it is discussed in detail in separate appendices: “Accounting for Dirichlet boundary conditions” and “Accounting for Neumann boundary conditions”.

For the nonstationary equation (implicit Euler scheme, see (I.6$\displaystyle \int_M \upsilon_n^2 \,dM$ $\displaystyle + \Delta t \cdot a^2 \cdot \int_M \nabla \upsilon_n \cdot \nabla \upsilon_n \,dM$ $\displaystyle - 2 \cdot \int_M \left[ \Delta t \cdot f + \upsilon_{n-1} \right] \cdot \upsilon_n \,dM$ $\displaystyle - \Delta t \cdot a^2 \cdot \int_S \upsilon_n \cdot \nabla \upsilon_n \,dS$ $\displaystyle \rightarrow \min.$)) the system at the time step takes the form