Variational formulation of the problem

This chapter is the key one: here we lay the foundation of the finite element method. Having mastered it, you will be able to apply the same approach to equations of other types containing a linear operator.

We have considered the analytical approach to solving the heat conduction equation; however, for more complex geometries it quickly becomes cumbersome or altogether inapplicable, so now we proceed to our main topic — the numerical solution of the heat conduction equation. Let us write equation (1.8$\displaystyle \frac{\partial T(M, t)}{\partial t}$ $\displaystyle = a^2 \cdot \Delta T(M, t)$ $\displaystyle + f(M, t),$) for convenience, replacing $T$ by $u$

We need to find a function $u(M, t)$ satisfying the equation and the boundary conditions; we first homogenize the boundary conditions, as was done in the analytical part: the inhomogeneity goes into the heat source density function and the initial condition, and is added back to the solution at the end. The problem itself can be represented as the Euler equation

where $L$ $= \frac{\displaystyle \partial}{\displaystyle \partial t}$ $- a^2 \cdot \Delta$ is the linear operator of the differential equation, $u$ is the unknown function, $f$ is the load function. As I said above, our equation also fits this form, which means we can apply the well-known approach to solving the Euler equation. That approach is the Bubnov–Galerkin method: in the weak form it is required that for any trial function $\upsilon(M, t)$ at any moment of time $t$ the following equality holds

Why is all this needed, and what are these trial functions? Above, when constructing the analytical solution, we saw how the solution is expanded into a Fourier series — the same thing happens here: if the functions are chosen orthogonal, the problem splits into several solutions independent of each other, which finally add up to the sought one. The freedom of choice is very large, so the functions can be taken as simple and convenient for computation as possible. One glance at the equation above makes it clear that the solution can be represented by the sum $u(M, t)$ $= \sum_{i=1}^n a_i \cdot \upsilon_i(M, t)$, where $a_i$ are the coefficients to be found, and $\upsilon_i(M, t)$ are the basis functions; in the Bubnov–Galerkin method the trial functions are taken from the same family. The approach to the numerical solution is analogous to the analytical one and essentially also reduces to finding the coefficients of known basis functions. I would say more: there is no fundamental distinction between analytical and numerical methods — they are essentially the same thing.

Writing integrals is inconvenient, so let us write (5.2$\displaystyle \int_M L[u(M, t)] \cdot \upsilon(M, t) \,dM$ $\displaystyle = \int_M f(M, t) \cdot \upsilon(M, t) \,dM.$) in the form of a scalar product

It is known that the Euler equation is related to the functional as follows

Minimizing this functional brings us back to the Euler equation, and the function $\upsilon$ on which the minimum is attained is precisely the sought solution $u$. Let us show this.

If on the function $u$ the functional attains its minimum, let us give it the increment $u$ $+ \epsilon \cdot \upsilon$, where $\epsilon$ is a small number, and $\upsilon$ is an arbitrary function, and carry out the transformations taking into account the obvious relation $I(u)$ $\leq I(u + \epsilon \cdot \upsilon)$

When expanding the brackets we used the self-adjointness of the operator: $(L[\upsilon], u)$ $= (L[u], \upsilon)$. For the time derivative this is, strictly speaking, not the case — a careful treatment of this point is given in the appendix “The Bubnov–Galerkin functional”. Since $\epsilon$ can be negative, for the minimum condition to hold it is necessary that the term linear in $\epsilon$ vanish, that is, that (5.3$\displaystyle (L[u], \upsilon)$ $\displaystyle = (f, \upsilon).$) holds. Thus, to solve the Euler equation by the Bubnov–Galerkin method, one needs to minimize the functional (5.4$\displaystyle I(\upsilon)$ $\displaystyle = (L[\upsilon], \upsilon)$ $\displaystyle - 2 \cdot (f, \upsilon).$).

Solving equations analytically showed us the approach of separation of variables, so let us try to apply that knowledge here. Logically, nonstationarity adds no complexity related to the geometry: geometry separately, time separately. Ideally, solving the nonstationary equation should reduce to successively solving $N$ stationary problems. In the appendix “The Bubnov–Galerkin functional” we obtained the functionals whose minimization yields the solution of the nonstationary (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.$) and the stationary (I.5$\displaystyle a^2 \cdot \int_M \nabla \upsilon \cdot \nabla \upsilon \,dM$ $\displaystyle - 2 \cdot \int_M f \cdot \upsilon \,dM$ $\displaystyle - a^2 \cdot \int_S \upsilon \cdot \nabla \upsilon \,dS$ $\displaystyle \rightarrow \min.$) problems.