PDE of real-world system, integral boundary condition

Question

I've stripped all the physical-significance for clarity, but I know that u[x,t] will be everywhere positive and continuous.

here are the equations in Mathematica code:

eqs = {D[u[x, t], {x, 2}] == D[u[x, t], t],
u[x, 0] == c1,
(D[u[x, t], x] /. x -> 1) == 0,
u[0, t] == 3 c2 - c2 Integrate[u[x, t], {x, 0, 1}]} /. {c1 -> 1, c2 -> 1}

The dependent variable u represents pressure, x represents distance and t time.

The last item in eqs represents a material-balance on the gas in the system - the integral is the amount of gas distributed in the region of interest - it accounts for gas in the region of interest, stating that

 gas in region-of-interest +
 gas in the isobaric region at the surface (not of interest) 

   =
 amount of gas in the entire system at time t=0 (represented by u[x,0]==c1 and the addition of 3 C2 - I stripped numerous physical symbols out so while the equations appear flaky, it's the shape that is of interest, not the internal consistency of the simplified equations).

When I try to solve this in Mathematica, regardless of whether I use Integrate or NIntegrate I get errors:

NDSolve[eqs, {x, 0, 1}, {t, 0, 1}] 

When using Integrate.... I get:

Equation or list of equations expected instead of .... (integral here)

When using NIntegrate I get:

The integrand u[x,t] has evaluated to non-numerical values 

(this error message is no surprise - just including it for completeness).

I modeled this using years ago for some graduate research, and got good solutions (I built the physical system in a lab and measured transient pressure profiles etc. to validate model results, so it's a real-world problem. Now I want to see if I can solve it again with Mathematica).

The actual BC is this:

P[x,t] is the dependent variable; x & t are independent variables, everything else is a known constant.

No postings I can find here or on Wolfram.com indicate how or if MMA can do this.

I would like to know if anyone has solved a similar problem (boundary condition containing an integral of the solution itself) using MMA, and how it might be done.

Answer 1

Because NDSolve cannot accommodate the x=0 boundary condition, it is necessary to perform this computation by discretizing the PDE in x. The resulting do-it-yourself procedure is discussed in Introduction to Method of Lines.

For illustrative purposes, assume that x is divided into five equal segments.

n = 5; h = 1/n;

with a variable defined at each node, at {0, .2, .4, .6, .8, 1}

U[t_] = Table[u[i][t], {i, 0, n}]
(* {u[0][t], u[1][t], u[2][t], u[3][t], u[4][t], u[5][t]} *)

The PDE then decomposes into n ODEs plus one algebraic equation. (The first and last equations are modified to take account of the boundary conditions at x = 0 and x = 1.)

Thread[D[U[t], t] == Join[{0}, ListCorrelate[{1, -2, 1}/h^2, U[t]], 
    2 {u[n - 1][t] - u[n][t]}/h^2]];
eqns = ReplacePart[%, 1 -> u[0][t] == (3 - (Total[U[t]] - u[0][t]/2 - u[n][t]/2) h) c2]
(* {u[0][t] == c2 (3 + 
      1/5 (-(1/2) u[0][t] - u[1][t] - u[2][t] - u[3][t] - u[4][t] - 1/2 u[5][t])), 
    Derivative[1][u[1]][t] == 25 u[0][t] - 50 u[1][t] + 25 u[2][t], 
    Derivative[1][u[2]][t] == 25 u[1][t] - 50 u[2][t] + 25 u[3][t], 
    Derivative[1][u[3]][t] == 25 u[2][t] - 50 u[3][t] + 25 u[4][t], 
    Derivative[1][u[4]][t] == 25 u[3][t] - 50 u[4][t] + 25 u[5][t], 
    Derivative[1][u[5]][t] == 50 (u[4][t] - u[5][t])} *)

Initial conditions at t = 0 are given by

initc = ReplacePart[Thread[U[0] == Table[c1, {n + 1}]], 1 -> eqns[[1]] /. t -> 0]
(* {u[0][0] == c2 (3 + 
      1/5 (-(1/2) u[0][0] - u[1][0] - u[2][0] - u[3][0] - u[4][0] - 1/2 u[5][0])), 
    u[1][0] == c1, u[2][0] == c1, u[3][0] == c1, u[4][0] == c1, u[5][0] == c1} *)

Note that u[0][0] is given by u[0][t]/.t -> 0 for consistency. NDSolve objects, if u[0][0] is set to anything else.

Finally, the n + 1 equations are solved and plotted, below with n = 250 for accuracy.

lines = NDSolve[{eqns, initc} /. {c1 -> 1, c2 -> 1}, U[t], {t, 0, 1}];
ParametricPlot3D[Evaluate[Table[{i h, t, First[u[i][t] /. lines]}, {i, 0, n}]], 
  {t, 0, 1}, PlotRange -> All, AxesLabel -> {"x", "t", "u"}]

The steady state converges on 1.5 as n becomes large.