# Numerically solving a system of PDEs using change of variables.

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?

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.
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}