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.
Plot
theAbs
of your function. Or justClip
the function. $\endgroup$