Exploiting Maximum Principles to Save Time and Resources

When your simulations consume significant memory, do you buy a bigger computer? When they take too long to solve, do you just run them overnight? Often, you don’t have another option. But sometimes, if you have the right tools, you can find a better approach by exploiting the mathematical structure. Today, we will show you how to use the so-called maximum principles to save computational resources and time in the COMSOL Multiphysics® software.

Introduction to Maximum Principles

We often want to know the maximum or minimum values of some quantities. When the number of degrees of freedom is large, searching for extreme values can be computationally expensive. The cost is higher if we use the extreme values as inputs, such as in boundary conditions or material properties. One example is the design of controllers.

Imagine having prior information about possible critical locations. Now, we can focus the search to those locations. In time-dependent problems, we can save memory if we could store data only on such suspected areas.

The challenge is that for many problems, we do not have such information. But for some problems, we can prove that the extreme values are on the boundaries. In mathematical literature, such properties are called maximum principles.

COMSOL Multiphysics provides tools to look for extreme values just on boundaries. The software also lets us store solutions on a limited set of domains or boundaries.

Today, we will show you how to exploit maximum principles for efficient problem solving and postprocessing. First, we will discuss conditions under which maximum principles exist and sketch mathematical proofs. Next, we will demonstrate the tools that COMSOL Multiphysics provides to use these principles, including maximum and minimum operators on boundaries, physics interfaces defined on surfaces, and the boundary element method, among others. Finally, we will highlight common situations where the maximum principles do not hold.

Before we begin, note that while maximum principles are useful, we should only apply them in situations where the premises in them are valid. Otherwise, we will be led astray by the streetlight effect.

The streetlight effect occurs when you can’t find something because you are only searching where it is easy to look. It comes from a well-known joke in which a person looks for their wallet under a streetlight, even though they know they left it in a nearby park.

Deriving Maximum Principles

Let’s start with a simple one-dimensional problem. Consider the second-order ordinary differential equation

over an interval.

This equation can represent several 1D stationary boundary value problems such as heat conduction or chemical transport in the absence of advection or sources/sinks. Irrespective of the physics, the general solution is a straight line

where the constants a and b can be determined from boundary conditions.

As the solution is a straight line, its maximum and minimum value will be at the left or right end of the domain (interval). As such, if we are interested in the extreme values, we just need to check the values at the left and right ends. If both boundary conditions are of the Dirichlet type, prescribing u, we do not even need to solve the equation to find its maximum and minimum. The smaller of the Dirichlet boundary conditions is the minimum value and the larger is the maximum value. This obviously saves us some time.

Convection-Diffusion Problems with Uniform, Isotropic Properties

In the above analysis, we used our knowledge of the general form of the solution to make conclusions about the locations of extreme values. This does not generalize to a practical method for more complicated source terms or spatial dimensions higher than one. Instead of deriving a general form for the solution, let’s use properties of local maxima or minima from elementary calculus to make similar conclusions about the following linear partial differential equation.

where \Delta is the Laplace operator \Delta u = \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}+\frac{\partial^2u}{\partial z^2}; \nabla is the gradient operator \nabla u = (\frac{\partial u}{\partial x},\frac{\partial u}{\partial y},\frac{\partial u}{\partial z}); k is a positive scalar (diffusion or conduction coefficient); v is a convective velocity; and f is a source term. Notice that the source term is negative everywhere in the domain.

Can the solution have a maximum at an interior point of the domain? From optimality conditions for functions of several variables, we know that if a function is smooth, the gradient should be zero at critical points. Additionally, at a local maximum, all of the second partial derivatives \frac{\partial^2u}{\partial x^2}, \frac{\partial^2u}{\partial y^2}, \frac{\partial^2u}{\partial z^2} should be negative or zero. These two conditions result in a positive or zero left-hand side of the partial differential equation.

This contradicts with the source term

on the right-hand side of the same equation.

This means that the maximum of our solution can only be found on the boundary. What we just proved is called the strong maximum principle.

There is also the weak maximum principle, which holds when the source term is allowed to be zero; i.e., f\le0 is the condition. In this case, the maximum is either at the boundary or the solution is constant throughout. Either way, if we are looking for the maximum, it is enough to search just the boundary. We skip the proof of the weak maximum principle here. We can derive analogous principles for the minimum when f\geq0.

Finally, if there are no sources or sinks, we can combine the weak maximum and minimum principles to conclude that both the maximum and minimum of our solution should be on the boundary. It is customary to use the term maximum principles to refer to both maximum and minimum principles, with the distinction being understood from the context.

Notice that we require smoothness of the solution to make the above conclusions.

