2
$\begingroup$

I'm trying to solve a set of PDEs and am running into problems. The system is defined by: System of PDEs

As can be seen, there are two dynamical equations for $\rho$ and $\phi$, but the equation for $v$ has no time derivative on it. I take Dirichlet boundary conditions for all variables, $\rho(x=-L) = \rho(x=L) = \rho_0$, $\phi(x=-L) = 1, \phi(x=L) = 0$ and $V(x=-L)=v(x=L)=0$. Initial conditions are $\rho(x, t=0)= \rho_0$, $\phi(x, t=0)$ is a sigmoid function going from 1 to 0 and $V(x, t=0)=0$.

Somehow, NDSolve is unable to solve this system under certain conditions. In that case, when I run NDSolve it seems to generate an infinite loop and never terminates, while using up a lot of CPU/memory on my PC.

For reference, here is the current working code thanks to @Alex Trounev (sorry for the poor formatting):

kA := 1/τ ((ρdA - ρ[x, t])/ρdA);
kB := 1/τ ((ρdB - ρ[x, t])/ρdB); 
e := (eM (1 - ϕ[x, t]) + eO ϕ[x, t]) ; 
kD := α e; params = {ρdA -> 0.3, ρdB -> 
   1, ρh -> 1, τ -> 0.5, 
  D1 -> 0.01, η -> 0.1, ξ -> 0.1, kF -> 1, z -> 1, 
  e0 -> 1, α -> 0, eM -> 1, eO -> 1, 
  fc -> 0, χ -> 1, κ -> 
   1}; tmax = 1; L = 1; iv = {ρ[x, 
     0] == (Tanh[8 x] + 1)/2 (ρdA - ρdB) + ρdB, ϕ[
     x, 0] == (1 - Tanh[8 x])/2, v[x, 0] == 0} /. 
  params; bc = {ρ[-L, 
     t] == (Tanh[-8 L] + 1)/
       2 (ρdA - ρdB) + ρdB, ρ[L, 
     t] == (Tanh[8 L] + 1)/
       2 (ρdA - ρdB) + ρdB, ϕ[-L, 
     t] == (1 + Tanh[8 L])/2, ϕ[L, t] == (1 - Tanh[8 L])/2, 
   v[-L, t] == 0, v[L, t] == 0} /. params;
PDEsys = {D[ρ[x, t], t] + 
     D[ρ[x, t] (v[x, t]), 
      x] == (kA (1 - ϕ[x, t]) + kB ϕ[x, t]) ρ[x, t], 
   D[ϕ[x, t], t] + (v[x, t]) D[ϕ[x, t], x] == 
    D1 D [ϕ[x, t], {x, 2}] + 
     z D1 D[ρ[x, t], 
       x] D[ϕ[x, t], x]/ρ[x, t] + (kB - kA + kF) ϕ[x, 
       t] (1 - ϕ[x, t]) + 
     kD (1 - ϕ[x, t]), η D[v[x, t], {x, 2}] - ξ v[x, 
       t] == e D[ρ[x, t], x]/ρ[x, t] + 
     D[e, x] Log[ρ[x, t]/ρh]};
vars = {ρ, ϕ, v}
fullsys = Simplify@Join[PDEsys /. params, bc, iv];
ndsol = NDSolve[fullsys, vars, {x, -L, L}, {t, 0, tmax}, 
  Method -> {"EquationSimplification" -> "Residual", 
    "IndexReduction" -> Automatic}]

Now I'm simplifying a few conditions, and the question is how to get it working when eM is not equal to eO, and when alpha is not equal to zero.

$\endgroup$
7
  • $\begingroup$ When you say "unable to solve the system", what exactly do you receive back? A solution that you believe is incorrect? An error, and if so, which error? $\endgroup$
    – MarcoB
    Sep 29, 2021 at 12:14
  • $\begingroup$ When I run NDSolve it seems to generate an infinite loop and never terminates, while using up a lot of CPU/memory on my PC. $\endgroup$ Sep 29, 2021 at 13:19
  • $\begingroup$ It also generates a warning "NDSolve::pdord" with "Some of the functions have zero differential order, so the equations will be solved as a system of differential-algebraic equations" $\endgroup$ Sep 29, 2021 at 13:20
  • $\begingroup$ You did not specify "params" and "bcDiriclet" and "iv " $\endgroup$ Sep 29, 2021 at 14:50
  • 2
    $\begingroup$ @PianoEntropy If you don't provide complete code and definitions, then we cannot run the code. If we cannot run the code, then we can't help you. I am not sure how you imagine we could help otherwise. $\endgroup$
    – MarcoB
    Sep 29, 2021 at 15:39

