A Brief Introduction to the Weak Form

This is an introduction to the weak form for those of us who didn’t grow up using finite element analysis and vector calculus in our daily lives, but are nevertheless interested in learning about the weak form, with the help of some physical intuition and basic calculus.

About the Weak Form, Briefly

For many of the different types of physics simulated with COMSOL Multiphysics, a weak formulation, or weak form, is used behind the scenes to construct the mathematical model. Understanding the weak form will help us gain insight into how the COMSOL software works internally as well as enable us to write our own equations when there is no built-in interface available for the particular physics involved in our model.

You may also be interested in reading my colleague Bettina Schieche’s blog post “The Strength of the Weak Form“.

A Simple Example

Let’s consider a concrete example of 1D heat transfer at steady state with no heat source. Specifically, we are interested in the temperature, T, as a function of the position x in the domain defined by the interval $1\le x\le 5$. For simplicity, we assume the heat conductivity is unity. Then, the heat flux, q, in the positive x-direction is given by the gradient of the temperature T:

(1)

and the conservation of the heat flux (with no heat source in the domain) simply says

(2)

This is the main equation that we want to solve. Its solution will give us the temperature profile within the domain. Equations of this form appear in many different disciplines. For example, in electrostatics, T is replaced by the electric potential and q by the electric field, while in elasticity, T becomes the displacement and q becomes the stress.

Here, we start to see why COMSOL Multiphysics can solve coupled multiphysics problems in a breeze: No matter what physical mechanisms are involved, they are modeled by equations and once the equations are written down, they can be discretized and solved straightforwardly by the core algorithms of the COMSOL software.

Some readers may ask why we chose this seemingly too simple of an example, whose analytical solution can be easily obtained by very simple math or physical arguments. The reason is two-fold:

1. We want to focus on the central idea of the weak form and not be distracted by the math of a complicated physical system.
2. In subsequent posts, we will expand the example to more than one domain to demonstrate the coupling between two equation systems through their boundary conditions. Starting with a more complicated example now will almost certainly obscure the central theme later when the example is expanded.

The Weak Formulation

Equation (2) involves the first derivative of the heat flux, q, or the second derivative of the temperature, T, which may cause numerical issues in practical situations where the differentiability of the temperature profile may be limited. For example, at a boundary where the adjacent materials have different values of thermal conductivity, the first derivative of the temperature T becomes discontinuous and the second derivative of T can not be evaluated numerically. The main idea of the weak form is to turn the differential equation into an integral equation, so as to lessen the burden on the numerical algorithm in evaluating derivatives.

To turn the differential equation (2) into an integral equation, a naive first approach may be to integrate it over the entire domain $1\le x\le 5$:

This asks that the average value of $\partial_x q(x)$ in the entire domain is zero. Indeed, it seems “too weak” as compared to the original differential equation, which asks that $\partial_x q(x)$ should be zero everywhere in the interval $1\le x\le 5$. To improve upon it, we can ask instead that the average value of $\partial_x q(x)$ in a very narrow domain is zero, say,

The integral only involves the value of $\partial_x q(x)$ in the vicinity of x=3.5. Thus, the relation above requires it not to be too far away from zero: $\partial_x q(3.5) \approx 0$. Extending the same idea to all locations in the entire domain $1\le x\le 5$, we see that the original differential equation may be approximated by a set of integral equations, like this:

(3)

The higher the number of integral equations in the set, the better the approximation. In the limit of the infinite number of such integral equations, we recover the original differential equation. It is cumbersome or even impossible to write out all the integral equations in the set, but we can apply the same idea in a different way.

The main idea is to sample the value of $\partial_x q(x)$ in a narrow range. This is done by integrating it over a narrow range in Eq. (3) above. The same kind of sampling can be done by multiplying the integrand by a weight function \tilde{T}(x) that is non-trivial only in a narrow range, as shown pictorially below:

Then, we can integrate the product $\partial_x q(x)\tilde{T}(x)$ over the entire domain $1\le x\le 5$ for a variety of weight functions \tilde{T}(x). Each weight function limits the contribution of the integrand to a narrow range centered around a different x value, thus achieving the same effect as the collection of integral equations in Eq. (3). This leads us to the weak formulation, which states that the relation

(4)

should hold for a set of weight functions \tilde{T}(x), commonly called test functions. For every value of x, say x=3.5, we can choose a test function \tilde{T}(x) that is a narrow weight function centered around x=3.5. Plugging this test function into Eq. (4) would sample the value of $\partial_x q(x)$ in the vicinity of x=3.5 and so require it not to be too far away from zero: $\partial_x q(3.5) \approx 0$.

By plugging a large number of narrow weight functions as test functions into Eq. (4), each centered at a different location in the interval $1\le x\le 5$, the value of the function $\partial_x q(x)$ will be clamped down to zero everywhere within the domain.

Footnote: In the picture above, we intentionally plotted $\partial_x q(x)$ as an arbitrary curve, not the final solution to the equation, to emphasize the fact that we haven’t found the solution yet. Later on in the solution process, this arbitrary curve will be pushed up and down by a collection of test functions to reach the shape of the final solution.

Reducing Order of Differentiation

Note that the order of differentiation in the integrand of Eq. (4) is still the same as in Eq. (2) (after all, it’s the same function $\partial_x q(x)$), but it can be reduced using the method of integration by parts:

