Solving Laplace’s Equation

(This article is the 4^th in a series on differential equations. The previous installments were 1) Understanding the Heat Equation, 2) Understanding the Wave Equation, and 3) Understanding Laplace’s Equation.)

The previous article explained the meaning of Laplace’s Equation. We saw a couple examples of equilibrium heat distributions in a piece of heat-conducting material, including this one:

Those illustrations correspond to solutions of Laplace’s. Now, how can we find such solutions?

First, when one tries to find a solution to Laplace’s Equation, what does that actually mean? The equation describes the relationship between the derivatives of a (perhaps unknown) function. Solving the equation means finding a function whose derivatives satify the equation; but it is not very interesting to find just any old function which happens to fit. Usually we want to find one specific function which corresponds to some physical situation of interest; perhaps, as in the examples, it might be a function which gives the temperature at various points on some specific object. Or it might be a function which gives electrical potential or fluid flow at various points in space.

So say we have a particular physical situation in mind. We don’t know, but want to know, how some parameter (temperature, electrical potential, etc.) will vary over space (in that situation), and believe that its distribution will satisfy Laplace’s Equation. What do we need to know about that physical situation to identify the one specific function which gives the desired distribution?

Remember that Laplace’s is a second-order differential equation. When searching for a specific solution to a 2^nd-order equation, we generally need the value of the unknown function and its first derivative(s) at one or more points. (These are called boundary conditions.)

So, when solving Laplace’s, you need to know the value of the unknown function at some points. For concreteness, say the unknown function represents temperature; then we need to know what the temperature will be somewhere on the object in question. If certain parts of the object will be “held” at a specific temperature, we need to know exactly which parts those are and what their respective temperatures will be.

(You might wonder: Why would some parts of the object be held at a specific temperature? Take the example of the metal bar thrust into a fire, which was mentioned in an earlier article in this series. Presumably, the entire end which is sitting in the fire will be at the same temperature as the fire. Often, there are similar constraints which allow you to estimate the value of the unknown parameter at some locations.)

The details of the physical situation may constrain the first derivative(s) of the unknown function as well. When solving for a temperature distribution, the simplest such constraint is to imagine, as an approximation, that the surface of the object is perfectly insulated and that no heat will enter or leave through the surface. That means, at every point on the surface, the directional derivative in a direction tangent to the surface will be zero.

So we need to find a function which satisfies Laplace’s, we know the value of the function for some parts of its domain, and we may also have some information about its first derivative(s), at least for some parts of its domain. In some cases, that is enough to solve the problem analytically, meaning you can find a formula which defines the target function. (I won’t tell you how to find analytic solutions; that is out of scope for this article.)

However, the unfortunate reality is that situations where an analytic solution can be found are usually contrived. If you are solving a real-life engineering problem, chances are high that you should completely skip looking for analytic solutions and use a computer to generate a numerical solution instead. This means you won’t get a formula which defines the target function, but will just get its value (or a close approximation) at certain points. The computer will not arrive at those values by algebraic reasoning, but will calculate them using an iterative algorithm.

To start computing a numerical solution, the first step is to discretize the problem, just like we did in the earlier article on the Heat Equation. We can divide space up into a grid, and also divide time into steps with a specific length.

In the earlier article, we learned that the Heat Equation implies each tiny piece of material will move towards the average temperature of its neighbours. For Laplace’s Equation to hold, meaning that none of the temperatures are changing, the temperature at each point must equal the average temperature of its neighbours. That is enough for us to express the discretized problem as a system of linear equations.

Here’s an example. Take a square, two-dimensional space, like the “plate” examples in the previous article, and discretize it into a grid with just four (arbitrarily-numbered) points:

For the temperature at each of those four points to equal the average of the neighbouring points, these linear equations must hold:

$$T_1=\frac{1}{2}(T_2+T_4)$$ $$T_2=\frac{1}{2}(T_1+T_3)$$ $$T_3=\frac{1}{2}(T_2+T_4)$$ $$T_4=\frac{1}{2}(T_3+T_1)$$

We can rewrite those as:

$$2T_1-T_4-T_2=0$$ $$2T_2-T_1-T_3=0$$ $$2T_3-T_2-T_4=0$$ $$2T_4-T_3-T_1=0$$

But better yet is to convert those equations to a single matrix equation:

$$\left[\begin{array}{cccc}2& -1& 0& -1\\ -1& 2& -1& 0\\ 0& -1& 2& -1\\ -1& 0& -1& 2\end{array}\right]\left[\begin{array}{c}{T}_1\\ T_2\\ T_3\\ T_4\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$$

We have four unknown variables and four equations, which gives some hope that we can solve for the unknown values, perhaps using a standard method like Gaussian Elimination. However, if you check the determinant of that matrix, you will find that:

$$\det \left[\begin{array}{cccc}2& -1& 0& -1\\ -1& 2& -1& 0\\ 0& -1& 2& -1\\ -1& 0& -1& 2\end{array}\right]=0$$

...which means the linear equations do not have one unique solution; they must have either infinite solutions, or no solution. And in fact, every vector of the form $\left[\begin{array}{cccc}c& c& c& c\end{array}\right]$ is a solution. In physical terms, that means: the entire object will be at the same temperature, but we don’t know what temperature that is.

Which makes sense! I told you earlier that to find a unique solution, we generally need to know what the temperature will be at one or more points. We did not provide any known temperatures as inputs. (We did implicitly assume that the plate is perfectly insulated, because our linear equations did not include any terms for heat entering from or escaping to the outside.)

What if we say that the temperature at one point is known; for example, maybe it could be that $T_1=100$? Then we can drop the equation for $T_1$, and the remaining three equations can be written as:

$$2T_2-T_3=100$$ $$2T_3-T_2-T_4=0$$ $$2T_4-T_3=100$$

Or:

$$\left[\begin{array}{ccc}2& -1& 0\\ -1& 2& -1\\ 0& -1& 2\end{array}\right]=\left[\begin{array}{c}{T}_2\\ T_3\\ T_4\end{array}\right]\left[\begin{array}{c}100\\ 0\\ 100\end{array}\right]$$

That matrix equation does have a unique solution: $\left[\begin{array}{ccc}100& 100& 100\end{array}\right]$. Which, again, makes sense! We have a perfectly insulated plate, and one corner is held at 100°, while the temperatures of all other parts are allowed to vary naturally. It’s fairly obvious that equilibrium will be attained when the entire plate warms up to 100°.

Let’s just stay with this example a little bit more. What will happen if we know the temperature for two opposing corners? Say $T_1=100,T_3=0$? (That could mean that one corner of the plate is submerged in boiling water, while the opposite corner is submerged in ice.)

We can drop the equation for $T_3$, leaving us with:

$$\left[\begin{array}{cc}2& 0\\ 0& 2\end{array}\right]\left[\begin{array}{c}{T}_2\\ T_4\end{array}\right]=\left[\begin{array}{c}100\\ 100\end{array}\right]$$

OK, so each of the remaining corners, which are halfway between the hot corner and the cold corner, will reach equilibrium at 50°. That’s just what we would expect.

Now here is the fun part. I’ll draw a grid below, which represents a plate of heat-conducting material. You can choose the grid size, i.e. how fine or coarse the discretization should be. Then click on any nodes that you like in the grid and enter their fixed temperatures. I will convert that to a matrix equation and solve for the equilibrium temperatures which all other nodes will reach through heat conduction.

What if the plate has a hole in the middle?

Thanks for reading! In the next installment of this series, I want to show you the algorithm which I am using to solve those matrix equations.