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:


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




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 = 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
MultilineFunction->None]\)[\[Xi], \[Eta], 
     t] + \[Xi] (\[Xi] u[\[Xi], \[Eta], t] - 2 
RowBox[{"1", ",", "0", ",", "1"}], ")"}],
MultilineFunction->None]\)[\[Xi], \[Eta], t] + 4 
RowBox[{"2", ",", "0", ",", "0"}], ")"}],
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"}]}


| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.