# Solving a non-linear heat equation with the galerkin method gives negative values

I am trying to solve a non-linear time-dependent heat equation $$\partial_tT=\nabla \left(k_T(T)\nabla T\right) + f$$ using the galerkin method, with neumann boundary conditions. For linearization of the nonlinear part I am using the newton approach. The source term $f$ always stays positive, thus I would expect that $$T\geq T_0\forall t\geq t_0$$ but now I notice that at some nodes my value for $T$ drops below the initial value, and even going into negative ranges. I do not understand why that happens, after it is independent of the time step, but not independent of the value for $f$.
What could be the reason for this behaviour?
I notice that the range for the negative values is approximately five to six orders smaller than the maximum value for $T$ at that time, thus could it be an accuracy problem?

• What timestepping scheme are you using? A common mistake I've seen is using the Crank-Nicholson scheme for problems where preserving positivity is important. The CN scheme, while 2nd-order accurate, doesn't have a discrete maximum principle. Implicit Euler is less accurate, but it's a better choice where this property is concerned. – Daniel Shapero Jun 4 '18 at 18:09
• Implicit euler, I tested with Crank-Nicholson, but was not satisfied with the results – arc_lupus Jun 5 '18 at 7:09

## 1 Answer

It turns out that a variation of this question was recently asked on the deal.II mailing list, so I'll copy my answer here:

Your misunderstanding is that you believe that a property of the exact solution of the equation is supposed to also hold for the discrete one. But there is really no reason for this to be the case other than the fact that you know it must be true for the limit of your numerical solutions as $h\rightarrow 0$ in some sense.

To give you just one example, if you take the Heaviside function $$H(x) = \left\{ {0 \qquad\text{if x\le0}} \atop {1 \qquad \text{if x>1}} \right.$$ then it is obviously non-negative. But if you take the Fourier transform and truncate after a finite number of basis functions, then you get Gibbs phenomenon -- there are numerous examples of the over- and undershots here: https://en.wikipedia.org/wiki/Fourier_transform So the continuous property is also not satisfied for the finite-dimensional case. In fact, it is not even satisfied in a pointwise sense if you let $N\rightarrow \infty$: the Fourier transform only converges in $L_2$, not in $L_\infty$!

• I will get a problem if other calculations are done using that result which will fail if $T$ goes into the negative region. How bad can it be (in terms of "what kind of error will I get") when limiting the calculations to certain values, i.e. resetting every value below $T_0$ to $T_0$? – arc_lupus Jun 5 '18 at 7:18
• I understand why you don't want $T$ to be negative. There are similar problems when you solve hyperbolic problems and the density becomes negative, for example. You then have two options: Either you use a (generally much more complicated and nonlinear) method that guarantees a positive solution, or you just clip to $T_0$ as you suggest. How large an error that introduces depends on how exactly you do it and what equation you solve; I would expect the error introduced by this procedure to be at least $O(h^2)$. – Wolfgang Bangerth Jun 5 '18 at 15:18
• Which method would you suggest? Thanks! Concerning the error: At least means $>O(h^2)$ or $<O(h^2)$? – arc_lupus Jun 5 '18 at 15:20
• What method to suggest depends on what you want to do. I don't think there is a general answer. In hyperbolic methods, positivity-preserving methods are generally nonlinear even if the equation is linear, with all attendant problems. I suspect that the same would be true for your equation. – Wolfgang Bangerth Jun 7 '18 at 2:08
• "I would expect the error to be at least $O(h^2)$" -- I mean this to be "I would expect the error to be at least $O(h^s)$ with $s\le 2$". – Wolfgang Bangerth Jun 7 '18 at 2:09