Anisotropic, Heterogeneous Convection-Diffusion Problems

Our conclusions so far can be used for stationary convection-diffusion-type physical problems such as heat transfer, current flow, or chemical transport, provided the conductivity (diffusivity) is homogenous and isotropic and the convection velocity v is homogenous as well. Let’s lift those restrictions now.

Typically, a stationary convection-diffusion problem has the form

where D is a conductivity (diffusion) tensor, which can be anisotropic and heterogeneous; v is a potentially nonuniform convective velocity; and f is a heat (chemical/current) source.

The requirements for establishing a weak maximum principle for a smooth solution are as follows:

1. The diffusion tensor D is positive definite
2. The convective velocity v is divergence free; i.e., \nabla \cdot \mathbf{v} = 0
3. The source term is nonpositive; i.e., f(\mathbf{x}) \leq 0, \quad \forall \mathbf{x} \in \Omega

Physically speaking, this covers transport problems with incompressible convective velocities. Heat transfer, reacting flows, and charge transport problems are some examples. In any case, if we have f \geq 0, we can conclude that the minimum should be on the boundary.

Consider the Joule heating problem (shown in detail in this video) in electric heaters, fuses, and other conductors where resistive losses in the electrical part lead to heating. The heat transfer in solids equation is given by

where k and σ are the thermal and electrical conductivity, and T and V are the temperature and electric potential.

The source term \sigma|\nabla V|^2 is always nonnegative. Based on our discussion so far, we can conclude that the minimum temperature has to be at the boundaries of our model. See the images below, where the locations of the volume and surface minimum have been annotated. The two locations overlap.

Temperature distribution in a Joule heating problem (left), location of the minimum overall temperature (center), and location of the minimum boundary temperature (right).

Transient Problems

Transient problems can broadly be divided into parabolic problems and hyperbolic problems. The first group contains diffusion-type phenomena such as heat transfer, chemical reactions, and groundwater flow. The second group contains wave-type problems such as acoustics, electromagnetic waves, and structural dynamics.

For parabolic problems, we can show that the overall maximum or minimum values of solutions are encountered either at the initial time or at the boundary. As in stationary problems, the sign of the source term plays a decisive role here. Hyperbolic problems in general do not obey maximum/minimum principles.

Consider the following parabolic equation

with appropriate initial and boundary conditions.

If the diffusion tensor D is positive definite, the convective velocity v is divergence free and the source term is such that f \leq 0. We can show that the maximum over all times and points is either at the boundaries or at the initial time. This is not a statement about the maximum at any given time. It is about the overall maximum.

To simplify the discussion here, let’s consider the strong maximum principle for a 1D version of the above parabolic equation.

Expand this equation with the product rule for derivatives to get

To prove the strong maximum principle, we will assume the maximum value of u to be at an interior point of the space-time domain, [0,L] \times [0,T], and show how that leads to a contradiction.

At an interior maximum, a smooth solution will have \frac{\partial u}{\partial x}=\frac{\partial u}{\partial t}=0. Additionally, an incompressible advection velocity will have \frac{\partial v}{\partial x}=0 everywhere for a 1D problem. Thus, at an interior maximum, the above equation reduces to

Since D is positive and f is negative, this means that at an interior maximum, \frac{\partial^2 u}{\partial x^2} >0. This contradicts the usual requirement from calculus that at an interior maximum, \frac{\partial^2 u}{\partial x^2} \leq 0.

This leaves us with the boundaries of the space-time domain as possible locations for the maximum solution. The boundary of this composite domain includes the initial time t = 0, the end time t = T, and the spatial boundaries at x = 0 and x = L.

A space-time domain for a 1D time-dependent PDE.

We narrowed down the location of the maximum solution to the four boundaries shown above. Our original statement was that the maximum is either at the spatial boundaries or at the initial condition. This excludes points that are spatially interior but at the terminal time t = T; i.e., the top of the rectangle in the above space-time domain. Let’s show why that is the case.

Say the maximum is at a spatially interior point \hat{x} \in (0,L) at the temporal boundary t = T. If we go left or right from that point, i.e., move in space but not in time, the solution should not increase. Therefore, \frac{\partial u}{\partial x}=0, \frac{\partial^2 u}{\partial x^2} \leq 0 . However, the first partial time derivative does not have to be zero, as we can only move in one direction in time from the terminal time. The restriction here is that for a maximum to happen at the terminal time, the time partial derivative cannot be negative there. That is, we need \frac{\partial u}{\partial t} \ge 0.. Does this agree with the partial differential equation? Setting the first spatial derivatives to zero, we have

