Is this bifurcation or numerical error?

Question

I encountered a definite integral (with parameter $k$) that produces the "chaotic" figures below.

$$F(k) \equiv \int_{\theta_k}^{\theta_{k+1}} f(\theta, J \sin\theta - k)\,\mathrm{d}\theta~, \qquad f(\theta,u) \equiv \left( \frac{2u - L\sin\theta }{1 - L \sin\theta} \right)^2 + \left( \frac{2 - 2u + L \sin\theta }{1 - L \sin\theta} \right)^2\\ \theta_k \equiv \sin^{-1}\bigl( \frac{k}{J} \bigr)~, \qquad |L|<1,\quad L\neq0,\quad J>k+1 \geq1,\quad J \in \mathbb{R},~k \in \mathbb{N}$$

It is a line integral along a sinusoidal path $u = J \sin\theta - k$, going from $u = 0$ to $u = 1$, where the angles are $\theta_k$ to $\theta_{k+1}$.

My Question:

Is this a bonafide bifurcation to chaos? Is it some other mathematical phenomenon? Or is it likely an implementation flaw?

Each of the plots is $F(k)$ versus $k$ for a distinct large $J$. Basically my intention was to use $\frac{k}J$ to approximate 0 to 1.

Roughly speaking, the discrete points start out nicely following a curve and then spreads out for perhaps $J > 10^4$.

There is a parameter $L$ which the distribution of the points depends on.

The actual parameters of these 4 plots are not important cuz this bifurcating-like behavior is very common. Following the figures is the code block of a Manipulate module that one can play with and generate similar plots

With[{Mn = 0.01, Mx = 16, SpL = Spacer@15, SpS = Spacer@3, Fz = FieldSize -> 2.5},
DynamicModule[{Opt, indef, Jp = 5(* Log@J *), Ns = 3(* Log@ sample points *), ShowAll = False, ImgSz = 200, PtSz = 1/100},
indef[t_, J_, L_, k_] := 1/(4 \[Pi]) (3 (1 - L^2)^(3/2))/(1 + 3 (1 - L^2)^(3/2))Divide[ 4/Sqrt[1 - L^2] ( (2 J - L (1 + 2 k)) ((2 J - L) (3 L^2 - 2) - 2 L^3 k) ArcTan[(L - Tan[t/2])/Sqrt[1 - L^2]] + (2 J + L (1 + 2 k)) (2 L + 2 J (3 L^2 - 2) + L^3 (2 k - 1)) ArcTan[(L + Tan[t/2])/Sqrt[1 - L^2]]) + 1/(-1 + L^2 Sin[t]^2) (-4 (2 J - L) L (2 - 3 L^2 + L^4) (1 + 2 k) t - L (-4 J L (4 - 5 L^2 + L^4) + 4 J^2 (8 - 5 L^2 + L^4) + L^2 (-5 L^2 + L^4 + 8 (1 + 2 k (1 + k)))) Cos[t] + L^3 ((-2 J + L) (-1 + L^2) (-4 (1 + 2 k) t Cos[2 t] + (-2 J + L) Cos[3 t]) + 8 J (1 + 2 k) Sin[2 t]))
, 2 L^3 (-1 + L^2)];
Manipulate[ (* lazy foolproofing *) If[L == 0, L = Mn (-1)^RandomInteger[]];
Opt = {PlotStyle -> PointSize -> Dynamic@PtSz, ImageSize -> Dynamic@ImgSz, Joined -> False, Filling -> False};
Row@{ DiscretePlot[{   indef[ArcSin@((k + 1)/J), J, L, k] - indef[ArcSin@(k/J), J, L, k], ArcSin@((k + 1)/J) - ArcSin@(k/J)  }, {k, 0, -1 + Floor@J, Max[1, Floor[J]/10^Ns]}, Evaluate@Opt, PlotRange -> If[ShowAll, All, Automatic]], 
DiscretePlot[ Divide[indef[ArcSin@((k + 1)/J), J, L, k] - indef[ArcSin@(k/J), J, L, k]
   , ArcSin@((k + 1)/J) - ArcSin@(k/J) ], {k, 0, -1 + Floor@J, Max[1, Floor[J]/10^Ns]}, Evaluate@Opt]}
, Column@{ Row@{"Max Log@J", SpS, InputField[Dynamic@Jp, Fz], SpL, 
   "ShowAll", SpS, Checkbox@Dynamic@ShowAll}, Row@{"Log@Ns", SpS, SetterBar[Dynamic@Ns, Range@4]}, Row@{"ImgSz", SpS, InputField[Dynamic@ImgSz, Fz], SpL, "PtSz", SpS, InputField[Dynamic@PtSz, Fz]} }
, {{J, 29000}, Dynamic[10^(Jp - 2)], Dynamic[10^Jp], 
Dynamic[10^(Jp - 2)]}, {{L, -0.17}, -1 + Mn, 1 - Mn, Mn}, TrackedSymbols :> {Jp, J, L, Ns, ShowAll, PtSz}, ControlPlacement -> Left]  ]   ]

