I'm interested in solving a first-order linear PDE with 2 dependent variables in 3 dimensions by the method of characteristics. Something of this general form:


$$
A \frac{\partial u}{\partial x} + B \frac{\partial u}{\partial y} + C \frac{\partial u}{\partial z} + D \frac{\partial v}{\partial x} + E \frac{\partial v}{\partial y} + F \frac{\partial v}{\partial z} + G = 0
$$


I've found a couple of approaches for extending the method of characteristics to 2 dependent variable in 2 dimension problems. For example:


* <https://web.archive.org/web/20220704095859/https://www2.whoi.edu/staff/lpratt/wp-content/uploads/sites/120/2019/07/PrattAppendixB_13509.pdf>
* <https://kyleniemeyer.github.io/gas-dynamics-notes/compressible-flows/method-characteristics.html>


There is a brief discussion of this in the answer to a question about 2x2 systems here: <https://mathoverflow.net/a/318860>



> 
> Lax calculus has been generalized to n×n systems of the rich class. A system is rich if it can be put in a conservative form and it can also be put in a diagonal form. This is called semi-hamiltonian system by the russian school.
> 
> 
> 


This seems like it might be the right direction, but I'm having trouble tracking down any references.


My questions are:


1. Are there particular conditions under which characteristic curves exist for this type of equation?
2. Are there any good references I could look at to better understand how (if it's possible) to find the characteristic curves?


Thanks!



Method of characteristics with 2 dependent variables in 3 dimensions
The method of characteristics is a bit strange here because the equation is underdetermined, so one can't expect to be able to specify a solution by fixing initial data for $u$ and $v$ along a surface in $xyz$ space. It would be more efficient to put the equation in normal form and use the integration method appropriate to that normal form.


It actually helps to consider a slightly more general form of a linear, first-order PDE for $u$ and $v$ as functions of $x$, $y$, and $z$ as follows:
$$
A \frac{\partial u}{\partial x} + B \frac{\partial u}{\partial y} + C \frac{\partial u}{\partial z} + D \frac{\partial v}{\partial x} + E \frac{\partial v}{\partial y} + F \frac{\partial v}{\partial z} 
 + G + Hu+Iv = 0,\tag0
$$
where $A,B,C,D,E,F,G,H,I$ are given functions of $(x,y,z)$. The key to understanding these equations is to understand the properties of the two vector fields
$$
U = A\,\frac{\partial\phantom{u}}{\partial x} +B\,\frac{\partial\phantom{u}}{\partial y} + C\, \frac{\partial\phantom{u}}{\partial z}
\quad\text{and}\quad
V = D\,\frac{\partial\phantom{u}}{\partial x} +E\,\frac{\partial\phantom{u}}{\partial y} + F\, \frac{\partial\phantom{u}}{\partial z}.
$$


In the generic situation, $U$ and $V$ will be linearly independent in the domain $\mathcal{D}$ in $xyz$-space that is being considered, and, moreover, their Lie bracket $W = [U,V]$ will be linearly independent from $U$ and $V$ in $\mathcal{D}$ as well. In this case, it is not hard to show that, by making a change of variables in $xyz$-space together with a change of variables of the form $(u,v)\mapsto (Pu+Qv+M,Ru+Sv+N)$, where $P,Q,R,S,M,N$ are functions of $xyz$ with $PS-RQ\not=0$, the equation can, locally, be put in the normal form
$$
u\_x - v\_y - x \,v\_z = 0.\tag1
$$
In this case, a general solution is constructed by choosing $v(x,y,z)$ arbitrarily and then integrating $v\_y+x\,v\_z$ with respect to $x$ to find $u$.


In the special case that $U$ and $V$ are linearly independent in the domain $\mathcal{D}$ in $xyz$-space that is being considered, and their Lie bracket $W = [U,V]$ is a linear combination of $U$ and $V$ in $\mathcal{D}$, the local normal form of the equation under coordinate changes as above simpifies to
$$
u\_x - v\_y = 0,\tag2
$$
which has, as, general solution, $(u,v) = (P\_y,P\_x)$, where $P$ is any function of $xyz$.


Finally, if $U$ and $V$ do not simultaneously vanish but they are linearly dependent everywhere in $\mathcal{D}$, then, using the above change of coordinates, one can put the equation in the normal form
$$
u\_x - P\,v = 0\tag3
$$
where $P$ is a function of $xyz$. On the open set where $P$ is nonvanishing, the normal form can be simplified to $u\_x-v=0$, while, if $P$ vanishes identically, the normal form simplifies to $u\_x=0$.


In each case, finding the normal form is a standard ODE problem in 3 variables in the domain $\mathcal{D}$. Any invocation of the method of characteristics would have to work in the case of these three or four normal forms, and its implementation would probably involve solving the ODE that put the equations in one of the above normal forms anyway.


Of course, there are intermediate cases, such as when $U$ and $V$ are linearly dependent along a subset of $\mathcal{D}$ or $U$, $V$, and $W = [U,V]$ are linearly dependent along a subset of $\mathcal{D}$, etc.


