NDSolve issue with initial and boundary conditions

Question

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?

Answer 1

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:

Show[
 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}]

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

Answer 2

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

Method->{"PDEDiscretization"->{"MethodOfLines",{"SpatialDiscretization"->"FiniteElement"}}}