While solving the heat equation in one spatial variable $u_t = u_{xx} $ (x goes from 0 to L) with the initial temperature distribution $T_0 \frac{x(L-x)}{L^2}$ , and with neumann boundary conditions $u_x(0,t) = u_x(L,t) = 0$, I got some really weird behaviour from NDSolve.

My code looks like this:

h[x_] := x*(30 - x)/900;
pde = D[u[t, x], t] == D[u[t, x], x, x]
begin = 0;
end = 30;
bc = {u[0, x] == 100*h[x], (Derivative[0, 1][u])[t, begin] == 
    0, (Derivative[0, 1][u])[t, end] == 0};
finaltime = 100

s = NDSolve[{pde}~Join~bc, u, {t, 0, finaltime}, {x, begin, end}];

Since heat cannot flow out through the ends, continuing this in time should yield a smoothening until it reaches the average everywhere. Instead, I get a very weird time evolution which when plotted seems to be of the form $u_x(x,t) = u_x(x,0) - kt$. This is particularily infuriating because the problem seems to be intermittent. Taking the square of h does not cause any trouble.

The problem seems to magically fix itself if I instead feed in a truncated cosine series into the code:

rule = t -> FourierCosSeries[t*(2*Pi - t), t, 35];
f[x_] := (t /. rule) /. (t -> x)
g[x_] := f[2*Pi*x/30]/(4*Pi*Pi)

and insert g instead of h into the code above. Trying a function interpolation gave only errors.

Is there a fix which is more general than this quick hack?

  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$
    – bbgodfrey
    May 16, 2015 at 13:49

2 Answers 2


You can use the finite element method with the method of lines as @toadatrix suggested, but for the FEM method to work, you need to do a little more. The Neumann boundary conditions need to be specified using NeumannValue.

h[x_] := x*(30 - x)/900;
op = D[u[t, x], t] - D[u[t, x], x, x];
begin = 0;
end = 30;
bc = {u[0, x] == 100*h[x]};
neumann = NeumannValue[0, x == begin] + NeumannValue[0, x == end];
finaltime = 100;

s = NDSolve[{op == neumann, bc}, 
   u, {t, 0, finaltime}, {x, begin, end}, 
   Method -> {"PDEDiscretization" ->
    {"MethodOfLines", {"SpatialDiscretization" -> "FiniteElement"}}}

The default boundary condition is NeumannValue[0,...], so you could just simply use op == 0.

Check the result:

 Plot3D[u[t, x] /. s // Evaluate, {t, 0, finaltime}, {x, begin, end}, 
  MeshFunctions -> {#2 &}, PlotRange -> All, AxesLabel -> Automatic],
 Graphics3D[{Red, Opacity[0.5], 
   InfinitePlane[{0, 0, Integrate[100*h[x], {x, 0, 30}]/30}, {{1, 0, 0}, {0, 1, 0}}]}]
Show[%, ViewPoint -> {Infinity, 0, 0}]

Mathematica graphics

Mathematica graphics
The graphs shows that the temperature evens out to the average over time.


I believe that all you need to do is select the following method option in NDSolve


  • $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this this meta Q&A helpful $\endgroup$
    – Michael E2
    Aug 14, 2015 at 18:20
  • $\begingroup$ Did you try this? I get an error. (V10.2.) $\endgroup$
    – Michael E2
    Aug 14, 2015 at 18:22

Your Answer

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