I solved the 2D heat equation using the finite element method. It all went well first with the adiabatic case, however problems occured when I introduced cooling with the enviroment. I modeled the process with newton's law of cooling: $$ q = h (T - T_{\infty}) $$ Using the term $$ q = k \nabla T $$ Assuming that $k=1$: $$ \nabla T = h(T - T_{\infty}) $$ For the heat equation I used $u$ as the function for temperature $$ \int{v\Delta v}dV = \oint{v\nabla u \cdot dS} - \int \nabla u \cdot \nabla v dV $$ Ending up with the boundary condition $\frac{\partial u}{\partial n} = h(T_{\infty} - u) $ During the testing I realized that this only works if $T_{\infty}$ is zero. I suspect that this is due to the fact that if $T_{\infty}$ is not zero, the right hand side is no longer bilinear. Am I on the right track here and if so what are ways to fix non linearity while using the fem?


Here is the equation which I am trying to solve: $$ \Delta u = \frac{du}{dt} $$ Multiplying with a test function and applying greens identity: $$ \oint{v\nabla u \cdot dS} - \int (\nabla u \cdot \nabla v) dV = \frac{d}{dt}\int (uv) dV $$

  • $\begingroup$ The equation you wrote, $\partial_t u = \Delta u$, is just the standard heat diffusion equation. This is a linear equation and is discretised in space as $\partial_t M u = -W u$, where $M$ is the mass matrix and $W$ is the stiffness matrix. You also didn't specify any boundary conditions. I am also unclear why you are replacing $\nabla u$ with $h(T_{\infty} - u)$, I am assuming that your boubdary condition is supposedly $n\cdot \nabla u + h u = h T_{\infty}$. You'll have to clarify these points. $\endgroup$
    – lightxbulb
    Mar 19 at 19:57
  • $\begingroup$ My thought process was the heat leaving the body is $h(T_{\infty} - u)$ by newtons law of cooling. A small change in temperature times the heat capacity is the heat added $c * dT = q$. So in my case it would result in $dT = h(T_{\infty} - T)$ where T is just the function u. The idea was to replace the gradient at the boundary with $dT$. If that doesn't work let me know. $\endgroup$
    – Boiler4562
    Mar 19 at 20:15
  • $\begingroup$ I need an initial boundary value problem formulation. Your current formulation is missing boundary conditions. $\endgroup$
    – lightxbulb
    Mar 19 at 20:18
  • $\begingroup$ Yes, my boundary condition is $n \cdot \nabla u + hu = hT_{\infty}$ and I thought that the way of applying it is plugging it into the integral over the surface of the domain $\endgroup$
    – Boiler4562
    Mar 19 at 20:40
  • $\begingroup$ This question is hard to read (as evidenced by the several questions in comments and answers) in part because it appears not to have been proof-read. For example, in the first equation the variable in question is $T$ whereas in the others the solution is $u$. You never explain what $q$ is supposed to be, and it doesn't ever appear again in the rest of the post. In the equation that shows the integration by parts, you have $v$ twice on the left when one should have been $u$. Sometimes, multiplication uses no sign at all and at other times it is $*$. Please make it easier to understand. $\endgroup$ Mar 20 at 17:15

2 Answers 2


From your formulation it's still unclear what are the boundary conditions. But from the sentence about replacing $n\cdot \nabla u$ with $h(T_{\infty} - u)$ and your comment it sounds like your boundary conditions are $hT_{\infty} u + n\cdot \nabla u = hT_{\infty}$ so I will just handle the general case.

Suppose you have the initial boundary value problem: \begin{alignat}{2} \partial_t u(t,\vec{x}) &= \Delta u(t,\vec{x}), &\quad& \vec{x}\in\Omega, \, t\in (0,\infty) \\ \alpha(\vec{x})u(t,\vec{x}) + \partial_n u(t,\vec{x}) &= g(\vec{x}), &\quad& \vec{x}\in\partial\Omega, \, t\in(0,\infty), \\ u(0,\vec{x}) &= f(\vec{x}), &\quad& \vec{x}\in\Omega. \end{alignat}

Multiplying the first equation by a test function $v$ and using the divergence theorem leads to:

\begin{equation} \int_{\Omega} v \partial_t u = \int_{\Omega} v\Delta u = -\int_{\Omega}\nabla v \cdot \nabla u + \int_{\Omega} \nabla \cdot (v\nabla u) = -\int_{\Omega}\nabla v \cdot \nabla u + \int_{\partial\Omega} v\partial_n u. \end{equation}

You can substitute the boundary condition:

$$\int_{\partial\Omega} v\partial_n u = \int_{\partial\Omega} v (g-\alpha u).$$

Then reordering some terms you get:

$$\int_{\Omega} v \partial_t u +\int_{\Omega}\nabla v \cdot \nabla u + \int_{\partial\Omega} v \alpha u = \int_{\partial\Omega} v g.$$

