4
$\begingroup$

(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 enter image description here with enter image description here 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 enter image description here

Attempt 1

I approximate the PDE on the square $x,y \in [0,c]^2$ where the boundary conditions become enter image description here 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)

Plot of u(x,t) Plot of u(x,x)

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 enter image description here

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 🙃

Plot of u(x,y) Plot of u(x,x)

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.

$\endgroup$
1
  • $\begingroup$ For some reasons, I was not able to submit the question with equations so I had to replace them with screenshots... $\endgroup$
    – pp.ch.te
    Oct 30, 2023 at 13:34

1 Answer 1

1
$\begingroup$

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

enter image description here

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!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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