I'm trying to solve the steady state of a heat equation problem in 2D $\Delta u = 0$ (3D also), with the method of solving the huge system of equations that arises from the discretization of the domain. My boundary conditions are different from any example that I've been able to find.

If my domain is a 2D square, I want to put a fixed temperature (say $T_e$) all around the square edges. However, I want to force a point in the center of the square to be at a different temp from the edges, $T_c, T_c > T_e$. Think of a square metal plate being heated at the center by a soldering iron. I want to calculate the steady state heat in such a scenario.

To the best of my knowledge it seems that this problem is not really a heat equation problem: having a point at the center of the domain with a local maxima temperature goes against the properties of the heat equation...

I'm obtaining results, but it seems that heat does not propagate linearly from the center to the edges, but rapidly decaying. The problem is the same if moved into the 3D domain (a sphere/cube that is heated from it's center to its exterior edges, and my results are the same.

Maybe I'm using a wrong model? or a wrong tool to solve this steady state problem? Any hints or directions will be highly appreciated.

Best regards

  • $\begingroup$ What makes you think the the decay should be linear? The fundamental solution of the Laplace equation in 2D behaves as $-\log |x|$ (and as $|x|^{-1}$ in 3D). $\endgroup$ Nov 29, 2012 at 15:43
  • $\begingroup$ I said linear because when using more standard boundary conditions ($T_1$ in one side of the square, and $T_2$ in the opposite side) the final heat distribution between the sides decays linearly. By fundamental solution you mean applying heat at only one point (in the boundary)? $\endgroup$
    – asmatic
    Nov 30, 2012 at 13:42
  • $\begingroup$ Yes, by an internal source at only one point. $\endgroup$ Nov 30, 2012 at 13:55
  • 2
    $\begingroup$ You can't expect a linear decay because you are injecting heat at only one point and you are taking it out along an entire, extended boundary. There is no reason to expect that this should lead to a linear decay -- in fact, I think it should be easy to see that it shouldn't be linear since the heat needs to spread out to an ever larger area as you move from the source to the boundary. $\endgroup$ Nov 30, 2012 at 14:27

4 Answers 4


You need to add the condition that $u(c) = T_c$ as an explicit constraint via a Lagrange multiplier. For simplicity, assume that $T_e=0$ (otherwise construct a $u_e$ such that $u_e|_{\partial\Omega} = T_e$ and add it to the solution with homogeneous boundary conditions, shifting the constraint by $T_c-u_e(c)$).

You know that the solution of $-\Delta u = 0$ solves the minimization problem $$ \min_{u\in H^1_0(\Omega)} \int_\Omega |\nabla u|^2 \,dx. $$ To enforce the constraint $u(c)=T_c$, you consider the Lagrange functional $$ \max_{\lambda\in\mathbb{R}} \min_{u\in H^1_0(\Omega)} \int_\Omega |\nabla u|^2 \,dx + \lambda(u(c)-T_c). $$ The necessary optimality conditions for this problem are obtained (formally) by differentiating with respect to $u$ and $\lambda$: \begin{align} \int_\Omega \nabla u\cdot \nabla v + \lambda v(c)&= 0 \quad\text{for all } v\in H^1_0(\Omega)\\ (u(c) - T_c)\mu &=0 \quad\text{for all }\mu\in\mathbb{R}. \end{align} (Note here that except in 1D, functions in $H^1(\Omega)$ are not continuous and thus do not have point values. You probably should replace the point value with an average value over a small subdomain. For a finite element discretization, you move to a subspace of piecewise polynomials, where point values make sense.)

