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:

dudt=cd2udt2dudt=cd2udt2\frac{du}{dt}=c\frac{d^2u}{dt^2}

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:

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…

– Albert Retey

Dec 18 ’12 at 19:32

=================

2 Answers

2

=================

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.

– 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