3
$\begingroup$

I am trying to solve the PDE below:

xInit = 1;
xMax = 10;
wid = 0.01;
T = 1000000;
Cc = 0.5

tMax = 1;
PDE = {D[c[x, t], t] == D[c[x, t], {x, 2}] + f[x, t], 
   D[f[x, t], t] == -f[x, t]/T, 
   f[x, 0] == (-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, 
   c[x, 0] == 0, (D[c[x, t], x] /. {x -> 0}) == 0, 
   WhenEvent[c[x, t] - Cc == 0, f[x, t] -> 1]};

When I do so with:

PDESol = NDSolve[PDE, {c[x, t], f[x, t]}, {t, 0, tMax}, {x, 0, xMax}, 
   MaxStepFraction -> 0.001];

I get the following error:

NDSolve::nbnum1: The function value -0.5+InterpolatingFunction[{{0.,10.}},{5,7,0,{1001},{6},0,0,0,0,Automatic,{},{},False},{{0.,0.01,0.02,<<45>>,0.48,0.49,<<951>>}},{Developer`PackedArrayForm,{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,<<952>>},{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,<<951>>}},{Automatic}][x]==0 is not True or False when the arguments are <<1>>. >>

I understand that this has something to do with the method of lines and/or how Mathematica is solving this PDE, but I don't understand what is going wrong and how to fix my specific case.

Thanks!

$\endgroup$
3
  • $\begingroup$ What version are you using? I also get an NDSolve::bcart error (insufficient boundary conditions) as well. $\endgroup$
    – SPPearce
    Nov 28, 2018 at 21:06
  • $\begingroup$ Take a look at mathematica.stackexchange.com/questions/133528/… $\endgroup$
    – SPPearce
    Nov 28, 2018 at 21:08
  • $\begingroup$ I'm using 10.0.2.0. The boundary condition error is not a problem; it's a common error in all similar problems I have tried and doesn't cause an issue with the solution so long as the boundary is placed sufficiently far. The difference between my question and the post above is that my replacement rule has two arguments in the function while the posited solution in that post has only 1... I don't see an obvious connection between the two. $\endgroup$ Nov 28, 2018 at 21:31

2 Answers 2

4
$\begingroup$

Too long for a comment, but also an answer to how to use WhenEvent in the "MethodOfLines". WhenEvent may not be the tool to use here.

WhenEvent is a time event. It occurs at a time $t=t_0$ for all $x$. If you want an event to happen depending on the spatial location $x$ and time $t$, such as implied by the condition c[x, t] - Cc == 0, then you cannot do this with WhenEvent.

Here's how one might use WhenEvent in the OP's PDE. I doubt it is the solution the OP wants, but the point of presenting it is to show how to use WhenEvent in a PDE.

As pointed out in the linked answer, the instances of c[x, t] in the WhenEvent code are replaced by an InterpolatingFunction with the following properties:

  1. The "Grid" of the InterpolatingFunction is the spatial grid being used in NDSolve.
  2. The "ValuesOnGrid" of the InterpolatingFunction are the values of c[x, t0] at t0, where t0 equals the time of the event, and for x ranging over the spatial grid.

In an action such as f[x, t] -> val, the expression val needs to evaluate to a value for f[x, t0] for x ranging over the spatial grid. Consequently in the case at hand, val needs to be vector of the same length as the spatial grid.

One further complication is that WhenEvent, while having the attribute HoldAll, evaluates some of its arguments in preparation for integrating the differential equation. So it seemed to me that I had to create functions that would not evaluate until the InterpolatingFunction for c[x, t] was substituted. One can see in the plots below that this works. Again, it might not be what the OP wants; however, it might show the OP whether or not WhenEvent will address her or his needs.

maxval[if : InterpolatingFunction[___][x_]] := 
  Max[if /. x -> "ValuesOnGrid"];  (* could be Min[] or some other function *)
fnew[if : InterpolatingFunction[___][x_]] := 
  ConstantArray[1, Length[if /. x -> "Grid"]];
PDE = {D[c[x, t], t] == D[c[x, t], {x, 2}] + f[x, t], 
   D[f[x, t], t] == -f[x, t]/T, 
   f[x, 0] == (-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, 
   c[x, 0] == 0, (D[c[x, t], x] /. {x -> 0}) == 0, 
   WhenEvent[maxval[c[x, t]] == Cc, f[x, t] -> fnew[c[x, t]]]};

sol = NDSolve[PDE, {c[x, t], f[x, t]}, {t, 0, tMax}, {x, 0, xMax}, 
  MaxStepFraction -> 0.001]

NDSolve::bcart: Warning: an insufficient number of boundary conditions...

Plot3D[c[x, t] /. sol, {x, 0, 10}, {t, 0, 1}, PlotRange -> All, 
 AxesLabel -> Automatic]

enter image description here

We can see the event happened around t = 0.61:

Plot3D[f[x, t] /. sol, {x, 0, 10}, {t, 0, 1}, PlotRange -> All, 
 AxesLabel -> Automatic]

enter image description here

IIRC, WhenEvent in PDEs is not extensively documented. The example in the docs for WhenEvent does not depend on the spatial variable(s). For this reason, I thought this explanation would helpful for the site. (Maybe it deserves its own question?)

$\endgroup$
1
  • $\begingroup$ Hi Michael, Thank you for the lucid and useful response. This information is indeed missing from the documentation. $\endgroup$ Dec 1, 2018 at 18:55
3
$\begingroup$

The PDF can be simplified by integrating the equation for $f (x,t)$ independently of $c (x,t)$, and then completing the definition of $f (x,t)$ on the line $c(x,t)=Cc$. The result is shown below.

xInit = 1;
xMax = 10;
wid = 0.01;
T = 1000000;
Cc = 0.5;

tMax = 1;

PDE = {D[c[x, t], t] == 
    D[c[x, t], {x, 2}] + 
     If[c[x, t] == Cc, 1, 
      Exp[-t/T]*(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi], c[x, 0] == 0, 
   Derivative[1, 0][c][0, t] == 0, c[xMax, t] == 0};
sol = NDSolveValue[PDE, c, {t, 0, tMax}, {x, 0, xMax}, 
   MaxStepFraction -> 0.001];

f[x_, t_] := 
 If[sol[x, t] == Cc, 1, 
  Exp[-t/T]*(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi]
{ContourPlot[sol[x, t], {x, 0, 4}, {t, 0, 1}, Contours -> 20, 
  PlotRange -> All, ColorFunction -> Hue, PlotLegends -> Automatic, 
  PlotLabel -> "c(t,x)", FrameLabel -> {"x", "t"}], 
 ContourPlot[f[x, t], {x, 0, 2}, {t, 0, 1}, Contours -> 20, 
  PlotRange -> All, PlotLegends -> Automatic, 
  ColorFunction -> "Rainbow", FrameLabel -> {"x", "t"}, 
  PlotLabel -> "f(x,t)"], 
 Plot[f[x, 0], {x, 0, 2}, AxesLabel -> {"x", "f(x,0)"}]}

fig1

In the statement of the author can not even see the line c[x,t]-Cc==0. I added code to solve this problem.

xInit = 1;
xMax = 10;
wid = 0.01;
T = 1000000;
Cc = 0.5;

tMax = 1;
PDE = {D[c[x, t], t] == D[c[x, t], {x, 2}] + f[x, t], 
   D[f[x, t], t] == -f[x, t]/T, 
   f[x, 0] == (-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, c[x, 0] == 0, 
   Derivative[1, 0][c][0, t] == 0, c[xMax, t] == 0, 
   WhenEvent[c[x, t] - Cc == 0, f[x, t] -> 1]};
sol = NDSolve[PDE, {c[x, t], f[x, t]}, {t, 0, tMax}, {x, 0, xMax}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 80, "MaxPoints" -> 100, 
       "DifferenceOrder" -> "Pseudospectral"}}] // Quiet
{ContourPlot[c[x, t] /. sol, {x, 0, 4}, {t, 0, tMax}, Contours -> 20, 
  PlotRange -> All, ColorFunction -> Hue, PlotLegends -> Automatic, 
  PlotLabel -> "c(t,x)", FrameLabel -> {"x", "t"}], 
 ContourPlot[f[x, t] /. sol, {x, 0, 2}, {t, 0, tMax}, Contours -> 20, 
  PlotRange -> All, PlotLegends -> Automatic, 
  ColorFunction -> "Rainbow", FrameLabel -> {"x", "t"}, 
  PlotLabel -> "f(x,t)"], 
 Plot[{(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, f[x, t] /. sol /. t -> 0,
    f[x, t] /. sol /. t -> tMax}, {x, 0, 2}, AxesLabel -> {"x", "f"}]}

fig2 Generally speaking the optionWhenEvent[c[x, t] - Cc == 0, f[x, t] -> 1] does not give the desired effect, so I wrote another code in which on the line $c (x, t) = Cc$, the function is defined as $t = t0 (x)$.Then the solution is determined. So that we can clearly see the effect I put T = 10.

xInit = 1;
xMax = 10;
wid = 0.01;
T = 10;
Cc = 0.5;
tMax = 1;
PDE1 = {D[c[x, t], t] == 
    D[c[x, t], {x, 2}] + 
     If[t >= t0[x], Exp[-(t - t0[x])/T], 
      Exp[-t/T]*(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi], c[x, 0] == 0, 
   Derivative[1, 0][c][0, t] == 0, c[xMax, t] == 0};
PDE = {D[c[x, t], t] == D[c[x, t], {x, 2}] + f[x, t], 
   D[f[x, t], t] == -f[x, t]/T, 
   f[x, 0] == (-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, c[x, 0] == 0, 
   Derivative[1, 0][c][0, t] == 0, c[xMax, t] == 0};
sol1 = NDSolveValue[PDE, c, {t, 0, tMax}, {x, 0, xMax}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 800, "MaxPoints" -> 1000, 
       "DifferenceOrder" -> "Pseudospectral"}}];

t0 = Interpolation[
  Table[{x, t /. FindRoot[sol1[x, t] == Cc, {t, .5}]}, {x, 0, 2, .1}]]
Sol2 = NDSolveValue[PDE1, c, {t, 0, tMax}, {x, 0, xMax}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 800, "MaxPoints" -> 1000, 
       "DifferenceOrder" -> "Pseudospectral"}}];

f1[x_, t_] := 
 If[t >= t0[x], Exp[-(t - t0[x])/T], 
  Exp[-t/T]*(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi]
{ContourPlot[Sol2[x, t], {x, 0, 4}, {t, 0, tMax}, Contours -> 20, 
  PlotRange -> All, ColorFunction -> Hue, PlotLegends -> Automatic, 
  PlotLabel -> "c(t,x)", FrameLabel -> {"x", "t"}], 
 ContourPlot[f1[x, t], {x, 0, 2}, {t, 0, tMax}, Contours -> 20, 
  PlotRange -> All, PlotLegends -> Automatic, 
  ColorFunction -> "Rainbow", FrameLabel -> {"x", "t"}, 
  PlotLabel -> "f(x,t)"], 
 Plot[{(-ArcTan[(x - xInit)/wid] + Pi/2)/Pi, f1[x, t] /. t -> 0, 
   f1[x, t] /. t -> tMax}, {x, 0, 2}, AxesLabel -> {"x", "f"}]}

fig3

$\endgroup$
1
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – Kuba
    Dec 1, 2018 at 21:33

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.