For details on saddle point problems for partial differential equations, see Dietrich Braess, Finite Elements, Cambridge University Press, 2007, Chapter III.4.

  • 3
    $\begingroup$ This is a good answer, but a simpler alternative is preferable for use with common discretizations: instead of using a Lagrange multiplier, build the condition $T(c) = T_c$ into the discrete solution space, eliminating the associated degree of freedom (either by "lifting" or by algebraically decoupling it). This will retain the SPD (and M-matrix, if applicable) properties of your original discretization, thus simplifying the solution methods. $\endgroup$
    – Jed Brown
    Nov 29, 2012 at 15:49
  • $\begingroup$ True; this is very focused on Galerkin discretizations (and slightly overkill even for those). For finite differences (and finite elements if it is convenient to have $c$ as a node), your suggestion (which should probably be an answer) is the way to go. $\endgroup$ Nov 29, 2012 at 16:07

From a mathematical viewpoint you can't do what you are trying to do: the problem $$ \Delta u = 0 \qquad \text{in } \Omega, \\ u = T_e \qquad \text{on } \partial\Omega, \\ u(c) = T_c \qquad \qquad\qquad\text{} $$ does not, in general, have a solution. This is because, mathematically speaking, the space of solutions of the Laplace equation, $H^1$, consists of functions that may not be continuous and so it doesn't make sense to take point values such as $u(c)$.

An alternative explanation of this problem would be to realize that if $T_c>T_e$ then there will be heat flowing from $c$ to the boundary, and consequently someone will have to put that heat into the body at $c$. But putting a finite amount of heat into a single point (of area zero) means to have a device that can produce an infinite heat flux -- hardly realistic.

The solution to what you are trying to do is to assume that the area at which you are able to impose a temperature $T_c$ is finite, i.e., to consider a domain that has a hole of finite size at the center.


Based on some of the previous answers, I think you are faced with a few choices in your modeling. If you are forced to use a cartesian unskewed grid (not clear from your question what discretization you are using), it's going to be a little tougher because as said previously, you need to cut off a finite area of a few cells in the center and you might get corner issues once you try to solve this.

If you can do arbitrary grids, may I suggest a circular one (not a square) which means you can use the symmetry properties to reduce the dimensionality of your problem. You then can compare your results with an analytical solution of the form:

$T(r) = (T_1 - T_0)\frac{\ln{\frac{r}{r_0}}}{\ln{\frac{r_1}{r_0}}}+ T_0$

And you can see that as $r_0 \rightarrow 0$, the temperature at the origin becomes undefined because your heat flux (positive or negative) would become infinite.

One solution

  • $\begingroup$ Currently we are working with Cartesian grid. In fact I work over the voxels of binary volumes. The example I gave was indeed simplified. I think I would prefer to stay away from highly specific grids as this should work for arbitrary volumes / shapes. $\endgroup$
    – asmatic
    Dec 11, 2012 at 9:27
  • $\begingroup$ But what I was trying to say was that if you are looking for a validation case for your method, you can use that specific grid and that analytical solution to compare your results. That should tell you if you are on the right track. $\endgroup$ Dec 11, 2012 at 12:51

I am not sure about interpreting $\nabla u=0$ as heat spreading, given that it corresponds to the steady state of the time dependent (parabolic) equation $u_t=\nabla u$. For me the meaning is rather closer to an interpolation of the values specified at the boundary of the domain of definition $\partial \Omega$. These values definitely determine the function $u$, so if they are set to zero on a ball of radius 1, the function is identically zero. This is consistent with the formula of the solution to $\nabla u=0$ with $u=f$ at the boundary, given by:

$u(x)=\frac{-1}{\omega_n} \int_{|s|=1} \frac{-1+|x|^2}{|s-x|^n} f(s) ds, x \in \Omega$

for $| \cdot |$ the norm of the Euclidean space.

Moreover, having an interior point with maximum value is against the maximum principle holding for all armonic functions.

Thus, I am not sure if a boundary problem on a unit ball with constant temperature at its boundary and another temperature at an interior point it is well formulated for $\nabla u=0$. At least the interior point should be removed, so the domain becomes a sort of spherical band.

In other words, given that the solution to $\nabla u=0$ depends on the function values set on the boundary of the domain and on the boundary shape itself, what are these in the case of the fundamental solution (given by a radial solution)?. Further, is there any categorization of generic (fundamental) solutions depending on boundary shapes and values?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.