Numerically solving a system of PDEs using change of variables.

Question

I am trying to numerically solve the system of PDE's written in matrix form as: $$ \begin{bmatrix} (\lambda + 2\mu)\partial_x & \lambda\partial_y \\ \mu\partial_y & \mu\partial_x \end{bmatrix}\begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} $$

using a finite difference discretization and multigrid Gauss-Seidel relaxation to solve the resulting linear system. In a course handout I read that this could be done efficiently by decoupling the equations using a change of variables using the transpose of the cofactor matrix of the system above:

$$ \begin{pmatrix} u \\ v \end{pmatrix} = \begin{bmatrix} \mu\partial_x & -\lambda\partial_y \\ -\mu\partial_y & (\lambda+2\mu)\partial_x \end{bmatrix} \begin{pmatrix} t_1 \\ t_2 \end{pmatrix} $$

so that the system decouples as:

$$ \begin{bmatrix} (\lambda\mu + 2\mu^2)\partial_x^2 - \lambda\mu\partial_y^2 & 0 \\ 0 & (\lambda\mu + 2\mu^2)\partial_x^2 - \lambda\mu\partial_y^2 \end{bmatrix} \begin{pmatrix} t_1 \\ t_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} $$

Now I get that I could discretize this, solve for $t_1$ and $t_2$ and then compute $u$ and $v$. However in the handout it is stated that a full change of variables is not actually necessary, and that "the solution process can be arranged such" that I am effectively solving the decoupled system, but in terms of the original variables $u$ and $v$.

Unfortunately the explanation ends there. How would this work in practice?

Answer 1

Rewrite the system as $AU_x + BU_y =0$, where $U =(u,v)^T$ and $$ A=\begin{bmatrix} \lambda +2\mu &0\\ 0& \mu \end{bmatrix} ,\qquad B=\begin{bmatrix} 0 &\lambda\\ \mu& 0 \end{bmatrix} . $$ For non-zero parameters, $A$ is invertible and the system reads $U_x + M U_y = 0$ with $M=A^{-1}B$. This system can be solved numerically and analytically by using diagonalization, and then transforming back to $U$. Setting $V=P^{-1}U$ and $M = P\Lambda P^{-1}$ with $$ P=\begin{bmatrix} -c & c\\ 1 & 1 \end{bmatrix} ,\qquad P^{-1}=\tfrac12\begin{bmatrix} -1/c & 1\\ 1/c & 1 \end{bmatrix}, \qquad c = \sqrt{\tfrac{\lambda}{\lambda+2\mu}}\; , $$ the system is rewritten as $V_x + \Lambda V_y = 0$ where $\Lambda = \text{diag}\lbrace{-c}, c\rbrace$ is diagonal.

A standard numerical procedure consists in using finite-volume schemes, such as the Upwind scheme, the Lax-Friedrichs method or the Lax-Wendroff method.

The analytical procedure is based on the method of characteristics (several related posts on this site), which provides the coordinates $(F(y+cx), G(y-cx))^T$ of $V$ where $F$, $G$ are arbitrary functions deduced from the boundary conditions. The coordinates of $U = PV$ are therefore \begin{aligned} u(x,y) &= c \big( G(y-cx) - F(y+cx)\big) \\ v(x,y) &= G(y-cx) + F(y+cx) \end{aligned}