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.

  • $\begingroup$ 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. $\endgroup$
    – Community Bot
    Dec 3, 2021 at 15:46
  • 2
    $\begingroup$ Also tell us something about what the equations represent and where they come from. $\endgroup$ Dec 3, 2021 at 18:39
  • $\begingroup$ I've updated the question to include more details. Does this provide sufficient information or do you need the full background? $\endgroup$ Dec 4, 2021 at 18:15

1 Answer 1


On a cursory look it seems you are solving a advection diffusion reaction equations. The first one seems to be conservation of mass while the second one is AD equation for variable phi. You can always discretize using finite differences and use upwind schemes for advection terms,central differences for diffusion and a second order implicit or a semi implicit scheme like crank Nicolson for transient term should work. Since we have cross terms and it looks nonlinear a stability analysis of linearized version can help you derive a CFL like condition. I have written a second order accurate axi symmetric code for some industrial heterogeneous catalytic systems and the above said method worked okay though the source terms were a tad simpler. My understanding of FE is that it has some issues in handling advection terms and hence nearly all commerical fluid codes are erstwhile of FD type and now of Finite volume which uses a conservative form of naiver stokes equations.

  • $\begingroup$ 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. $\endgroup$ Dec 6, 2021 at 14:30
  • $\begingroup$ So I've tried the finite difference method for an easier but related system, and it seems like the issue is that the equations quickly become too complicated for larger mesh sizes. Basically for e.g. 10 mesh points in space the system already needs to solve a 10x10 matrix and the expression become extremely. I used Mathematica as an algebraic solver, but maybe you know some way of dealing with this? $\endgroup$ Dec 6, 2021 at 14:46
  • $\begingroup$ 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). $\endgroup$ Dec 6, 2021 at 15:30
  • $\begingroup$ Hey @PianoEntropy, oh yes when we discretize using n points we typically be inversing a (n-2)*(n-2) sparse square matrix. You can either vectorize or code using a simple for loops for writing down the equations and forming the matrix but I have inversed large matrices of the order of 5000*5000,which i believe is relatively small compared to kind of matrices which a low level langauge handles in a commerical system. I have observed if we don't explicitly calculate the inverse but instead solve the equations it's way faster ~ 5-10 secs for the sizes mentioned (in matlab) $\endgroup$ Dec 6, 2021 at 17:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.