From this, we can see that the starting assumption about the source term and the restriction on the second spatial partial derivative render a negative time derivative. This conclusion, based on the PDE, contradicts the statement we made based on requirements of a local maximum. Using a similar reasoning at the initial time, in contrast, shows that the time derivative at a possible interior maximum at the initial condition is a possibility.

In a computational work, this means that if we are looking for the maximum solution, we should limit our search to the boundaries, except at the initial time, where we should include interior points as well. Beware, however, that this statement only applies to the maximum of all times and points. It is still possible for the maximum at a given time to be at an internal point.

Systems and Higher-Order Equations

Maximum principles are not generally available for equations with an order higher than two, such as the biharmonic equation or second-order systems such as two- and three-dimensional stress analysis. The basic idea in looking for maximum principles is to get some scalar function P(u) of the solution vector u to satisfy the Poisson equation. Such a scalar function is called a P-function, after L. E. Payne, who contributed significantly to this line of thinking. The challenge in engineering is to find a P-function that has practical significance.

Consider, for example, stationary stress analysis. The equilibrium equation is given by the system of equations

where σ_ij and f_i are components of the stress tensor and body force per unit mass, respectively.

In the case of linear elasticity, the stress and strain are related by the constitutive relation

where \epsilon_{ij} are components of the strain tensor, which in turn are related to displacements u_i by the strain-displacement relationship

If the material properties λ and μ are constant and the body force f is solenoidal (divergence free), such as the case with gravitational force in nonrelativistic applications, we can show that the volume change \epsilon_{kk}=\epsilon_{11}+\epsilon_{22}+\epsilon_{33}, also called the dilatation, satisfies both the maximum and minimum principles.

We can show this by substituting the constitutive relation in the equilibrium equation to get

and taking the divergence of the above vector equation to get

From the strain-displacement relationship and interchangeability of the order of differentiation, we see that \epsilon_{ij,ji}=\epsilon_{jj,ii}. This finally leads to

Since the dilatation satisfies the Laplace equation, both its maximum and minimum must be achieved at the boundary, unless the dilatation is uniform over the object. Taking the trace of the constitutive relation reveals that the same is true for the mean stress p=\frac{1}{3}\sigma_{kk} under the assumption of linear elasticity, homogenous material properties, and solenoidal body forces.

While the above result is theoretically appealing, engineers are rarely interested in finding the extreme values of the dilatation or mean stress. On the other hand, the maximum von Mises stress and the maximum shear stress are scalars of significant importance. For example, elastic limits (yield criteria) are often given in terms of either of these two quantities. Unfortunately, we do not have broadly applicable maximum principles for these and other quantities of practical interest. But we can still derive maximum principles for special cases.

A good example is the analysis of homogenous beams under torsion. If we pick the z-axis to be the axis of the beam, the stress analysis problem can be reformulated in terms of the stress function \phi(x,y), defined over the cross section Ω. In the absence of body forces, the stress function satisfies

and the nonzero components of stress are

where θ is the angle of twist per unit length of the beam.

Torsional stress distribution calculated with the Beam Cross Section interface.

We are interested in the maximum shear stress. Note that since the torsion problem does not have normal stresses, the magnitude of the shear stress is proportional to the von Mises stress. For magnitude τ of the shear stress, we have

Let’s take the divergence of |\nabla\phi|^2 to investigate the possibility of a maximum principle. Using the index notation

If we exchange the order of differentiation in the last term of the above equation, we get

The term in parentheses is the divergence of the stress function. According to the equation we started from, this value is -2. Thus, the second right-hand term above is zero. The first right-hand term is nonnegative because it is a sum of squares of components of the Hessian (matrix of second partial derivatives) of \phi. Additionally, it cannot be zero because the diagonals of the Hessian add up to a nonzero value, -2.

As a result, we have

This means that the magnitude of the shear stress achieves its maximum value at the boundary of the cross section of the beam. Many textbooks on linear elasticity contain analytical solutions to the Saint-Venant torsion problem for simple cross sections and boundary conditions. The explicit solutions there show that the maximum shear stress is encountered at the boundary. What we have done here is arrive at that conclusion for arbitrary cross sections and boundary conditions, without explicitly solving the boundary value problem. See the comparison below for an I-beam subjected to torsion.

Comparing the locations of the overall maximum shear stress (left) and boundary maximum shear stress (right) in a beam under torsion.

How to Exploit Maximum Principles in the COMSOL Multiphysics® Software

Say you are interested in checking extreme values only on boundaries of your problem. This may be because you have a maximum principle for the quantity of interest or you simply want to focus on the boundary. What functionality does COMSOL Multiphysics have to help you do so?

Maximum and Minimum Evaluations in Postprocessing

