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

Question

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

Answer 1

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.

Answer 2

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?