Return to Question

1. The functions $g_1$ and $g_2$ are exponential - this makes it hard to implement with the FEM package (as suggested here ) .

2. This is a set of a differential-algebraic equations (I am using the "IndexReduction" method to solve the system).

3. As explained, the last equation is derived from the third one (taking spatial derivative, which is very much like this discussion ).

1. The functions $g_1$ and $g_2$ are exponential - this makes it hard to implement with the FEM package (as suggested here ) .

2. This is a set of a differential-algebraic equations (I am using the "IndexReduction" method to solve the system).

3. As explained, the last equation is derived from the third one (taking spatial derivative, which is very much like this discussion ).

The way I implement this in Mathematica is, I allow for $f$ to vary also in space, but the solution of the last equation ensures that $f$ is spatially independent. I impose periodic B.C. for all the fields (note that even though in principle I don't need to impose B.C. for $f$, because I let $f\equiv f\left(x,t\right)$ for simulation purposes, I have to put B.C. for this field as well).

I would like to impose periodic boundary conditions for $\rho,\chi$ and $f$, but demand that $v$ would vanish at the boundaries. This makes much more sense for the system I am dealing with, and I also suspect that the $v$ dynamics with periodic B.C. generate some convergence problems.

Naive NDSolve with the mixed boundary conditions yields the following error message:

NDSolve`ProcessEquations::pdnbc: Mixed periodic and Dirichlet/Neumann boundary conditions are not allowed.

    \[Alpha] = 1; \[Beta] = 2; \[Gamma] = 2;
g[\[Zeta]_, \[Xi]_] = 
  E^(-(1/\[Xi])) (E^ (-1 - \[Alpha] \[Zeta]) + 
     E^ (-1 + \[Alpha] \[Zeta])) (1 - \[Beta]/\[Zeta]); 
eqns[S_, v_, f_, \[Xi]_] := {D[S, t] + v D[S, x] == -D[v, x] S,
   D[\[Xi], t] + v D[\[Xi], x] == 
    f/S  g[f/S, \[Xi]] (\[Gamma] - \[Xi]), 
   D[v, x, x] == 
    1/(f - S) (D[v, x] D[S, x] + D[S, x] - D[S g[f/S, \[Xi]], x]), 
   D[f, t] == S (1 + D[v, x] - g[f/S, \[Xi]]) - f D[v, x]}; 
initc[\[Rho]init_, finit_, \[Chi]init_, 
   vinit_] := {\[Rho][0, x] == \[Rho]init, 
   f[0, x] == finit, \[Chi][0, x] == \[Chi]init, v[0, x] == vinit};
boundc = {\[Rho][t, -1] == \[Rho][t, 1], 
   f[t, -1] == f[t, 1], \[Chi][t, -1] == \[Chi][t, 1], 
   v[t, -1] == v[t, 1]};
sol = First@
   NDSolve[Join[eqns[\[Rho][t, x], v[t, x], f[t, x], \[Chi][t, x]], 
     initc[1, \[Beta], 1 - 10^-4 Cos[x], 0], boundc], {\[Rho], v, 
     f, \[Chi]}, {x, -1, 1}, {t, 0, 0.05}, 
    Method -> {"IndexReduction" -> Automatic}];
\[Rho] = \[Rho] /. sol[[1]];
v = v /. sol[[2]];
f = f /. sol[[3]];
\[Chi] = \[Chi] /. sol[[4]];

In Mathematica I allow for $f$ to vary also in space, but the solution of the last equation ensures that $f$ is spatially independent. I impose periodic B.C. for all the fields (note that even though in principle I don't need to impose B.C. for $f$, because I let $f\equiv f\left(x,t\right)$ for simulation purposes, I have to put B.C. for this field as well).

I would like to impose periodic boundary conditions for $\rho,\chi$ and $f$, but demand that $v$ vanish at the boundaries. This makes much more sense for the system I am dealing with, and I also suspect that the $v$ dynamics with periodic B.C. generate some convergence problems.

Naive NDSolve with the mixed boundary conditions yields the following error message:

NDSolve`ProcessEquations::pdnbc: Mixed periodic and Dirichlet/Neumann boundary conditions are not allowed.

α = 1; β = 2; γ = 2;
g[ζ_, ξ_] = E^(-(1/ξ)) (E^ (-1 - α ζ) + E^ (-1 + α ζ)) (1 - β/ζ);  

eqns[S_, v_, f_, ξ_] := {
   D[S, t] + v D[S, x] == -D[v, x] S,
   D[ξ, t] + v D[ξ, x] == f/S  g[f/S, ξ] (γ - ξ), 
   D[v, x, x] == 1/(f - S) (D[v, x] D[S, x] + D[S, x] - D[S g[f/S, ξ], x]), 
   D[f, t] == S (1 + D[v, x] - g[f/S, ξ]) - f D[v, x]
 }; 

initc[ρinit_, finit_, χinit_, vinit_] := {
  ρ[0, x] == ρinit, f[0, x] == finit, χ[0, x] == χinit, v[0, x] == vinit
 }; 

boundc = {
   ρ[t, -1] == ρ[t, 1], f[t, -1] == f[t, 1], 
   χ[t, -1] == χ[t, 1], v[t, -1] == v[t, 1]
 }; 

sol = First@
   NDSolve[
    Join[
     eqns[ρ[t, x], v[t, x], f[t, x], χ[t, x]], 
     initc[1, β, 1 - 10^-4 Cos[x], 0],  
     boundc
    ], 
    {ρ, v, f, χ}, {x, -1, 1}, {t, 0, 0.05}, 
    Method -> {"IndexReduction" -> Automatic}
   ];

ρ = ρ /. sol[[1]];
v = v /. sol[[2]];
f = f /. sol[[3]];
χ = χ /. sol[[4]];