NDSolve for Coupled Second Order PDEs

Question

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.

Answer 1

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"}]}