# 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:

$\frac{du}{dt}=c\frac{d^2u}{dt^2}$

$u$ is pore water pressure

$z$ is depth

$t$ is time

$c$ is a constant

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

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

with boundary conditions as follows:

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

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

u[0, t] = 0 ($u$ 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:

DSolve::litarg:

"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... Dec 18, 2012 at 19:32

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)])

Remove[sol]
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. Dec 29, 2012 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. Dec 29, 2012 at 9:42
• Oh, I see..ok thx. Dec 30, 2012 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). Dec 18, 2012 at 19:38
• back in the old days (V4) I used to get around discontinuities or inconsistencies with an If in the conditions. Dec 18, 2012 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}] Dec 29, 2012 at 3:04