# NDSolve for Coupled Second Order PDEs

I have the following set of equations that I have reduced down from a larger third-order partial differential equation

$f(x,y,t) = p(x,y,t) - (x-y)\left(\frac{\partial}{\partial x} - \frac{\partial}{\partial y} \right)p(x,y,t)$

$\frac{\partial}{\partial t}f(x,y,t) = -(x-y)^{2}p(x,y,t) - (x-y)\left(\frac{\partial^{2}}{\partial x^{2}} - 2\frac{\partial^{2}}{\partial y \partial x} + \frac{\partial^{2}}{\partial y^{2}} \right)p(x,y,t)$

My only stipulation for initial conditions is that the functions begin the same so I have chosen arbitrarily:

$f(x,y,0)=p(x,y,0)=e^{-(x^{2}+y^{2})}$.

They also need to dissipate at some boundary $L$ so:

$f(L,y,t)=f(-L,y,t)=f(x,L,t)=f(x,-L,t)$

and

$p(L,y,t)=p(-L,y,t)=p(x,L,t)=p(x,-L,t)$.

My attempt at solving this is as follows:

L = 3
T = 3
sol1 = NDSolve[{{
v[x, y, t] == u[x, y, t] - (x - y) (D[u[x, y, t], x] - D[u[x, y, t], y]),
D[v[x, y, t], t] == -(x - y)^2 (u[x, y, t]) - (x - y) (D[u[x, y, t], x, x] - 2*D[u[x, y, t], y, x] + D[u[x, y, t], y, y])},
{u[x, y, 0] == v[x, y, 0] == If[-L < x < L && -L < y < L, Exp[-(x^2 + y^2)], 0],
u[L, y, t] == u[-L, y, t] == u[x, L, t] == u[x, -L, t] == 0,
v[L, y, t] == v[-L, y, t] == v[x, L, t] == v[x, -L, t] == 0}},
{u, v}, {t, 0, T}, {x, -L, L}, {y, -L, L}, Method -> Automatic]


Unfortunately this produces the error

 NDSolve::pdord: Some of the functions have zero differential order, so the equations will be solved as a system of differential-algebraic equations.
NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.


If I remove the latter terms in both PDEs then Mathematica doesn't produce such errors.

I'm really struggling to find a way to get Mathematica to solve this PDE. It really had trouble with a third order pde that dissipated to zero at the boundaries, and I'm still finding problems after reducing it into a form that I had hoped would be easier for Mathematica to handle.

We make a change of variables $\xi =x-y, \eta=x+y$, we set $p=p[\xi,\eta,t]$, then the system is simplified and reduces to a single wave equation that does not depend on $\eta$. We solve this equation with periodic boundary conditions. The figures show the final state and motion of the initial pulse, reflection from the wall, etc.

L = 3; T = 1;
eq =
\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[\[Xi], \[Eta],
t] + \[Xi] (\[Xi] u[\[Xi], \[Eta], t] - 2
\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[\[Xi], \[Eta], t] + 4
\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"2", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[\[Xi], \[Eta], t]) == 0;
ic = u[\[Xi], \[Eta], 0] == Exp[-(\[Xi]^2 + \[Eta]^2)];
bc = {u[-L, \[Eta], t] == u[L, \[Eta], t] == 0,
u[\[Xi], -L, t] == u[\[Xi], L, t]};
sol = NDSolveValue[{eq, ic, bc},
u, {\[Xi], -L, L}, {\[Eta], -L, L}, {t, 0, T}]
{ContourPlot[sol[\[Xi], \[Eta], T], {\[Xi], -L, L}, {\[Eta], -L, L},
Contours -> 20, ColorFunction -> "TemperatureMap",
PlotLegends -> Automatic, PlotRange -> All,
FrameLabel -> {"\[Xi]", "\[Eta]"}],
ContourPlot[sol[\[Xi], 0, t], {\[Xi], -L, L}, {t, 0, T},
Contours -> 20, ColorFunction -> "TemperatureMap",
PlotLegends -> Automatic, PlotRange -> All,
FrameLabel -> {"\[Xi]", "t"}]}