(5)

Now, there is no derivative of the heat flow, q, in the equation, or in terms of the temperature, T, the order of differentiation is reduced from two to one. What about the first derivative of the test function \tilde{T}(x), which just now appeared in the equation?

As we have seen in the previous section, the test function serves as a tool for us to find the solution to the equation. Thus, we have the freedom to choose any conveniently differentiable form for it.

Natural Boundary Condition

The first two terms of Eq. (5) involve the heat flux and test function at the domain boundaries x=1 and x=5, with the heat flux, q, defined in the positive x-direction. We can rewrite them in terms of the flux going out of the domain and move them to the right-hand side:

(6)

Here, $\Lambda$ is the outgoing flux, and the subscripts 1 and 2 represent the domain boundaries x=1 and x=5, respectively \Lambda_1\equiv -q(x=1) \mbox{, } \Lambda_2\equiv +q(x=5) \mbox{, } \tilde{T}_1\equiv \tilde{T}(x=1) \mbox{, } \tilde{T}_2\equiv \tilde{T}(x=5).

Also, we have used the heat flux relation (1) to write the integrand in terms of the temperature, T, and its test function, \tilde{T}. The right-hand side of the equation provides a natural way to assign boundary conditions in terms of the heat flux. The simplest is to set both \Lambda_1 and \Lambda_2 to zero to get insulating boundary conditions (no heat flux through the boundaries).

This is exactly the reason why in COMSOL Multiphysics the default boundary condition for heat transfer is “Thermal Insulation” and the one for solid mechanics is “Free (no boundary force)”. This kind of boundary condition, which specifies the flux or force (the first derivative of the variable being solved), is commonly called the natural boundary condition or the Neumann boundary condition.

Fixed Boundary Condition

Another type of boundary condition, commonly called the fixed boundary condition or the Dirichlet boundary condition, specifies the value of the variable being solved. In our current example, it specifies the value of the temperature at a point on the boundary. This kind of boundary condition is usually needed to set up a well-posed problem with a unique solution. For example, in fluid flow, we need to specify the pressure somewhere (not just the flow); and in solid mechanics, we need to specify the displacement somewhere (not just the force).

As we have seen in the example, the weak formulation provides a natural way to specify the heat flux at a boundary. How do we specify a fixed temperature at a boundary, then?

The trick is to take advantage of the mathematical structure of the natural boundary condition and apply the same idea of using test functions to clamp down the solution. Conceptually, to maintain a fixed temperature at a boundary point, a certain heat flux coming from the outside of the boundary is needed to compensate for the heat flux inside of the boundary. The weak formulation poses the problem as this: Find the heat flux needed to maintain the fixed temperature at the boundary point.

For example, if we want to specify the outgoing flux \Lambda to be 2 at x=1 and the temperature T to be 9 at x=5, then we introduce a new unknown variable \lambda_2 and its corresponding test function \tilde{\lambda}_2, and write Eq. (6) as:

(7)

Here, on the right-hand side, the first term specifies the outgoing flux of 2 at x=1 and the second term specifies the unknown flux at x=5; both terms come directly from the natural boundary condition terms on the right-hand side of Eq. (6).

The new variable \lambda_2 represents the unknown heat flux to be determined at the boundary x=5. The third term is added to the equation to force the solution to be T=9 at x=5 by means of the test function \tilde{\lambda}_2 in the same fashion as in the earlier discussion about the test function \tilde{T}(x).

Comment on Higher Dimensions

So far, we have been discussing a very simple one-dimensional example. In higher dimensions, such as a 2D surface domain or a 3D volume domain, the equations become more complicated, but the basic idea remains the same.

The weak formulation turns a differential equation into an integral equation. Integration by parts reduces the order of differentiation to provide numerical advantages, and generates natural boundary conditions for specifying fluxes at the boundaries. In the simple 1D example, the boundary is the two end points and the flux is a single value at each point.

In 2D and 3D, the boundary is the closed curve and the closed surface enclosing the domain, respectively. The right-hand side of Eq. (6) becomes the line or surface integral of the incoming flux density, in other words, the total incoming flux. In essence, the process of integration by parts in 2D and 3D uses the divergence theorem to obtain the line or surface integral of the flux at the boundary of the modeling domain.

In this blog post, we chose the simple 1D example so that the central idea wouldn’t be obscured by the complexity of the math.

Summary and Next Up

Today, we learned about the idea of the weak formulation of using test functions to clamp down the solution. Integrating the weak form by parts provides the numerical benefit of reduced differentiation order. It also provides a natural way to specify boundary conditions in terms of the fluxes or forces (the first derivatives of the variables being solved), the so-called natural boundary condition or the Neumann boundary condition.

When solving a practical problem, it’s often necessary to specify the variable being solved — not just its derivative — via the so-called fixed boundary condition or the Dirichlet boundary condition. We saw that the weak formulation uses the same mechanism of test functions and its natural boundary conditions to construct additional terms for the fixed boundary conditions.

So far, we have left the equations in their original analytical forms without any numerical approximation. In the next blog post, we will implement the weak form equation (7) in COMSOL Multiphysics to solve it numerically. After that, we will discuss how the numerical approximation is done internally, how the same problem can be solved in different ways, and how different boundary conditions can be set up for different types of problems.