If you go to Results > Derived Values from the Model Builder or the menu in the GUI, you find maximum and minimum operators defined on volumes, surfaces, or lines. When you are interested in extreme values on the boundaries, you can use the surface or line operators as shown below, depending on the spatial dimension of your problem.

Maximum and minimum evaluations can be done on different geometric entity levels.

Maximum and Minimum Component Coupling Operators

The use of the above operators is limited to postprocessing. If you have to use the extreme values in your physics, say to implement a controller, you can define maximum or minimum component coupling operators in the Definitions section of your model. These operators also provide you with a choice of geometric entity between volumes, surfaces, curves, or points. After adding the operator, choose a geometric entity level in the settings window as shown below.

Geometric entity selection for component coupling operators.

In addition to saving time and memory because of fewer degrees of freedom, the boundary coupling operators cause lesser disturbance to the structure of the matrices involved than domain coupling operators. When you know that the equation you are solving has a maximum or minimum principle and you need to use an extreme value in your problem formulation, COMSOL Multiphysics offers you an economic alternative of using boundary coupling operators over domain coupling operators.

Storing Solutions on Selected Parts of Your Geometry

In a time-dependent problem where you are interested in postprocessing only boundary data, it can be economical or necessary — when you have bigger models and memory constraints — to store solutions from only the critical parts of your geometry. There is more than one way to do this in COMSOL Multiphysics, which we discuss in a previous blog post.

Surface Fatigue

If you are performing a fatigue analysis and you know that the critical behavior is at the boundaries of your geometry, COMSOL Multiphysics provides you with fatigue analysis features limited to just boundaries.

Fatigue analysis on domains vs. boundaries.

The Boundary Element Method

COMSOL Multiphysics is largely based on the finite element method. For some problems, it supports the boundary element method. In the latest version of the software (version 5.3), the boundary element method is available for electrostatics, corrosion, and user-defined equations.

Problems solved with the boundary element method have no source terms. As such, these problems are likely to obey both the maximum and minimum principles. Exploiting the maximum principles is not a legitimate motivation for using the boundary elements method. But when you can use the method, the problem most likely has extreme values on boundaries.

Caveats to Keep in Mind

The maximum principles are very useful when applied to problems where their premises apply. However, we should refrain from applying them to problems where they are not applicable. Discontinuities and changes in the sign of the source terms, for example, are flagrant violations of these premises.

Discontinuities

In our mathematical analysis, we assumed smoothness of the coefficients. If there are discontinuities, the principles do not apply to the whole domain. However, it is possible to subdivide the domain into regions with no discontinuities and apply the maximum principles to each region. Practically, this means you have to include internal boundaries as well as external boundaries. The interface between two materials with different material properties is one example.

We mentioned earlier that in COMSOL Multiphysics, fatigue analysis can be done either on domains or on surfaces. In many problems, fatigue failure starts at the surface. However, in contact mechanics, it is known that the fatigue failure starts in subsurface points. For more on this, including the implications of the finite element meshing strategy, please see the following blog post on modeling contact fatigue.

Changing Signs of the Source Term

In problems with nonzero source terms, the maximum principles require nonnegativity or nonpositivity. As such, they do not apply if the source term changes signs in the domain. For example, in a Joule heating problem, the electromagnetic heat source is always positive. In a chemical reaction, on the other hand, the same chemical species can be produced in parts of the domain and consumed in other parts. Unless we are sure that the reaction will be producing a species everywhere or consuming a species everywhere, we cannot use the maximum principles.

Effects of Discretization

Our derivation of the maximum principles uses the partial differential equation directly. It is possible that discretizations in the finite element or other numerical methods destroy this property. If a problem you expect to obey one of the maximum principles does not exhibit that when you solve it numerically, check if the numerical method obeys the so-called discrete maximum principles. Keep numerical errors in mind as well.

Concluding Thoughts on Using Maximum Principles

The maximum principles are traditionally used in the analysis of partial differential equations for existence/uniqueness type proofs, establishing error bounds, and other similar ends. In this blog post, we have shown how to use these results for efficient computational analyses.

Note that the mathematical discussion provided in the first section of this post was meant for a broader audience. For a more rigorous discussion, please refer to standard texts in the analysis of partial differential equations. Here are my personal favorites for this subject:

• Lawrence C. Evans, Partial Differential Equations, American Mathematical Society, 2010
• Rene Sperb, Maximum Principles and Their Applications, Academic Press, Inc., 1981

If you have any questions related to this topic or using COMSOL Multiphysics, please contact us.

Further Reading

• Browse more blog posts related to postprocessing functionality in the COMSOL® software
• Learn about new numerical studies and solvers in COMSOL Multiphysics version 5.3