How Does StreamPlot Handle Singularity or Indeterminates?

Question

I will like to know how StreamPlot finds its way around singular/indeterminate values. Actually, I was thinking of using ListStreamPlot to generate a phase portrait for a huge dataset that contains some singularities. To illustrate my point here is a simple example that was previously dealt with in this post

Clear[f]
f[x_, y_] := {1, x/(-x - 2 y)}
StreamPlot[f[x, y], {x, -2, 2}, {y, -2, 2}, Frame -> False,
Axes -> True, AspectRatio -> 1/GoldenRatio,StreamStyle -> "PinDart"]

Obviously there was a singularity at y=-x/2. However, StreamPlot was still able to generate a plot. Dataset showed that there are singular values in the range of plot (except that StreamPlot seemed to have used some subtle methods to generate the values plotted while avoiding values associated withy=-x/2) - though the discontinuity can be seen in the plot.

exampleDataSet=Table[f[x, y], {x, -2, 2}, {y, -2, 2}]//N produces

{{{1., -0.333333}, {1., -0.5}, {1., -1.}, {1., ComplexInfinity}, {1., 1.}},
{{1., -0.2}, {1., -0.333333}, {1., -1.}, {1., 1.}, {1., 0.333333}}, 
{{1., 0.}, {1., 0.}, {1., Indeterminate}, {1., 0.}, {1., 0.}}, 
{{1., 0.333333}, {1., 1.}, {1., -1.}, {1., -0.333333}, {1., -0.2}}, 
{{1., 1.}, {1., ComplexInfinity}, {1., -1.}, {1., -0.5}, {1., -0.333333}}}

ListStreamPlot[exampleDataSet] will generate a blank plot.

In my situation, the singularities have special physical meaning in the model and should be captured on the phase plot. Interpolating the data won't work either since there are some singularities involved. ListStreamPlot won't work with exampleDataSet. Removing the singular values will change the problem and the phase portrait completely. I am thinking if I know how StreamPlot works its way around those singular values, I can possibly adapt the knowledge to ListStreamPlot or any other phase plot options in Mathematica.

How do I go about generating a phase plot with a dataset (containing singularities) while capturing those singularities in the plot? Are there alternatives that I could explore? Any help will be greatly appreciated.

Answer 1

One idea is to inverse scale the interpolation of the scaled vector field ${\bf F}(x,y)$, in which the scaling factor $\alpha(x,y)$ is zero at the singularities; or might describe the process as conjugating the interpolation operator $I$ by the scaling $\alpha$: $$ \alpha^{-1} \circ I \circ \alpha $$ However, after computing $\alpha \,{\bf F}$, one needs to replace the singular points by an appropriate vector. Here one has to fudge if the only information is the given vector field data.

The vectors of the form $(1, \infty)$ should be replaced by vector of the form $(0,v_y)$, but the magnitude and sign cannot be determined. So I chose $(0,1)$. In the vector field at hand, the field reverses direction at all these points. One can see from the construction, knowledge of which is presumably unknowable in the actual use-case, the singularity does not consist of isolated points but a continuous submanifold. Ideally the scaling function $\alpha(x,y)$ would vanish on the submanifold, and not just at the singular points in the data. (Ideally, we would know the values of $v_y$, too.)

The vectors of the form $(v_x, \hbox{??})$, where "??" denotes Indeterminate, are impossible to replace with an accurate vector, unless the limit of the vector field can be evaluated at the point. I just treated them like the $(1, \infty)$ points.

In the real use-case, if the vector field comes from measured data, I expect that some amount of judgment and manual adjustment of the replacement values might be needed. In the code below, it is all done blindly and automatically, because it's convenient for a rough illustration. It may also be possible that whatever steps in the data generation that led to ComplexInfinity and Indeterminate might be altered to produce appropriate replacement values.

ClearAll[x, y, sc, nix, vf];
exampleDataSet = Table[{1, x/(-x - 2 y)},
     {x, -2, 2, 1/32}, {y, -2, 2, 1/32}] // N // Quiet;
{xdim, ydim} = Dimensions[exampleDataSet][[;; 2]];
xrange = {-2, 2};
yrange = {-2, 2};
xsc = Evaluate@Rescale[#, {1, xdim}, xrange] &;
ysc = Evaluate@Rescale[#, {1, ydim}, yrange] &;

sing = Most /@ 
   Position[exampleDataSet, _?(Not@*NumericQ), {3}, Heads -> False];
(* nix[p] has a scaling factor of zero at a singular point p *)
nix[{x0_, y0_}] := 2/Pi*ArcTan[(x - xsc[x0])^2 + (y - ysc[y0])^2];
(* sc[x, y] is the scaling function $\alpha(x,y)$ *)
sc = Compile @@ {{x, y},
    nix /@ sing // Apply[Times],
    RuntimeOptions -> "EvaluateSymbolically" -> False};
data2 = Table[
    {{xsc[a], ysc[b]}, exampleDataSet[[a, b]]*sc[xsc[a], ysc[b]]},
    {a, xdim}, {b, ydim}] // Quiet;

vf[x_, y_] = 1/sc[x, y]*Interpolation[
     Flatten[ReplacePart[
       data2,
       PadRight[sing, {Automatic, 3}, 2] -> {0., 1.}
       ], 1],
     InterpolationOrder -> 1
     ][x, y];

StreamPlot[vf[x, y], {x, -2, 2}, {y, -2, 2}]

Undersampling does mean that the flow computed by StreamPlot can leak around the singular points Technically, all discrete approximations to the vector field will have leaks through the singular submanifold, but if discretized on a fine enough scale, StreamPlot will not show it. The following shows the above technique on the OP's exampleDataSet:

Interpreting the flow with such coarse data poses some difficulties, I should think.