If $v$ is from a finite-dimensional space then it can be written as a linear combination of some basis vectors $\psi_1,\ldots,\psi_n$, similarly if $u$ is from a finite-dimensional space then it can be written as a linear combination of some basis vectors $\phi_1,\ldots,\phi_n$. Setting $v$ to $\psi_1,\ldots,\psi_n$ consecutively and $u(t,\vec{x}) = \sum_{j=1}^n a_j(t) \phi_j(\vec{x})$ you get for $1\leq i \leq n$: \begin{gather} \sum_{j=1}^n \left(\int_{\Omega} \psi_i \phi_j\right) \partial_t a_j(t) + \sum_{j=1}^n \left(\int_{\Omega} \nabla \psi_i \cdot \nabla \phi_j\right) a_j(t) + \sum_{j=1}^n \left(\int_{\partial\Omega} \alpha \psi_i \phi_j\right) a_j(t)= \int_{\partial\Omega} \psi_i g. \end{gather}

The mass matrix is $M_{ij} = \int_{\Omega} \psi_i\phi_j$, the stiffness matrix is $W_{ij} = \int_{\Omega}\nabla \psi_i \cdot \nabla \phi_j$, the boundary correction matrix is $B_{ij} = \int_{\partial\Omega} \alpha \psi_i \phi_j$, the right-hand side is $b = \int_{\partial\Omega} \psi_i g$. Then you have: $$M\partial_t a + W a + B a = b \implies \partial_t a = -M^{-1}(W+B)a + M^{-1}b.$$

Formally the solution is:

$$a(t) = \exp(-tM^{-1}(W+B))a(0) + \int_0^{t} \exp((t-s)M^{-1}(W+B)) b\,ds,$$

although the exponential may not be easy to compute (the integral can be rewritten so that it disappears and you only have a matrix exponential but I will not go into that unless you need it).

You can instead apply e.g. explicit Euler to the above: \begin{align} \frac{a^{k+1}-a^k}{\tau} = -M^{-1}(W+B)a^k + M^{-1}b \implies a^{k+1} = (I-\tau M^{-1}(W+B))a^k + \tau M^{-1}b, \end{align}

or implicit Euler: \begin{align} M\frac{a^{k+1}-a^k}{\tau} = -(W+B)a^{k+1} + b \implies a^{k+1} = (M+\tau (W+B))^{-1}Ma^k + \tau b. \end{align}

With implicit Euler you can take arbitrary large step sizes and it will remain stable (if the matrices $(W+B)$ and $M$ are positive semi-definite). Generally you will also have that $\psi_i = \phi_i$ which would yield symmetric positive semi-definite matrices (I am not sure how the $B$ affects that, maybe there are some conditions on the $\alpha$) allowing you to use the conjugate gradient method.

In either case from $u(0,\vec{x}) = f(\vec{x})$ you also get:

$$\int_{\Omega} v u = \int_{\Omega} v f \implies a(0) = M^{-1} \begin{bmatrix} \int_{\Omega} \psi_1 f \\ \vdots \\ \int_{\Omega} \psi_n f\end{bmatrix}.$$

So you know what $a^0$ is. Also one often makes $M$ diagonal, also known as a lumped mass matrix (see this post https://scicomp.stackexchange.com/a/19714/34261), this can for instance save resources in the inversion if you use explicit Euler. If you use lumping there you will not need to solve linear systems. On the other hand, for explicit Euler you will have some constraints on the step size $\tau$ in order to guarantee stability, e.g. $\tau < 2/\rho(M^{-1}(W+B))$.

  • $\begingroup$ So oversimplified the solution is to split the integral over the boundary into two seperate ones? And what is missing for the boundary conditions when I say that the gradient in normal direction is equal to something? I thought specifying the the values or gradient at the border is all that is needed. $\endgroup$
    – Boiler4562
    Mar 19 at 21:17
  • $\begingroup$ @MatthiasFriedel You want to have the unknowns all on one side, that's why we split the integral yes. What was missing in your question was a proper formulation, see for example how I formulated the initial boundary value problem - that's what people expected to see in order to be able to help you. Everything starts with your (initial) boundary value problem, all weak formulations and discretisations come after that. $\endgroup$
    – lightxbulb
    Mar 19 at 21:20
  • $\begingroup$ @MatthiasFriedel Also probably worth noting, but you can't always plug in boundary conditions, e.g. if you had Dirichlet boundary conditions you would have had to handle those in a special manner. $\endgroup$
    – lightxbulb
    Mar 19 at 21:33

You can move the term with $h T_\infty v$ to the right hand side so it becomes part of the linear form and the load vector.

  • $\begingroup$ Maybe I misunderstood your answer, however I try to solve the non - stationary heat equation, that means I have a bilinear form on both sides. The equation I posted in my question is only one side of the equation. $\endgroup$
    – Boiler4562
    Mar 19 at 15:09
  • 1
    $\begingroup$ Please post the equation you are trying to solve. $\endgroup$
    – knl
    Mar 19 at 17:10
  • $\begingroup$ What @knl is pointing out is that if you have inhomogenous Neumann or Robin boundary conditions, you end up with a right hand side term that leads to a linear (not a bilinear) form. $\endgroup$ Mar 20 at 17:17

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.