I have a plot which suffers from poor numeric accuracy, any tips on how I can improve the quality of this plot?
loss1 = (
4 s^2 - 3 \[Pi] s BesselI[0, 2 s] + 3 \[Pi] BesselI[1, 2 s] -
3 \[Pi] s BesselI[2, 2 s] + 3 \[Pi] s StruveL[0, 2 s] -
3 \[Pi] StruveL[1, 2 s] + 3 \[Pi] s StruveL[2, 2 s])/(8 s^2);
loss2 = (Gamma[0, (2 s)/5] - Gamma[0, 2 s])/Log[5];
LogPlot[{loss1, loss2}, {s, 1, 20}, PlotLabel -> "loss after s steps",
AxesLabel -> {"s"},
PlotLegends -> {"random quadratic", "harmonic decay"}]
Based on asymptotics, the first graph should eventually overtake the second, but I can't plot it far enough to get a sense of how far this is.
This formula was obtained programmatically by calling Laplace transform (background here), so perhaps I should be doing something different upstream
interval[var_, min_, max_] :=
HeavisideTheta[var - min] - HeavisideTheta[var - max];
{xmin, xmax} = {0, 1};
xvals = Range[xmin, xmax, (xmax - xmin)/10];
icdf[y_] =
InverseCDF[WignerSemicircleDistribution[1], y] //
Refine[#, 0 < y < 1] &;
g[x_] = icdf[(1 - x/2)]; (* reverse sort, map 0.5->1 range to 0..1*)
ymin = Limit[g[x], x -> xmax];
ymax = Limit[g[x], x -> xmin, Direction -> "FromAbove"];
gi[y_] =
First@SolveValues[{g[x] == y, xmin < x < xmax}, x,
Assumptions -> {ymin < y < ymax}];
arg = -D[gi[y], y]*y;
fwd = LaplaceTransform[arg*interval[y, ymin, ymax], y, s] /. s -> 2 s;
loss0 = Asymptotic[fwd, s -> 0];
loss = fwd/loss0 // Simplify;
Print["normalized loss after s steps=", loss];
LogPlot[loss, {s, 1, 20}, PlotLabel -> "loss after s steps",
AxesLabel -> {"s"}]