NDSolve unable to solve system of PDEs

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

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.

• 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? Sep 29, 2021 at 12:14
• 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. Sep 29, 2021 at 13:19
• 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" Sep 29, 2021 at 13:20
• You did not specify "params" and "bcDiriclet" and "iv " Sep 29, 2021 at 14:50
• @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. Sep 29, 2021 at 15:39

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"]}


• 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! Sep 30, 2021 at 16:34
• 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. Oct 1, 2021 at 4:58
• 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? Oct 1, 2021 at 14:59
• @PianoEntropy Do you mean that e=(eM (1 - \[Phi][x, t]) + eO \[Phi][x, t])=1 at eM=1;eO=1;? Oct 1, 2021 at 15:15
• @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. Oct 1, 2021 at 16:18