# How to ensure for a solution of NDSolve to be positive?

I'am trying to solve kinetic equation for test particles in the Maxwellian plasma. It is Vlasov equation appended with Landau integral of collisions. It is proven that solution of such an equation is everywhere non-negative. However, NDSolve gives a solution which has small negative values in some small regions where distribution function is expected to be small positive.

MWE below demonstrates the problem. Equation for NDSolve contains an integral of collisions in the form of simplified Rosenbluth's potentials for a case modelling injection of fast ions into a target Maxwellian plasma under the assumption that the density of fast ions is small as compared to the target plasma density. In this case, the stationary kinetic equation reduces to PDE of the 1st order over velocity v and the second order over the angle variable $$\xi=\cos\theta$$. However, similar problem occurs also for the full equation which is second order over v as well.

HeavisidePar[x_] /; Abs[x] <= 1 := 3/4 (1 - x^2)
HeavisidePar[x_] /; Abs[x] > 1 := 0

ClearAll[eq6LC, f, fvLC, \[Tau]s, M, \[Beta]m, \[Xi]0, v0, vc, \
\[Xi]0]
With[{v0 = 10, \[CapitalDelta]\[Xi] = 0.04, M = 8.,
vc = 1, \[Tau]c = 1./100, \[Tau]s = 1,
S0 = 1, \[Beta]m = 1/2, \[Xi]0 = Cos[45. Degree],
pltRange = {Automatic, Automatic, {0, 1}}
},
lbl = StringForm[
"M=, \!$$\*SubscriptBox[\(v$$, $$0$$]\)=, \
\!$$\*SubscriptBox[\(\[Tau]$$, $$c$$]\)=", M, v0, \[Tau]c];
eq6LC =
D[(v^3 + vc^3) f[v, \[Xi]], v] + \[Beta]m/\[Tau]s vc^3/
v^3 \[Tau]s v^2 D[(1 - \[Xi]^2) D[f[v, \[Xi]], \[Xi]], \[Xi]] ==
0;
fvLC = NDSolveValue[{
eq6LC,
f[v0, \[Xi]] == (
S0*\[Tau]s)/(v0^3 +
vc^3) (HeavisidePar[(\[Xi] - \[Xi]0)/\[CapitalDelta]\[Xi]]/\
\[CapitalDelta]\[Xi])
, Derivative[0, 1][f][v, 0] == 0
, f[v, Sqrt[1 - 1/M]] == 0
}, f, {v, 0.0001, v0}, {\[Xi], 0, Sqrt[1 - 1/M]}
];
fsLC[v_, \[Xi]_] /; (v < v0 && \[Xi] <= Sqrt[1 - 1/M]) =
fvLC[v, \[Xi]];
fsLC[v_, \[Xi]_] /; (v >= v0 || \[Xi] > Sqrt[1 - 1/M]) = 0;

plt3DLC =
Plot3D[fsLC[Sqrt[vx^2 + vy^2], vx/Sqrt[vx^2 + vy^2]], {vx, 0,
v0}, {vy, 0.0001, v0},
AxesLabel -> {"\!$$\*SubscriptBox[\(v$$, $$z$$]\)",
"\!$$\*SubscriptBox[\(v$$, $$\[UpTee]$$]\)", f},
PlotLabel -> lbl
, BaseStyle -> {12}, PlotRange -> pltRange]

]


There are some posts in this forum asking similar question for ODE. One of the advises given there is to smooth out sharp corners, e.g. use something smoother than UnitStep. In the MWE, is possible to smooth custom function HeavisidePar that imitates source of fast ions, but it does not help.

Although the problem of negative values of the distribution function seems not critical for me, I would appreciate any advice how to solve it.

Attached picture demonstrates presence of negative regions in the distribution function. They are shaded in gray.

• Somewhat related: mathematica.stackexchange.com/q/10055/1871 Jan 19, 2023 at 8:21
• Since you know that the actual solution is everywhere positive and that the problem is just low precision misrepresenting small values, Plot the Abs of your function. Or just Clip the function. Jan 19, 2023 at 17:26

It is proven that solution of such an equation is everywhere non-negative.

I assume the above is theoretical result. i.e. in a perfect world using exact solver.

But NDSolve is not an exact solver. It is numerical. So there will always be very small approximation in the results. One can tweek the hundreds of options for NDSolve to try to improve the numerical result accuracy.

I tried few options, and found at the end that setting AccuracyGoal->12 for example seems to eliminate the negative output, at least in the plot. But I am sure one can come up with other combinations of setting to do that or improve the accuracy more. I just changed this one line

fvLC = NDSolveValue[{eq6LC,
f[v0, ξ] == (S0*τs)/(v0^3 +
vc^3) (HeavisidePar[(ξ - ξ0)/Δξ]/Δξ), Derivative[0, 1][f][v, 0] == 0,
f[v, Sqrt[1 - 1/M]] == 0},
f, {v, 0.0001, v0}, {ξ, 0, Sqrt[1 - 1/M]}, AccuracyGoal -> 12];


Notice it gives warning

• Yes, positivity of solution is proved analitically. Jan 19, 2023 at 8:42
• Unfortunately, this method is insufficient for the full collision integral. Jan 19, 2023 at 9:18
• @IgorKotelnikov I hope someone can find better way to improve the accuracy. As I mentioned there are many ways to do this and it just a matter of trying different settings and options to find what works best for each case. Jan 19, 2023 at 9:29
• @IgorKotelnikov If $a$ is the AccuracyGoal, then round-off errors less than $10^{-a}$ are treated as acceptable. If a solution satisfies $|u| < 10^{-a}$, then round-off error will probably make it change signs. Making AccuracyGoal greater than the $-\log_10$ of a lower bound on the solution $u$ theoretically would fix the problem. However, there may be other problems as there seem to be here. I don't have a suggestion for solving the other problems, sorry. Jan 21, 2023 at 18:02

Insufficient numerical accuracy leads to small negative function values. For the sake of the plot, a cheap hack is to add a small positive constants to get ride of the negative patches:

... Plot3D[fsLC[Sqrt[vx^2 + vy^2], vx/Sqrt[vx^2 + vy^2]] + 10^-4, ....


• Problem not in getting nice picture. The solution should be used in further calculation, and I want to be sure that it is suffuciently accurate. Jan 22, 2023 at 7:16