I am interested in the numerical solution of the following system of non-linear partial-differential algebraic equations, where the independent variables are $X$ and $T$, representing non-dimensional space and time, respectively. The fields are:

$$ f = f(X,T)$$

$$ g = g(T)$$

$$ \xi = \xi(T)$$

$$ \psi = \psi(T)$$

The system of 5 equations can be defined by:

$$ \frac{\partial f}{\partial X} = \frac{\partial F}{\partial T}-AH\left( I \psi + Jf \right) + B \frac{\partial f^2}{\partial T^2} - K q ~~~~~~~~~~(1) $$

$$ \frac{\partial g}{\partial T} = q ~~~~~~~~~~(2) $$

$$ \frac{\partial q}{\partial T} = C \left( f-g \right) + L \left( 2 D q^2 - E \left|q \right| q \right) - G H \left( I \xi + Jq \right) ~~~~~~~~~~(3) $$

$$ \frac{\partial \xi}{\partial T} = I \xi + Jq ~~~~~~~~~~(4) $$

$$ \frac{\partial \psi}{\partial T} = I \psi + Jf ~~~~~~~~~~(5) $$

where $F$ is the Burgers conservation form, defined by:

$$ F = f^2/2 $$

and $A-E,G-L$ are constants and $q=q(T)$ is introduced to reduce the system to first order equations only, excluding equation (1).

I am interested if anyone has suggestions for how to approach numerically solving the above system of equations in a fully coupled manner?

As a note, I have already managed to solve the above system, but using a de-coupled approach, i.e: solve equations (2) - (4) using a standard ODE solver with an initial condition for $f = f(0,T)$ and $g = g(0)$. The newly calculated value of $q$ can then be substituted into equation (1) to find $f$ for the next space-step. This is repeated until the final space-step is reached.

Thanks in advance

  • $\begingroup$ Why are g, $\xi$, and $\psi$ not functions of X since they depend on f? $\endgroup$ Jun 21, 2016 at 18:08
  • $\begingroup$ Yes you are right, they are sort of functions of X. The reason I didn't state this is that these variables only depend on X at the current space point $\endgroup$
    – james506
    Jun 30, 2016 at 16:55
  • $\begingroup$ OK, thanks for the clarification. So, if I understand right, you have to solve the ODE at each spatial location? $\endgroup$ Jun 30, 2016 at 17:18
  • $\begingroup$ Yep thats right. The system of ODE's needs to be solved at each spatial-step $\endgroup$
    – james506
    Jun 30, 2016 at 17:31

2 Answers 2


If there are only continuous solutions in your problem, the finite difference method is enough. However, if there exists shock, the finite difference/finite volume WENO/ENO suit for your purpose.

Generally speaking, you should first discretize the systems in space and make it to a systems of ODEs. Then, use Runge-Kutta methods to evolve in time.

  • $\begingroup$ I actually use the WENO finite difference approach for equation (1), as there is the possibility of shocks forming in the solution. However, I have not considered this for the other equations $\endgroup$
    – james506
    Jun 21, 2016 at 17:20
  • $\begingroup$ Yes, there is no need to use WENO for other equations, since they are ODEs. $\endgroup$
    – Michael
    Jun 21, 2016 at 17:25

So I have managed to find two approaches for solving the coupled system:

1) Solve equations (2) - (4) using a standard ODE solver with an initial condition for $f_0=f(0,T)$ and $g_0=g(0)$. The newly calculated value of $q_0$ can then be substituted into equation (1) to find $f$. Use fixed point iteration with under-relaxation until the residuals for $f$ and $g$ converge and then progress to the next space step

2) Convert the system to a DAE system by moving everything to the RHS in (2-5) and discretising all terms in time using finite differences. This can be solved using a DAE solver, and specifying a mass matrix with zeros along the diagonal for (2-5).

Both methods work, however method (1) takes many iteration to converge, even with under-relaxation and method (2) becomes extremely demanding on memory when there are many time points

Following suggestions in the comments, another method is to discretise in $X$, rather than in $T$ and time-march, leading to the following ODE system (after some re-arrangement):

$$\frac{\partial f}{\partial T} = \gamma$$ $$\frac{\partial}{\partial T}(\frac{f^2}{2} + B \gamma) = \frac{\partial f}{\partial X} + ...$$ $$ ... $$

where $\gamma$ is a new variable to remove the 2nd order derivative in (1). However, I am unsure how to treat the term: $\frac{f^2}{2} + B \gamma$ on the LHS. Does anyone have any suggestions?

  • $\begingroup$ I was thinking of a fairly standard method-of-lines approach where you would discretize all five equations in the spatial dimension and then pass the resulting system of ODE to a standard solver for solution in time. $\endgroup$ Jun 30, 2016 at 18:01
  • $\begingroup$ This is essentially what method (2) is. The equations are discretised in $T$ only, using finite differences. However, my understanding is you can only use an ODE solver which is compatible with DAEs for this system as there is no $\partial / \partial X$ terms in (2-5) $\endgroup$
    – james506
    Jun 30, 2016 at 18:37
  • $\begingroup$ No, my suggestion is to discretize in space (X) so that you end up with a standard system of ODE to solve. The ODE system in T is a standard initial value problem that should be solvable with many black box ODE solvers. $\endgroup$ Jun 30, 2016 at 18:44
  • $\begingroup$ @BillGreene, Would you be able to clarify how you would do this discretising in space ($X$) instead of $T$? If you discretise in $X$ I can't see how you would put the equations in the form $d/dT = rhs(X)$, which as I understand is needed for an ODE solver? $\endgroup$
    – james506
    Jun 30, 2016 at 19:19
  • $\begingroup$ Well, I appreciate that solving your eqn 1, a nonlinear, hyperbolic equation with shocks has its own set of discretization issues. But if discretizing $\partial f/\partial X$ with a simple first-order backward difference (i.e. upwind) scheme, is acceptable, wouldn't that lead to an ODE system? The unknowns in the system of ODE would be $f, g, q, \xi, \psi$ at each of the FD points. You would also need a sixth unknown to convert eqn 1 to first order form. Doesn't that give a set of ODE with 6N equations (N=number of FD points)? $\endgroup$ Jun 30, 2016 at 19:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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