Piecewise BVP not working on interval for NDSolve

Question

My Mathematica skills aren't the best but this solution isn't making sense considering the boundary value I am supplying.

The leftmost boundary condition, i.e. at $x=0$, $c[t,x]$ should be $1$ for $t\leq 1$ and $0$ for $t>1$. i.e. $$c[t\leq1,0]=1$$ $$c[t>1,0]=0$$ However the solution from NDSolve when evaluated at any $t>1$, say, $c[2,0]$ is non-zero. Why? Is it the way I am solving the PDEs or? I have tried other methods other than the piecewise but also no luck.

Clear[Derivative]
tfinal = 50;
L = 10;
difeqns1 = {D[c[t, x], {t, 1}] + D[c[t, x], {x, 1}] + 
     D[q[t, x], {t, 1}] == 0, 
   D[q[t, x], {t, 1}] == c[t, x] - q[t, x], 
   c[0, x] == If[x != 0, 0, 1], 
   c[t, 0] == Piecewise[{{1, t <= 1}, {0, t > 1}}, 0], 
   q[0, x] == If[x != 0, 0, 0]};
solv = NDSolve[
   difeqns1, {c[t, x], q[t, x]}, {t, 0, tfinal}, {x, 0, L}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> "TensorProductGrid"}];
c[t, x] /. solv /. t -> 2 /. x -> 0

Why these If[x != 0, 0, 1] & If[x != 0, 0, 0] code constructions? Which output do you expect? — Ulrich Neumann, Jan 18 '24 at 17:04
Hi @UlrichNeumann If[x != 0, 0, 0] Should return 0 everywhere initially for $q$. I constructed it that way because I wanted to have nonzero options as well, but yes for this example I could have simply have said q[0,x]==0. As for If[x != 0, 0, 1] that should return $c[0,0]=1$ and $c[0,x>0]=0$ unless I have misinterpreted how that works. As far as I could tell, $c[0,0]$ should return a value of $1$ in both those conditions? — Kendall, Jan 18 '24 at 18:56
To be more specific, setting Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100}} should resolve the problem. To learn more about DifferentiateBoundaryConditions, see also https://mathematica.stackexchange.com/a/127411/1871 — xzczd, Jan 19 '24 at 08:57
@xzczd Thank you for this, I love the amount of detail in your response. I realised, more specifically after seeing Ulrich's response, that the boundary conditions are not so much the problem but rather my understanding on how to implement NDSolve correctly for a given PDE and BCs/ICs in a correct manner. This is due to my lack of understanding in Numerical techniques/methods. The MethodOfLines gave a better fit than the FiniteElement in the end. Edit: This PDE has an analytic solution that I used for comparison, but the more complex version (That I will attempt myself) does not. — Kendall, Jan 19 '24 at 12:42

Ulrich Neumann · Accepted Answer · 2024-01-19T08:31:05.470

Assuming If[x != 0, 0, 0] equals zero, using "FiniteElement" for the complete system and changing the condition c[0, x] == If[x != 0, 0, 1] to c[0, x] == Exp[-x/eps] to ensure consistent boundary and inital conditions we get

tfinal = 50;
L = 10;
eps = .01;
difeqns1 = {Derivative[1, 0][c][t, x] + Derivative[0, 1 ][c][t, x] +Derivative[1, 0][q] [t, x] == 0,Derivative[1, 0][q][t, x] == c[t, x] - q[t, x],
c[0, x] ==  Exp[-x/eps] , 
c[t, 0] ==  Piecewise[{{1, t <= 1}, {0, t > 1}}, 0] ,
q[0, x] == 0 };
{cc, qq} = 
NDSolveValue[difeqns1, {c , q }, {t, 0, tfinal}, {x, 0, L},Method ->  "FiniteElement" ]
Plot[cc[t, 0], {t, 0, tfinal}, PlotRange -> All]

addendum

With smoothed bc c[t,0]==(1 - Tanh[(t - 1)/eps])/2 MethodOfLines works too

Clear[Derivative]
tfinal = 10;
L = 50;
eps = .01;
difeqns1 = {Derivative[1, 0][c][t, x] + Derivative[0, 1 ][c][t, x] + 
     Derivative[1, 0][q] [t, x] == 0, 
   Derivative[1, 0][q][t, x] == c[t, x] - q[t, x],
   c[0, x] == {Exp[-x/eps], Max[0, ( eps - x)/eps]}[[1]], 
   c[t, 0] == {(1 - Tanh[(t - 1)/eps])/2, 1 - UnitStep[t - 1], 
      Piecewise[{{1, t <= 1}, {0, t > 1}}, 0]}[[   1]],
   q[0, x] == 0 };
 {cc, qq} = 
 NDSolveValue[difeqns1, {c , q }, {t, 0, tfinal}, {x, 0, L}, 
  Method -> {  
     "FiniteElement", {"MethodOfLines", "TemporalVariable" -> t, 
      "SpatialDiscretization" -> "FiniteElement" }}[[-1]],MaxStepSize->eps   ] 
Plot[cc[t, 0], {t, 0, tfinal}, PlotRange -> All]

Thank you Ulrich, my boundary conditions work fine when I change the method to FiniteElement, and I get the same output as your solution. Why does method of lines not work here? I thought method of lines was good for time dependent PDEs. — Kendall, Jan 18 '24 at 19:09
@Kendall With smoothed bc MethodOfLines works too , see my modified answer — Ulrich Neumann, Jan 19 '24 at 08:23
Thank you, it appears it's not really a boundary problem but more so my own understanding of how to restructure the boundary conditions in such a way as to be consistent with a chosen method. I highly appreciate your feedback. — Kendall, Jan 19 '24 at 12:37

Piecewise BVP not working on interval for NDSolve

1 Answers1