I have a system of 2 non-linear, coupled PDEs that I would like to transform to a stiff ODE system to solve them using the method of lines. The 2 equations:

$$\begin{align} \frac{\partial \phi}{\partial t} &= - \phi^m \frac{P}{R - H(P)(R-1)}\\ \phi^m \frac{P}{R - H(P)(R-1)} &= \nabla \cdot [\phi^n(\nabla P + u_z)] \end{align} $$

$\phi$ and $P$ are my two unknowns, $R$ is a constant, $H(P)$ is the heaviside function of $P$ and $u_z$ is an upward directed unit vector.

The problem is that my second equation is an elliptic equation, and hence I have no time derivative in there. My idea is to transform them to Differential Algebraic Equations (DAE) using finite differences for the space derivatives and then to do index reduction to obtain an ODE system.

I have no experience with DAE, so I am stuck there. If we work in 1D in the z direction, I obtain with finite differences:

$$ \frac{\partial \phi_i}{\partial t} = - \phi_i^m \frac{P_i}{R - H(P_i)(R-1)} \\ \phi_i^m \frac{P_i}{R - H(P_i)(R-1)} = \frac{\bigg[\phi_{i+\frac{1}{2}}^n (\frac{Pc_{i+1}-Pc_{i}}{\Delta z} + 1) - \phi_{i-\frac{1}{2}}^n (\frac{Pc_{i}-Pc_{i-1}}{\Delta z} + 1)\bigg]}{\Delta z} $$

From what I read, it seems similar to a Hessenberg index-1:

$$\begin{align} y' &= f(t, y, z) \\ 0 &= g(t, y, z) \end{align} $$

This would only require 1 additional derivative.

May I have a bit of help on how to start from there? Pointers to do that using symbolic computation like Sympy could also work!

Update 1

On @Wolfgang Bangerth's advice, I've rewritten my second equation by isolating P, my unknown for this equation:

$$ \frac{\partial \phi_i}{\partial t} = - \phi_i^m \frac{P_i}{R - H(P_i)(R-1)} \\ (\frac{\phi_i^m \Delta z}{R - H(P_i)(R-1)} + \frac{\phi^n_{i+\frac{1}{2}}}{\Delta z} + \frac{\phi^n_{i-\frac{1}{2}}}{\Delta z}) P_i - \frac{\phi^n_{i-\frac{1}{2}}}{\Delta z} P_{i-1} - \frac{\phi^n_{i+\frac{1}{2}}}{\Delta z} P_{i+1}= (\phi_{i+\frac{1}{2}}^n - \phi_{i-\frac{1}{2}}^n) $$

  • $\begingroup$ This system doesn't quite seem right. The first equation describes the evolution of $\phi$, and it contains $P$ is one of the quantities. One would think that the second equation therefore tells us how we can compute $P$ from $\phi$, but at the least you don't write it that way. Are you thinking of the second equation as a nonlinear equation for $P$? $\endgroup$ Jan 18, 2022 at 2:48
  • $\begingroup$ Sorry if it is confusing. I wrote it that way because it makes sense physically, but less mathematically. Yeah the second equation is used to compute P. What I did before to solve the system was to use two nested Picard Iteration: A first one solving P implicitly and updating H(P) at each iteration, and a second one for calculating $\phi$ from the new values of P but it is quite inefficient because I wrote it myself. My goal would be to solve my system using DifferentialEquations.jl in julia. But I need an ODE system for that. $\endgroup$
    – Iddingsite
    Jan 18, 2022 at 9:23
  • $\begingroup$ Would it be more readable for you if I isolate P on 1 side? I actually don't know what to do with H(P) as it depends on P.. $\endgroup$
    – Iddingsite
    Jan 18, 2022 at 9:26
  • 1
    $\begingroup$ @lddingsite did you try discretizing the PDE into a DAE and using ModelingToolkit.jl structural simplify for the index reduction and tearing? It may need dummy derivatives to continue imposing the implicit constraint. Navier-Stokes is a classic example where custom solvers are not necessary due to the DAE transformations. $\endgroup$ Jan 20, 2022 at 0:44
  • 1
    $\begingroup$ MTK needs work scaling (and there's a lot of work ongoing in that front), but 200x200x2 PDE should be fine. Otherwise you could try doing some of the reduction by hand and find the pattern, though I find that's a bit tedious and error prone. $\endgroup$ Jan 20, 2022 at 15:11

1 Answer 1


The usual approach to solving these kinds of system is not to actually think of them as an ODE with an algebraic constraint, but to use specialized algorithms that makes use of the structure. The prototypical example from which you can learn is the Navier-Stokes equations, in which you have one equation that contains the time derivative of the velocity, and another than has the divergence constraint. The algorithms that solve these equations often alternate between solving the two equations, or they use one as a predictor of another, or use many other ideas. In most cases, one will not get better than second-order-in-time for the accuracy.

Your equation is not so different. The simplest approach to solving is to alternate: Use the previous $\phi$ to solve for $P$, then use $P$ in the equation for $\phi$ to advance by one time step. Then repeat. This is akin to a first-order (Lie) operator splitting, and in the porous media community is called the IMPES (implicit pressure, explicit saturation) method if you solve the equation for $\phi$ explicitly. You can make this a bit more accurate if you extrapolate from the previous values of $\phi$ to the next time step, and then solve for $P$ with that, instead of taking the old $\phi$. You can then take the previous and the new $P$ to solve the equation for $\phi$ with something like the Crank-Nicolson method.

  • $\begingroup$ Hi, thx for your help. I am aware that it would be the classical way to solve it. I did it for these two equations and it works. I was just looking for another method and the MOL seemed quite attractive because I will not build the solver myself if I can use an ODE solver, and my current code is quite slow for what I want to do. $\endgroup$
    – Iddingsite
    Jan 19, 2022 at 21:44
  • $\begingroup$ But if you tell me that it is a bad idea, I will believe you. $\endgroup$
    – Iddingsite
    Jan 19, 2022 at 23:15

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.