# How to solve 2D PDE with a constraint?

As related to the 2D+1 PDE problem,

if now I would like to solve the stationary solution which satisfies $\partial_t u(t,x,y)=0$.

Or equivalently,

$-y\partial_{x}u+\partial_{y}\left[a y+b sin(x)u+c\partial_{y}u\right]=0$

ωcb = -50;
ωct = 50;
ωb = -5;
ωt = 5;
A = 10;
γ = 0.1;
kT = 0.1;

With[{u = u[θ, ω]},
eq = -D[ω u, θ] - D[-A Sin[θ] u, ω] - γ kT D[u, ω] + γ  D[ω u, ω]==0
];

ufun = NDSolveValue[{eq, u[-π, ω] == u[π, ω], u[θ, ωcb] == u[θ, ωct]},u, {θ, -π, π}, {ω, ωcb, ωct}];

Plot3D[Abs[ufun[θ, ω]], {θ, -π, π}, {ω, ωb, ωt},
PlotRange -> All, AxesLabel -> Automatic, PlotPoints -> 50,
BoxRatios -> {Pi, ωb, 1}]


We see the message:

NDSolveValue::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified for {u}; the result may not be unique.
NDSolveValue::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.


As expected, the solution is not unique.

Since the function $u(x,y)$ represents the probability density function, namely, we have the normalization condition :$\int dx\int dy \text{ }u(x,y)=1$.

Question: Is it possible to include this additional constraint in NDSolve to make it unique?

Note:

For the boundary condition, we have $u(-\pi,y)=u(\pi,y)$. In $y$-direction, it is unbounded, but in the code I use a periodic boundary condition at some cutoff coordinates to mimic the case without boundary.

• Please post your code. – David G. Stork Aug 16 '18 at 19:01
• I am not sure my problem is precise or not. I guess my 2D PDE with one periodic boundary condition will not give an unique solution. But with an additional constraint can make the solution unique. – Bob Lin Aug 16 '18 at 19:59
• Are you sure the constraint can make the solution unique? – xzczd Aug 17 '18 at 5:54
• Physically It should be. Since this equation determines a stationary probability density. But mathematically i am not sure. – Bob Lin Aug 17 '18 at 6:55
• In time-dependent case, this normalization condition has been implicitly specified by the amplitude of initial condition function. – Bob Lin Aug 17 '18 at 6:57

Here we have an equation of parabolic type whose solution is defined in two squares $xy>0$. In two other squares where $xy<0$ the solution diverges under any conditions. To determine the solution, we must specify initial conditions for $x = 0$ and boundary conditions.

{a, b, c, k, L, y0, T} = {1, 10, 1, 1, 5, 10^-3, Pi};
u0[y_] := Exp[-(y^2)/2]/(2*Pi)
eq = k*y*D[u[x, y], x] -
D[a*y + b*Sin[x]*u[x, y] + c*D[u[x, y], y], y];
ic = u[0, y] == u0[y];
bc = {u[x, y0] == u0[y0], u[x, L] == u0[L]};
sol = NDSolveValue[{eq == 0, ic, bc}, u, {x, 0, T}, {y, y0, L}]
Plot3D[Abs[sol[x, y]], {x, 0, Pi}, {y, y0, L}, PlotRange -> All,
ColorFunction -> "LakeColors", Mesh -> None]


Example of a solution with symmetric boundary conditions on the edges

{a, b, c, k, L, y0, T} = {1, 10, 1, 1, 5, 10^-3, 4*Pi};
u0[y_] := Exp[-(y - L/2)^2/2]/(2*Pi)
eq = k*y*D[u[x, y], x] -
D[a*y + b*Sin[x]*u[x, y] + c*D[u[x, y], y], y];
ic = u[0, y] == u0[y];
bc = {u[x, y0] == u0[y0], u[x, L] == u0[L]};
sol = NDSolveValue[{eq == 0, ic, bc}, u, {x, 0, T}, {y, y0, L}]
Plot3D[Abs[sol[x, y]], {x, 0, T}, {y, y0, L}, PlotRange -> All,
ColorFunction -> "LakeColors", Mesh -> None, AxesLabel -> Automatic,
PlotPoints -> 50]


• Could I ask how did you include the constraint and why you took different bc? – Bob Lin Aug 17 '18 at 6:20
• The normalization condition does not change anything in the solution, it simply determines the wave scale. It is possible to put the same conditions on two edges, but these conditions must be consistent with the initial data. I added an example with equal conditions. But periodic boundary conditions can not be fulfilled, since the solution is not extended to the region $y<0$. – Alex Trounev Aug 17 '18 at 6:59