1
$\begingroup$

I am trying to numerically solve the following PDE

pde = 
  With[{d = .25, \[Rho] = .5, M = 10, T = 2, 
    L = 10, \[Alpha] = 0.2},
   {
    D[n[x, t], {t, 1}] ==
     Piecewise[
       Table[{(1 - x^2/(3 d (M T - i T)))* 
          n[x, t]* g / Sqrt[3 d (M T - i T)] , 
         i T < t < i T + \[Alpha] T && -\[Sqrt](3 d*(M T - i T)) < 
           x < \[Sqrt](3 d*(M T - i T)) }, {i, 0, M - 1}]
       ]
      +
      Piecewise[
       Table[{d * D[n[x, t], {x, 2}] - \[Rho] * n[x, t], 
         i T + \[Alpha] T < t < (i + 1) T}, {i, 0, M - 1}]
       ],
    n[x, 0] == InitialCondition,
    (D[n[x, t], x] /. x -> -L) == (D[n[x, t], x] /. x -> L) == 0
    }
   ];   

However, when I try to get the solution through

sln = With[{M = 10, L=10},
   ParametricNDSolveValue[
   pde /. InitialCondition -> 0.5  (*uniform*),
   n, {x, -L, L}, {t, 0, M T}, {g},
   AccuracyGoal -> 20, PrecisionGoal -> 10]
  ];  

and then write, for example, sln[.5][0,5], MMA v13.2 gives an error NDSolve::ibcinc stating the initial and boundary conditions are inconsistent, but it also produces a result. Could someone help me figure out what's happening?

I want to keep g as a parameter to then be able to manipulate the plots by varying it, say

Manipulate[Plot[{sln[g][0, t]}, {t, 0, 10}, PlotRange -> All],{g, 0.1, 10}]

However, in that case MMA throws out the error and then stops the evaluation.

PS: The reason I have added PrecisionGoal is that without it, at longer times the solution seems to have a lot of numerical error. For example, after removing PrecisionGoal and AccuracyGoal, MMA solves the pde with no error, but when I look at Plot[{sln[.1][x, 14]}, {x, -10, 10}, PlotRange -> All], the plot does not make sense to me.

$\endgroup$
5
  • $\begingroup$ PS2: I realized with only PrecisionGoal there is no error, but the final results/plots don't make sense. However I got things to work with AccuracyGoal -> 15, PrecisionGoal -> 10. This seems like an accident, since with a different initial condition n[x,0] = Sin[ Pi x /(2L)]^2 the error persists even with these accuracy and precision goals. $\endgroup$
    – SaMaSo
    May 24, 2023 at 20:01
  • $\begingroup$ @AlexTrounev the Neumann boundary condition is on the derivative of the function, not its value $\endgroup$
    – SaMaSo
    May 25, 2023 at 13:01
  • $\begingroup$ What do you try to solve? It looks like delay equation. $\endgroup$ May 25, 2023 at 15:18
  • $\begingroup$ @AlexTrounev It's a reaction-diffusion system, where at specified time intervals there is particle creation (the first piecewise), and in other intervals there is diffusion and particle annihilation $\endgroup$
    – SaMaSo
    May 26, 2023 at 3:59
  • $\begingroup$ Thank you, please, see my answer. $\endgroup$ May 26, 2023 at 4:12

1 Answer 1

1
$\begingroup$

This problem can be solved with using Mathematica FEM as follows

Needs["NDSolve`FEM`"]

{d = .25, \[Rho] = .5, M = 10, T = 2, L = 10, \[Alpha] = 0.2}; mesh = 
 ToElementMesh[ImplicitRegion[-L <= x <= L, {x}], 
  MaxCellMeasure -> 1/100];
pde = {D[n[x, t], {t, 1}] == (f1[t, x, g] + f3[t]) n[x, t] + 
     f2[t] D[n[x, t], x, x], n[x, 0] == InitialCondition};

f1[t_, x_, g_] := 
 g Piecewise[
   Table[{(1 - x^2/(3 d (M T - i T)))/Sqrt[3 d (M T - i T)], 
     i T < t < i T + \[Alpha] T && -\[Sqrt](3 d*(M T - i T)) < 
       x < \[Sqrt](3 d*(M T - i T))}, {i, 0, M - 1}]]; 
f2[t_] := 
 Piecewise[
  Table[{d, i T + \[Alpha] T < t < (i + 1) T}, {i, 0, M - 1}]]; 
f3[t_] := 
 Piecewise[
  Table[{-\[Rho], i T + \[Alpha] T < t < (i + 1) T}, {i, 0, M - 1}]];
sol = NDSolveValue[pde /. {InitialCondition -> .5, g -> .5}, 
  n, {t, 0, M T}, {x} \[Element] mesh];

Visualization

{Plot3D[sol[x, t], {t, 0, T}, {x, -L, L}, ColorFunction -> Hue, 
  AxesLabel -> Automatic, MeshStyle -> White], 
 Plot3D[sol[x, t], {t, 0, M T}, {x, -L, L}, ColorFunction -> Hue, 
  AxesLabel -> Automatic, MeshStyle -> White]}

Figure 1

Update 1. In a case of initial condition in a form of n[x,0] = Sin[ Pi x /(2L)]^2 we also have numerical solution without errors. It looks like this one Figure 2

$\endgroup$
4
  • $\begingroup$ Do we not need to specify boundary conditions? $\endgroup$
    – SaMaSo
    May 26, 2023 at 12:53
  • $\begingroup$ @SaMaSo In a case of Mathematica FEM when using (D[n[x, t], x] /. x -> -L) == (D[n[x, t], x] /. x -> L) == 0 there is automatic NeumannValue[0, True] applied. So you don't need to take into account these boundary conditions. $\endgroup$ May 26, 2023 at 16:31
  • $\begingroup$ Ok thank you. Just to make sure I understand correctly, with FEM I would only need to specify the BC if they are Dirichlet or mixed, or otherwise Neumann is automatically enforced. Is that right? $\endgroup$
    – SaMaSo
    May 30, 2023 at 13:14
  • 1
    $\begingroup$ Yes you are right, NeumannValue[0, True] is automatically enforced. $\endgroup$ May 30, 2023 at 15:15

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.