Numerical method to solve PDE system $\mathbf{A}\dot{\mathbf{U}} + \mathbf{B}\mathbf{U}^{'} + F(\mathbf{U}) = 0$

Question

I would like to numerically solve a system of PDEs of the form

$$\mathbf{A}\dot{\mathbf{U}} + \mathbf{B}\mathbf{U}^{'} + F(\mathbf{U}) = 0$$

Where, $\dot{\mathbf{U}}(s, t)$ represents the time derivative, $\mathbf{U}^{'}(s, t)$ represents the space derivative, $\mathbf{U}$ is the vector of unknowns, and $F(\mathbf{U})$ is a nonlinear function of $\mathbf{U}$.

Both $\mathbf{A}$, and $\mathbf{B}$ are constant singular matrices.

So, the given system is a coupled system of PDEs and ODEs.

How to numerically solve such system of PDEs given initial and boundary conditions?. Which numerical method will be suitable?

Answer 1

Have you looked at the Keller Box scheme? It's described in this paper http://dx.doi.org/10.1016/j.amc.2009.07.054. It's relatively straight forward in principle and suitable for PDE's and non-linear systems. Depending on the size of your system though it can be a bit of a hassle to work out.

I'm not sure what your problem is but I've used it for solving one dimensional Stefan problems, where I have time and space derivatives and non-linear terms. It essentially applies the finite difference operators defined in the paper (pg 1614) and then you can linearise the problem using Newton's method.

Answer 2

There are several related posts on this site where a possible and convenient method is presented. A possible approach known as operator splitting consists in writing two PDE systems \begin{aligned} &(S_1):\qquad {\bf A}\dot{\bf U} + {\bf B} {\bf U}' = 0 \\ &(S_2):\qquad {\bf A}\dot{\bf U} + F({\bf U}) = 0 \end{aligned} which are then integrated successively. The splitting error resulting from this process is often small enough compared to the numerical errors corresponding to the integration of $S_1$ and $S_2$, see this post. For the numerical resolution of $S_1$ a finite-difference method could be used, or a finite-volume scheme for instance. For the numerical resolution of $S_2$, a well-suited ODE solver will do. However, for most schemes the matrix $\bf A$ is required non-singular. One may need to introduce a change of variable or a perturbation of the initial system to rewrite it in a non-singular form. The most appropriate approach is very likely to depend on the form of the matrices and functions involved.