1 Answer 1

2
$\begingroup$

There are exist range of parameters where NDSolve getting solution, for example,

params = {kA -> .1, kB -> 1, D1 -> .1, \[Eta] -> 1, \[Xi] -> 1, 
  kD -> 1, kF -> .5, 
  z -> 1}; rh0 = 1; tmax = 1; L = 1; iv = {\[Rho][x, 0] == 
   rh0, \[Phi][x, 0] == (1 - Tanh[8 x])/2, 
  v[x, 0] == 0}; bc = {\[Rho][-L, t] == rh0, \[Rho][L, t] == 
   rh0, \[Phi][-L, t] == (1 + Tanh[8 L])/2, \[Phi][L, 
    t] == (1 - Tanh[8 L])/2, v[-L, t] == 0, v[L, t] == 0};

PDEsys = {D[\[Rho][x, t], t] + 
     D[\[Rho][x, t] (v[x, t]), 
      x] == (kA (1 - \[Phi][x, t]) + kB \[Phi][x, t]) \[Rho][x, t], 
   D[\[Phi][x, t], t] + (v[x, t]) D[\[Phi][x, t], x] == 
    D1 D[\[Phi][x, t], {x, 2}] + 
     z D1 D[\[Rho][x, t], 
       x] D[\[Phi][x, t], x]/\[Rho][x, t] + (kB - kA + kF) \[Phi][x, 
       t] (1 - \[Phi][x, t]) + 
     kD, \[Eta] D[v[x, t], {x, 2}] - \[Xi] v[x, t] == 
    D[\[Rho][x, t], x]/\[Rho][x, t]};
vars = {\[Rho], \[Phi], v}
fullsys = Join[PDEsys /. params, bc, iv];

ndsol = NDSolve[fullsys, vars, {x, -L, L}, {t, 0, tmax}, 
  Method -> {"EquationSimplification" -> "Residual", 
    "IndexReduction" -> Automatic}]

Visualization

{Plot3D[\[Rho][x, t] /. ndsol[[1]], {x, -L, L}, {t, 0, tmax}, 
  Mesh -> None, ColorFunction -> "Rainbow", AxesLabel -> Automatic, 
  PlotLabel -> "\[Rho]"], 
 Plot3D[\[Phi][x, t] /. ndsol[[1]], {x, -L, L}, {t, 0, tmax}, 
  Mesh -> None, ColorFunction -> "Rainbow", AxesLabel -> Automatic, 
  PlotLabel -> "\[Phi]"], 
 Plot3D[v[x, t] /. ndsol[[1]], {x, -L, L}, {t, 0, tmax}, Mesh -> None,
   ColorFunction -> "Rainbow", AxesLabel -> Automatic, 
  PlotLabel -> "v"]}

Figure 1

$\endgroup$
10
  • $\begingroup$ Very interesting! It seems that with a smaller choice for tmax and L this gives the right results. It still puzzles me why one needs specific parameters and not everything works. But this already helps a lot, thanks! $\endgroup$ Sep 30, 2021 at 16:34
  • 1
    $\begingroup$ There is region of stability for numerical method. Actually parameters $L=1; tmax=1$ are not small, but about unit. For some numerical methods it is optimal combination. In this example we use DAE solver. $\endgroup$ Oct 1, 2021 at 4:58
  • $\begingroup$ If you have a little more time, here's another question: if I replace the third PDE by \[Eta] D[v[x, t], {x, 2}] - \[Xi] v[x, t] == e D[\[Rho][x, t], x]/\[Rho][x, t], and define e:=(eM (1 - \[Phi][x, t]) + eO \[Phi][x, t]), then it doesn't solve it anymore, and outputs "The DAE solver failed at t = 0. The solver is intended for index 1 \ DAE systems and structural analysis indicates that the DAE is \ structurally singular." Any idea how to solve this? $\endgroup$ Oct 1, 2021 at 14:59
  • $\begingroup$ @PianoEntropy Do you mean that e=(eM (1 - \[Phi][x, t]) + eO \[Phi][x, t])=1 at eM=1;eO=1;? $\endgroup$ Oct 1, 2021 at 15:15
  • 1
    $\begingroup$ @PianoEntropy Use my code with 3rd equation in a form \[Rho][x, t] (\[Eta] D[v[x, t], {x, 2}] - \[Xi] v[x, t]) == (1 + kap \[Phi][x, t]) D[\[Rho][x, t], x] with -1<kap<1. $\endgroup$ Oct 1, 2021 at 16:18

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.