# Stability issue of time dependent PDE with variable dependent coefficients [closed]

I have the following PDE:

$$\qquad \frac{\partial P(x_1,x_2,t)}{\partial t} = -\frac{\partial}{\partial x_1}[F_1(x_1,x_2)P] - \frac{\partial}{\partial x_2}[F_2(x_1,x_2)P] + D(\frac{\partial^2 P}{\partial x_1^2} + \frac{\partial^2 P}{\partial x_2^2} )$$

where

$$\qquad F_1 = \frac{\epsilon^2 + x_1^2}{(1+x_1^2)(1+x_2)}-ax_1$$ and $$F_2 = \frac{1}{\tau_0}(b-\frac{x_2}{1+cx_1^2})$$.

The values of parameters $$\epsilon, a, b, c, \tau_0, D$$ are 0.1, 0.1, 0.1, 100, 5.0, 0.001 respectively.

Updated: I also include the codes after reading through various manuals in Mathematica. I try to implement FEM to solve the aforementioned PDE.

Updated2: Thanks to Andrew, I managed to get the code run as follows:
1. Mesh generation

    Needs["NDSolveFEM"]
mesh = ToElementMesh[Rectangle[{0, 0}, {10, 10}], "MaxCellMeasure" -> 0.1,
"MeshElementType" -> TriangleElement];
mesh["Wireframe"]

1. Boundary and initial condition

Γ =
{DirichletCondition[u[t, x, y] == 0, x == 0 || y == 0 || x == 10 || y == 10]};
ic =
u[0, x, y] ==
PDF[MultinormalDistribution[{5, 5}, {{1/100, 0}, {0, 1/100}}], {x, y}];

2. PDE formula

F1[x_, y_] =
(ϵ^2  + x^2)/((1 + x^2)(1 + y)) - a x  /. ϵ -> 0.1  /. a -> 0.1
F2[x_, y_] =
1/τ*1/(b - y/(1 + c*x^2))  /. τ -> 5  /. c -> 100 /. b -> 0.1
β = {F1[x_, y_],  F2[x_, y_]}
D = 0.001;
eq =
D[u[t, x, y], t] -
d Laplacian[u[t, x, y], {x, y}] +
D[F1[x, y] u[t, x, y], x] + D[F2[x, y] u[t, x, y], y] == 0;

3. Solve PDE

uif = NDSolveValue[{eq, Γ, ic},
u, {t, 0, 100}, {x, 0, 10}, {y, 0, 10}]
Plot3D[-Log[uif[100, x, y]], {x, 0, 10}, {y, 0, 10}, PlotRange -> All]


However, stability issue arises because probability becomes negative at some regions.

May I know how to resolve the issue? I've tried to modify mesh elements since it is more preferable compared to adding artificial diffusion terms.

However, I wanted to try artificial diffusion as well (by following Documentation) but I didn't know how to obtain the norm of smallest elements of convection for my equation.

For initial condition, I tried to make it look like:

(* Initial condition is u[0,1,1] = 1 and others places are 0 *)


so that total probability in the domain of interest is one. Other distribution may be ok.

My purpose is that I want to :

1. solve PDE to observe time evolution to reach steady state solution
2. take the log of the solution and plot it

Can you help me with purpose (1) first ? Then (2) later ? The solution is expected to be something like

## closed as off-topic by Szabolcs, Henrik Schumacher, m_goldberg, bbgodfrey, LCarvalhoJan 3 at 19:09

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Szabolcs, Henrik Schumacher, m_goldberg, bbgodfrey, LCarvalho
If this question can be reworded to fit the rules in the help center, please edit the question.

• Please show what you have tried. I assume you already searched the documentation for PDEs and found NDSolve. If you have specific issues after having looked at the documentation, please point them out. – Szabolcs Dec 26 '18 at 8:13
• Type "partial differential equations" into the Documentation Center's search box. This will reveal a wealth of information on solving PDEs with Mathematica. – m_goldberg Dec 27 '18 at 15:47
• @Szabolcs, m_goldberg: thanks for your comments, I've updated my attempts by including the codes. Hope you can help me figure out because I'm a newbie in this. – canhochoi Dec 29 '18 at 5:24
• The steady state means that the solution $u$ does not depends on the time variable $t$. And for the resulting elliptic problem the solution is zero. – Andrew Dec 29 '18 at 6:43
• That is a trivial solution, and it's more important to observe the time evolution of probability distribution to steady state solution. – canhochoi Dec 30 '18 at 4:10

Lots of syntax errors. Also the parameter b was not defined. It's always possible to look up the documentation to see the format in which arguments for a given function should be. Or copy an example and change as required. To get something out of it you can try the "find ten differences" game for the two codes.

bmesh = ToBoundaryMesh[Rectangle[{0, 0}, {10, 10}]];
bmesh["Wireframe"];
mesh = ToElementMesh[bmesh];
mesh["Wireframe"];
\[CapitalGamma] =
DirichletCondition[u[t, x, y] == 0,
x == 0 || y == 0 || x == 10 || y == 10];
F1[x_, y_] = (\[Epsilon]^2  + x^2)/((1 + x^2)*(1 + y)) -
a*x  /. \[Epsilon] -> 0.1  /. a -> 0.1;
F2[x_, y_] =
1/\[Tau]*1/(b - y/(1 + c*x^2))  /. \[Tau] -> 5  /. c -> 100;
\[Beta] = {F1[x, y],  F2[x, y]};
b = 100;
d = 0.001;
eq = \!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]$$u[t, x, y]$$\) -
d*Laplacian[u[t, x, y], {x, y}] + \[Beta].Grad[
u[t, x, y], {x, y}] == 0;

uif = NDSolveValue[{eq, \[CapitalGamma], u[0, x, y] == 1},
u, {t, 0, 100}, {x, 0, 10}, {y, 0, 10}]

Plot3D[Log[uif[x, y, 1]], {x, 0, 10}, {y, 0, 10}]

• Thanks for your help! However, your equation formula is not similar to the PDE. It misses u[t,x,y]*(grad[[Beta[[1]],{x}]+grad[[Beta][[2]],{y}]). However, when I include the missing term, it pops up "Inconsistent equation dim error". Do you know how to fix it? – canhochoi Jan 2 at 3:18
• D[F1[x,y] u[x,y],x]+D[F2[x,y] u[x,y],y] – Andrew Jan 2 at 7:15
• Thanks @Andrew, now the equation can be solved. However, stability issue arises because the probability becomes negative in some regions. I've tried to changed mesh sizes mesh = ToElementMesh[Rectangle[{0, 0}, {10, 10}], "MaxBoundaryCellMeasure" -> 0.005, "MeshElementType" -> TriangleElement]; but seems like it didn't work. – canhochoi Jan 5 at 14:41
• May be the problem is that solution tends to zero fast. As an option instead of a BVP for a parabolic equation $Lu=0$ one can consider BVP for the equation $Lv=\lambda v$. The solutions are closely related: $u=v e^{-\lambda t}$. Choosing $\lambda>0$ large enough one can get $v$ growing (exponentially) instead of vanishing. Of course in this case another instability could arise. – Andrew Jan 6 at 7:57
• Optimally to take $λ$ equal to the first eigenvalue $λ_1$ of the corresponding elliptic problem. Actually the solution for large t will be like $\varphi_1(x)e^{-\lambda_1t}$, where $\varphi_1(x)$ is the first eigenfunction. – Andrew Jan 6 at 12:45