numerical approach for system of non-linear partial-ordinary differential equations

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.

• Why are g, $\xi$, and $\psi$ not functions of X since they depend on f? Jun 21, 2016 at 18:08
• 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 Jun 30, 2016 at 16:55
• OK, thanks for the clarification. So, if I understand right, you have to solve the ODE at each spatial location? Jun 30, 2016 at 17:18
• Yep thats right. The system of ODE's needs to be solved at each spatial-step Jun 30, 2016 at 17:31

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.

• 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 Jun 21, 2016 at 17:20
• Yes, there is no need to use WENO for other equations, since they are ODEs. 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?

• 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. Jun 30, 2016 at 18:01
• 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) Jun 30, 2016 at 18:37
• 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. Jun 30, 2016 at 18:44
• @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? Jun 30, 2016 at 19:19
• 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)? Jun 30, 2016 at 19:57