Problem with boundary values for a partial differential equation

I have a PDE problem that I can’t solve. and I wonder if you could help me with it.

I have a Terzaghi Consolidation equation (similar to the heat equation) that is written as follows:


uuu is pore water pressure

zzz is depth

ttt is time

cc is a constant

pore water pressure is equal at all depth at t=0t = 0, but it becomes 00 at top and bottom directly after that.

u = u[z,t] (zz is depth and tt is time)

with boundary conditions as follows:

u[z, 0] = 100 (uu at all depth when t=0t = 0)

u[2, t] = 0 (uu at top of the layer; in this case, 22 is the height of the top layer)

u[0, t] = 0 (uu at bottom of the layer)

I have solved this kind of problem before using a finite difference technique in a spread sheet, I wanted to find out whether Mathematica can solve this problem directly, but unfortunately I failed.

This is what I entered into Mathematica:

NDSolve[{Derivative[0, 1][u][z ,t] == Derivative[2, 0][u][z, t],
u[0, t > 0] == 0, u[2,t > 0] == 0, u[2 > z > 0, 0] == 100}, u, {z, 0, 2},{t, 0, 100}]

Unfortunately, evaluating this expression leads gives the error:


“To avoid possible ambiguity, the arguments of the dependent variable in u[0,t>0]==0; should literally match the independent variables”

I dont understand why I’m getting this message.

I found a heat problem similar to this. I tried to use the that equation but changed the boundary problem, but that did not work.

Can you help me detect which part of equation is wrong, and how I can fix the boundary values?




@Martin Wijaya: the boundary conditions as you have given them are not valid Mathematica syntax. The are also obviously not consistent with the initial conditions at u[0,0] and u[2,0]. See Nassers answer for the correct syntax…
– Albert Retey
Dec 18 ’12 at 19:32


2 Answers


You can get rid of the warnings if you modify your initial condition close to the space boundaries :

limit = 3;
newIC[z_] = 100 (HeavisideTheta[z – 10^-limit] – HeavisideTheta[z – (2 – 10^-limit)])

sol[z_, t_] = NDSolve[{Derivative[0, 1][u][z, t] == Derivative[2, 0][u][z, t],
u[0, t] == 0, u[2, t] == 0, u[z, 0] == newIC[z]},
u[z, t], {z, 0, 2}, {t, 0, 100}][[1, 1, 2]]

Plot[{newIC[z], sol[z, 0]}, {z, 0, 2}, PlotRange -> All, Exclusions -> None]



Dear b.gatessucks, what does space boundaries means?what does newIC[z_] means?, I’m a bit confused with the new boundary condition.
– Martin Wijaya
Dec 29 ’12 at 3:06



By space boundaries I mean the points z=0 and z=2. The new initial condition I propose takes the value 100 everywhere except in the immediate vicinity of those two points and so there are no warnings.
– b.gatessucks
Dec 29 ’12 at 9:42



Oh, I see..ok thx.
– Martin Wijaya
Dec 30 ’12 at 13:12

I need to look more into this, the warning about

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.

As I’ve seen this many times. Are you sure the IC’s you are using are correct? (ie. consistent) ? Please double check the textbook or from the source you obtained this.

But just wanted to say that you still get a solution and can animate it just fine. A different method solver might be needed like MethodOfLines. But need to look more into this. Here is the solution

ClearAll[u, z, t]
eq = Derivative[2, 0][u][z, t] == Derivative[0, 1][u][z, t];
ic = {u[0, t] == 0, u[2, t] == 0, u[z, 0] == 10};

sol = u /. First@NDSolve[Flatten[{eq, ic}], u, {z, 0, 2}, {t, 0, 10}]

Plot3D[sol[z, t], {z, 0., 2.}, {t, 0., 10}, PlotRange->All, AxesLabel->{z, t, u[z, t]}]



the boundary conditions can be made consistent with the initial conditions by using something like Exp[-a*t] with an a of users choice. From comparing with giving these bc explicitly I suspect that Mathematica automatically chooses these boundary conditions with a==1 (or something very close to that).
– Albert Retey
Dec 18 ’12 at 19:38



back in the old days (V4) I used to get around discontinuities or inconsistencies with an If in the conditions.
– Mike Honeychurch
Dec 18 ’12 at 21:06



Dear Nasser, Why do I have to use u/.First@NDSolve and Flatten instead of directly typing it as Plot3D[NDSolve[{Derivative[0, 1][u][z, t] == Derivative[2, 0][u][z, t], u[0, t] == 0, u[2, t] == 0, u[z, 0] == 100}, u, {z, 0, 2}, {t, 0, 100}], {z, 0., 2.}, {t, 0., 100}]
– Martin Wijaya
Dec 29 ’12 at 3:04