PDE of real-world system, integral boundary condition

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.

• NDSolve cannot handle integro-differential equations. Learn more at Advanced Numerical Differential Equation Solving in the Wolfram Language. Mar 29 '15 at 5:36
• Also, you have two typos. u[0, t] = 3 c2 - c2 ... should be u[0, t] == 3 c2 - c2 ... and NDSolve[eqs, {x, 0, 1}, {t, 0, 1}] should be NDSolve[eqs, u, {x, 0, 1}, {t, 0, 1}], although fixing them leads to other error messages. Mar 29 '15 at 5:45
• Finally, setting c2 -> 0 instead of 1 allows NDSolve to produce an answer, although this is not the problem you wish to solve. In my experience, you need to discretize your integral boundary condition to make headway. Mar 29 '15 at 5:53
• thanks for comments. The typos are transcription errors; the equation in MMA is fully formed and correct Mar 29 '15 at 17:48

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[t], u[t], u[t], u[t], u[t], u[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[t] == (3 - (Total[U[t]] - u[t]/2 - u[n][t]/2) h) c2]
(* {u[t] == c2 (3 +
1/5 (-(1/2) u[t] - u[t] - u[t] - u[t] - u[t] - 1/2 u[t])),
Derivative[u][t] == 25 u[t] - 50 u[t] + 25 u[t],
Derivative[u][t] == 25 u[t] - 50 u[t] + 25 u[t],
Derivative[u][t] == 25 u[t] - 50 u[t] + 25 u[t],
Derivative[u][t] == 25 u[t] - 50 u[t] + 25 u[t],
Derivative[u][t] == 50 (u[t] - u[t])} *)

Initial conditions at t = 0 are given by

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

Note that u is given by u[t]/.t -> 0 for consistency. NDSolve objects, if u 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.