I have a question concerning how to solve a 1D PDE with inconsistent boundary conditions using mathematica. The function U is in x and t, with two derivatives in x and one in t. Here is what the PDE looks like (foam drainage): $$ c_1 U \frac{\delta U}{\delta x} + c_2 (\frac{1}{2\sqrt{U}} (\frac{\partial U}{\partial x})^{2} + \sqrt{U}\frac{\partial^2 U}{\partial x^2}) - \frac{\partial U}{\partial x} = \frac{\partial U}{\partial t} $$

Where C1 and C2 are constant parameters. (Yes, this could be simplified somewhat). We need three boundary conditions. For the IC

$$U[x,t=0] = z\:(a\:constant) $$

For the BCs there is first: $$U[x=0,t] = z\:(the\:same\:constant) $$

and the BC that is causing me to lose sleep:
$$U^{3/2}[x=L,t] = -\frac{c_2}{c_1}\frac{\partial U}{\partial x}[x=L,t] $$ that is, at U[x=L,t] the first derivative in x is non-zero (it is "z^(3/2)", in my code, z = 1), but U[x,t=0] is constant (and so a derivative in x for U[x=L,t=0] should be zero). Thus the inconsistent boundary conditions.

In order to clarify what I am doing, here is the code block I built up from the heat equation step by step. I didn't have any problems until adding the Neumann BC for this problem:

C1 = 0.1; 
C2 = 0.1; 
L = 1; 
tfinal34 = 1;  
GovEqFD = 
  C2 Sqrt[u[x, t]] D[u[x, t], {x, 2}] + 
  C2 (1/(2 Sqrt[ u[x, t]])) (D[u[x, t], x])^2 +  
  C1 (u[x, t] - 1) D[u[x, t], x] - D[u[x, t], t] == 0  
BC34 = 
  {u[0, t] == 1, 
   -C2/C1 (D[u[x, t], x] /. x -> L) == (u[L, t])^(3/2) /. x -> L)]}  
IC34 = u[x, 0] == 1  
EqsFD4 = Join[{GovEqFD}, BC34, {IC34}]

FDsol4 = NDSolve[EqsFD4, u[x, t], {x, 0, L}, {t, 0, tfinal34}]

In order to try to find a solution around the inconsistent boundary conditions I have looked into the Method of Lines, but to be honest, I haven't gone too far into trying to figure out how to code it yet. It's not clear to me from the documentation exactly how complicated it would be to use FDM on the above governing equation.

In attempting a work around, I have thought to use WhenEvent to define the boundary condition @ t= 0 and also for t>0. I didn't have any problem with using WhenEvent in terms of getting warnings.

I changed the list for BC34 to:

BC34 = 
  {u[0, t] == 1, 
   -C2/C1 (D[u[x, t], x] /. {x -> L,t -> 0}) == 0, 
     (t > 0), (-C2/C1 (D[u[x, t], x] /. x -> L) == 
       (u[L, t])^(3/2) /. x -> L)]}

As I said there were no warnings, and I saw that beloved output for the interpolating function. However, in plotting I see the solution is the same unstable plane it shouldn't be for the inconsistent solution. It seems that the BC remains -C2/C1 (D[u[x, t], x] /. x -> L) == 0, for times greater than t==0. Can someone help me understand if I am using WhenEvent incorrectly, or how I might implement it in this case to resolve the inconsistency in the BC @ x=L and IC? Or perhaps suggest other solution strategies?

If there is anyway I can clarify my question, formatting, or the information I am providing, please let me know.

  • $\begingroup$ Possible duplicate of Boundary condition with spatial derivative is ignored by NDSolve $\endgroup$
    – xzczd
    Jan 19, 2018 at 4:25
  • $\begingroup$ To be specific, adding Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}} is enough to solve your problem. (The warning will still be there, but don't worry about it. ) $\endgroup$
    – xzczd
    Jan 19, 2018 at 4:28

1 Answer 1


Do pretty much the same as what you have:

C1 = 1/10;
C2 = 1/10;
L = 1;
tfinal34 = 6;

GovEqFD = 
 C2 Sqrt[u[x, t]] D[u[x, t], {x, 2}] + 
   C2 (1/(2 Sqrt[u[x, t]])) (D[u[x, t], x])^2 + 
   C1 (u[x, t] - 1) D[u[x, t], x] - D[u[x, t], t] == 0;

BC34 = {u[0, t] == 1, (D[u[x, t], x] /. x -> L) == -C1/C2 ((u[L, t])^(3/2) /. x -> L)};

The trick is choosing an IC that closely matches the IC we want, but is not too inconsistent with our BC's. I generally try to tack on a parabola on the IC at the BC that is unmatched. Let the deviation from our original IC to a parabola occur at x = u1. In this case u0 is the constant value of the original IC. With the equation of a parabola with vertex at {u0,u1}, use a as the parameter to match the BC at x = L.

u0 = 1;
u1 = 99/100;

y[x_] = u0 - a (x - u1)^2

Find the value of a that matches the right slope at the boundary.

(y'[x] /. x -> L) == -y[L]^(3/2) // N

a = a /. Solve[%, a][[1]]

Transition from the constant value u0 to the parabola with a UnitStep.

IC34 = u[x, 0] == u0 - (UnitStep[x - u1]) a (x - u1)^2

This is what our adjusted IC looks like at the end.

Plot[IC34[[2]], {x, .98, 1}]

enter image description here

The transition is fairly smooth and the deviation from the desired IC is small.

EqsFD4 = Join[{GovEqFD}, BC34, {IC34}];

solu = FDsol4 = 
  NDSolve[EqsFD4, u[x, t], {x, 0, L}, {t, 0, tfinal34}, 
   PrecisionGoal -> Infinity, MaxStepFraction -> 1/200, 
   MaxSteps -> {50000, Automatic}];

u[x_, t_] = u[x, t] /. %[[1]]

gifs = Table[
   Plot[u[x, t], {x, 0, 1}, PlotRange -> {.5, 1}, PlotLabel -> t "t", 
    AxesLabel -> {"x", "u"}], {t, 0, tfinal34, .1}];

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.