Some progress can be made as follows. First, rescale X
by 10^7
to avoid large values of X
(at least at first). Then, because X
and x
appear only as a ratio in the integrands, define the functions,
f1[y_?NumericQ] := NIntegrate[Sqrt[t + (1817/2) t^2]/(1 + Exp[t - y]),
{t, 0, Infinity]}, MaxRecursion -> 100]
f2[y_?NumericQ] := NIntegrate[t Sqrt[t + (1817/2) t^2]/(1 + Exp[t - y]),
{t, 0, Infinity}, MaxRecursion -> 100]
Plots of these two functions are
and the ODE is solved by
xm = Rationalize[.001438, 0];
sol = NDSolveValue[{X'[x] + 4*π*x 10^-7 (3/(4*π) -
(Sqrt[2]/π^2*1817^(3/2)*(f1[X[x]/x] + (1817*(f2[X[x]/x]))))) == 0,
X[10^-8] == 0}, X[x], {x, 10^-8, xm}, WorkingPrecision -> 30];
The lower bound of integration is set to 10^-8
to avoid division by zero at x = 0
. The upper bound is set to 0.001438
, because the solution blows up at about that value of x
, as we see in to following plot.
Plot[sol, {x, 10^-8, xm}, PlotRange -> All, AxesLabel -> {x, X},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]
The explosive growth shown follows from the ODE and is not a numerical artifact, although accuracy deteriorates with this rapid growth.
Approximate expressions for f1 and f2
For t
large relative to 2/1817
, f1
can be approximated by
Integrate[Sqrt[1817/2] t/(1 + Exp[t - y]), {t, 0, Infinity}]
(* -Sqrt[1817/2] PolyLog[2, -E^y] *)
and similarly for f2
Integrate[Sqrt[1817/2] t^2/(1 + Exp[t - y]), {t, 0, Infinity}]
(* -2 Sqrt[1817/2] PolyLog[3, -E^y] *)
Unfortunately, this approximation does not help much in solving the ODE.