The codes above also generates a plot comparing the blue curve $F(k)$ with $\theta_{k+1}-\theta_k$ in yellow.

How $F(k)$ is implemented in the module above: hard-code the indefinite integral $\int f$, which is previously computed outside, then evaluate it at $\theta_{k+1}$ and $\theta_k$ then take the difference.

The code block below is how I compute the indefinite integral and FullSimplify to improve it slightly.

ClearAll[k, J, u, f, indef, indefSimp, L, myAss];
f[t_, u_, L_] :=((2 u - Jc Sin@t)/(1 - Jc Sin@t))^2 + ((2 - 2 u + Jc Sin@t)/(1 + Jc Sin@t))^2;
myAss = -1 < L < 1 < k <= J && K != 0; Timing[indef = Integrate[Sin@t f[t, J Sin@t - k, Jc], t, Assumptions -> myAss]]
Timing[indefSimp = FullSimplify[indef, myAss && 0 < t < Pi/2]]  (* this takes about 30 seconds up to 2 minutes *)

Just in case if this is relevant, I also have a Plot3D module showing the surface of $f(\theta,u)$ like these figures below (left: $L = -0.86$ and right: $L = 0.23$ arbitrary).

Note that $f(\theta,u)$ is a legitimate probability density over $\theta \in [0,2\pi)$ and $u \in [0,1)$ with a normalization constant that is actually implemented in the DiscretePlot Manipulate above (though irrelevant I think there's a multiplicative factor of 4 missing).

With[{ mM = Max[ #[[1]], Min[ #[[2]], Re@#[[3]]   ]   ] &, SpL = Spacer@15, SpS = Spacer@3 (* Spacer Large & Small *), Mn = 0.01, Mx = 10, Fz = FieldSize -> 2.5}, 
DynamicModule[{ImgSz = 250, Opt, intP = True}, Manipulate[{L, p} = 
mM /@ {{-1 + Mn, 1 - Mn, L}, {If[intP, 1, 0], Mx, If[intP, p, Ceiling@p]}};  Opt = {BoxRatios -> {Pi, 2, 2}, PlotPoints -> {4, 1} 10, 
 ClippingStyle -> None, ImageSize -> Dynamic@ImgSz};
Plot3D[ ((2 u - L Sin@t)/(1 - L Sin@t))^p + ((2 - 2 u + L Sin@t)/(1 + L Sin@t))^p , {t, 0, 2 Pi}, {u, 0, 1}, Evaluate@Opt], 
Row@{"ImgSz", SpS, InputField[Dynamic@ImgSz, Fz], SpL, "Integer p", SpS, Checkbox@Dynamic@intP}, {{L, 0.23}, -1 + Mn, 1 - Mn, Mn}, {{p, 2}, Dynamic@If[intP, 1, 0], Mx, Dynamic@If[intP, 1, Mn]}, TrackedSymbols :> {L, p, intP}, ControlPlacement -> Left]  ]  ]

This module of Plot3D is an overkill that covers the $f(\theta,u)$ in question, as a special case where the power is $p=2$, as well as general $p$.

Optional follow up Non-Mathematica${}^{\text{TM}}$ mathematical question

If this is a real mathematical behavior, then it seems that for example $\lim_{J \to \infty} F(k)/(\theta_{k+1}-\theta_k)$ cannot be taken. What would be the proper analysis to bound $F(k)$ for large $J$?

I actually have seen something like this before, but honestly I'm not sure if mine is an example of the ubiquity of chaos or if it's an artifact. Originally I was dealing with something that should converge (the blue dots should be a smooth curve). I cannot make the call if there's something overlooked in my analytic work or is there some unintentional effects due to how I code. For example, $F(k)$ should be non-negative but in the plots there's about half below zero. Any input will be appreciated.

Please let me know if this post can be improved (too few info, too much info, etc). Thank you for your time.

Answer 1

Some time after posting the question, I realized that there is a problem in the mathematical formulation. This is not to say that there isn't a numerical issue, but practically with the following obvious change of variable the computation (and graphics) is stable.

$$F(k) = \int_{u=0}^1 f\left( \sin^{-1}\bigl( \frac{u+k}J \bigr),\; u \right)\,\mathrm{d}u $$

Basically, the original parametrization is similar to calculating an arc length as $\int_{x=0}^{t} (1 + y'^2)^{1/2}\; \mathrm{d}x$ while taking a problematic limit of $\{ y' \to \infty,\; t \to 0 \}$ for infinite local steepness and vanishing transversed distance.

Although the sinusoidal path (of the line integral) can indeed be parametrized however one desires, the angle $\theta$ is a bad choice (among a class?) that does not work.

I'm not deleting the post here on StackExchange because there's still the numerical issue of implementing the "problematic" parametrization. The spreading out of the points is perhaps not an example of bifurcation (as usually seen in iterated mapping), but I'm still not sure that it should have happened even if the wrong parametrization of the integrand is not convergent.