I'm trying to solve very simple differential equation with perdiodic boundary conditions and non-constant initial condition, however the results I'm getting don't seem to be correct. Here's my code:

size = 10 (* size of the domain *)
tMax = 20000000 (* time of the experiment *)
pde = D[h[x, y, t], t] == -0.7
BConds = {h[0, y, t] == h[size, y, t], h[x, 0, t] == h[x, size, t]} (* periodic boundary conditions *)
topography[x1_, x2_] = Piecewise[{{400./size*x1 - 100., size/4. < x1 < size/2.}, {300. - 400./size*x1, size/2. < x1 < 3.*size/4.}, {100., x1 == size/2.}, {0, x1 <= size/4.}, {0, x2 >= 3.*size/4.}}]
ICond = h[x, y, 0] == topography[x, y] (* initial condition as a function of x *)
sol = NDSolve[Join[{pde}, BConds, {ICond}], h, {x, 0, size}, {y, 0, size}, {t, 0, tMax}]
fun[x_, y_, t_] := h[x, y, t] /. sol[[1]]

I wound expect, because of how the initial condition is specified, the solution to be symmetrical with respect to size/2. That's true for e.g. x=size/4 and x=3*size/4, meaning .fun[size/4,0,0]=fun[3*size/4,0,0].

However that's not true for x=size/3 and x=2*size/3, for which the solution in my opinion also should have the same values, but it doesn't:

fun[size/3, 0, 0]
fun[2*size/3, 0, 0]

Why the values aren't the same?

The other problem is when I try to plot the solution:

Manipulate[ContourPlot[fun[x, y, t], {x, 0, 10}, {y, 0, 10}, ColorFunction -> (ColorData["Rainbow"][1 #] &), PlotLegends -> Automatic], {t, 0, tMax}]

The plot, for t=0, displays regions of fun[x,y,t]<0 even though the minimum value for the initial condition is 0. I would guess it's because of how the solution is interpolated, but could someone explain it to me in more details?

  • $\begingroup$ I tried to use Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> "FiniteElement"}} and it gives much better results, but when I apply periodic boundary conditions the execution fails saying that FEM doesn't support that. $\endgroup$
    – lgaza
    Dec 2, 2015 at 21:15

1 Answer 1


If you change the discretization points to be odd, then you get a better result.

sol = NDSolve[Join[{pde}, BConds, {ICond}], 
  h, {x, 0, size}, {y, 0, size}, {t, 0, tMax}, 
  Method -> {"PDEDiscretization" -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> 35, "MinPoints" -> 35}}}]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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