Numerical solution of second-order linear hyperbolic PDE

Question

(I'm also searching for analytical solutions to this PDE; check the bountied questions here and here if you have any ideas)

---

I'm trying to find the numerical solution of the following 2D second-order linear hyperbolic PDE with where $k$, $s$, $t$, $h$, $m$, $M(0)$ are given constants.

The PDE is to be solved for all $x,y>0$ and the boundary conditions are

Attempt 1

I approximate the PDE on the square $x,y \in [0,c]^2$ where the boundary conditions become and use the following code for some choice of the parameters

pde = D[u[x,y], x, y] + k (D[u[x,y], x] + D[u[x,y], y]) + (k^2 - s^2 P[x-y]) u[x,y] == f[x] f[y];
bc1 = u[x,0] == M[x] M0;
bc2 = u[0,y] == M0 M[y];
bc3 = u[x,c] == M[x] Minf;
bc4 = u[c,y] == Minf M[y];

P[x_] := Exp[-Abs[x]/t]
f[x_] := h + m M[x]
M[x_] := h/(k-m) (1-Exp[-(k-m)x]) + M0 Exp[-(k-m)x]
Minf = h/(k-m);

c = 10;

Block[{k=2, s=1, t=1, h=1, m=0, M0=0}, NDSolve[{pde, bc1, bc2, bc3, bc4}, u, {x,0,c}, {y,0,c}]]

Mathematica more or less delivers (on the right I'm plotting $u(x,x)$ along the diagonal)

together with the warning

NDSolve: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

It's clear that some numerical instability is taking place at large $x$, $y$. Increasing $c$ does not improve the situation. I would like to get rid of this instability.

Attempt 2

The PDE is equivalently solved on the half-quadrant $x>y$ with the boundary conditions

Again I approximate the PDE on the triangle with vertices $(0,0), (c,0), (c,c)$ and running the following code (it is my understanding that unspecified boundary conditions are assumed to be Neumann $\partial u/\partial n=0$, correct me if I'm wrong)

pde = D[u[x,y], x, y] + k (D[u[x,y], x] + D[u[x,y], y]) + (k^2 - s^2 P[x-y]) u[x,y] == f[x] f[y];
bc1 = u[x,0] == M[x] M0;
bc2 = u[c,y] == Minf M[y];

P[x_] := Exp[-Abs[x]/t]
f[x_] := h + m M[x]
M[x_] := h/(k-m) (1-Exp[-(k-m)x]) + M0 Exp[-(k-m)x]
Minf = h/(k-m);

c = 10;
w = Triangle[{{0,0}, {c,0}, {c,c}}];

Block[{k=2, s=1, t=1, h=1, m=0, M0=0}, NDSolve[{pde, bc1, bc2}, u, {x, y} \[Element] w]]

I get the same warning as before and the situation now is even worse 🙃

Attempt 3

I tried increasing AccuracyGoal, PrecisionGoal, MachinePrecision, but didn't get any improvement. I also followed the instructions given in the warning at FEMDocumentation/ref/message/InitializePDECoefficients/femcscd. For both attempt 1 and 2 I tried to add artificial diffusion and tried refining the mesh. Nothing has worked in resolving the instabilities, quite the contrary, the situation got worse.

My question

How can I get a solution for this PDE without numerical instabilities? I would like to get the numerical solution both to the original PDE and the PDE on the half-quadrant $x>y$, if possible.

Answer 1

Here my actual best approach.

I made several assumptions (c>>1):

• If x==y is line of symmetry it is sufficient to consider only an eighth circle
• There is only one boundary condition bc1 left (->DirichletCondition)
• On the remaining boundary we have no flux (NeumannValue==0)

here the code

pde = D[u[x, y], x, y] +k (D[u[x, y], x] + D[u[x, y], y]) + (k^2 - s^2 P[x -y]) u[x,y] == f[x] f[y] ;

bc1 = DirichletCondition[u[x, y] == M[x] M0, y == 0];
P[x_] := Exp[- Abs[x] /t]
f[x_] := h + m M[x]
M[x_] := h/(k - m) (1 - Exp[-(k - m) x]) + M0 Exp[-(k - m) x]
Minf = h/(k - m);

c = 20;
reg = ImplicitRegion[x >= 0 && y >= 0 && x^2 + y^2 <= c ^2 && x >= y, {x, y}];


Block[{k = 2, s = 1, t = 1, h = 1, m = 0, M0 = 0}, 
U = NDSolveValue[{pde, bc1(*,bc2  ,bc3,bc4 *) }, u,Element[{x, y}, reg], 
Method -> {"FiniteElement" , "MeshOptions" -> { "MeshOrder" -> 2
,"MaxCellMeasure" -> Pi/4 c^2/100    }}  ]
Plot3D[{U[x, y], (f[x] f[y])/(k^2 - s^2 P[x - y]) }, Element[{x,y},reg], PlotRange -> {0, 1/2}, PlotStyle -> {Automatic,Opacity[.25]}, Mesh -> None ]
]

Plot shows good agreement of the solution U[x,y] with the stationary solution ((f[x] f[y])/(k^2 - s^2 P[x - y])) (transparent).

I tried to improve the result by changing MaxCellMeasure but NDSolve result gets unstable. Don't know why.

Adding some diffusion to the pde Lapacian[U[x,y],{x,y}](don't know wether that's the correct diffusion for this type of pde) gives a smoother result.

Hope it helps!