# Best method to solve this system of PDEs?

I have a system of PDEs constituting an initial value problem (IVP) consisting of three coupled PDEs:

\begin{align} \partial_t \rho + \partial_x(\rho v) &= \left(k_A (1-\phi) + k_B \phi \right)\rho\ \\ \partial_t \phi + v \partial_x{\phi} &= D \partial^2_x \phi + D \frac{\partial_x \rho \partial_x \phi}{\rho} + \left(k_B - k_A \right) \phi(1-\phi) + k_D (1-\phi) \\ \eta \nabla^2 v - \gamma v &= E \frac{1}{\rho} \partial_x \rho \end{align} Here,
\begin{align} E &= \left( E_A + (E_B - E_A) \phi \right) \rho \\ k_D (1-\phi) &= \alpha (E_B - E_A) \phi(1 - \phi) \\ k_i(\rho) &= \frac{1}{\tau} \frac{\rho^c_i-\rho}{\rho^c_i} \hspace{2cm} \text{(i=A, B)}, \end{align} and the rest of the variables are constants.

The equations come from continuum mechanics applied to biological tissues, in this case composed of two different cell types A and B. The unknown variables to solve for are $$\rho(x, t), \phi(x,t)$$ and $$v(x,t)$$. $$\rho$$ represents a density field, $$0 \leq \phi \leq 1$$ represents a fraction of cells belonging one of the two cell types and $$v$$ is the velocity of the medium. $$k_A, k_B, k_D$$ are reaction rates, $$D$$ is a diffusion constant, and $$\eta, \gamma, E$$ are material constants. $$\rho_i^{c}$$ represents a homeostatic density and $$\alpha, \tau$$ are proportionality constants.

I take Dirichlet boundary conditions for all variables, $$\rho(x=-L) = \rho(x=L) = \rho_0$$, $$\phi(x=-L) \approx 1, \phi(x=L) \approx 0$$ (such that it is consistent with $$\phi(x, t=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$$.

I've been looking into some high-level methods to solve this IVP, but none of them work well enough. The best so far has been Mathematica (link to mathematica.stackexchange post), which works to some extent. However, the solutions are highly unstable. For some parameters I get a nice solution, but when I change some parameters Mathematica gets stuck in some infinite loop and never finishes the computation.

Now I'm learning a finite element solver (https://fenicsproject.org/), but I'm not even sure whether in principle this system can be solved using finite elements. Could someone confirm whether this is the right approach?

Or should I be looking at a more low-level approach and program the numerics myself?

My guess is that some methods for solving the Navier-Stokes equations or those for elasticity might also be applicable here, but I'm totally inexperienced with numerically solving PDEs so correct me if I'm wrong.

• Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
– Community Bot
Dec 3, 2021 at 15:46
• Also tell us something about what the equations represent and where they come from. Dec 3, 2021 at 18:39
• I've updated the question to include more details. Does this provide sufficient information or do you need the full background? Dec 4, 2021 at 18:15

• Thanks! The equations are indeed analogous to advection diffusion reaction equations, but coupled to equations for linear elasticity. The third equation is force balance equation arising from some assumptions and is time-independent (i.e. there is no $\partial_t v$ anywhere). I hope this does not affect the rest of your answer? Yes, finite differences seem like a possible way forward. Dec 6, 2021 at 14:30
• The problem with not having a time derivative for $v$ is that one can only use implicit methods to solve for $v$, but these have the problem that they scale poorly with mesh size (quickly becomes too cumbersome). Dec 6, 2021 at 15:30