4
$\begingroup$

I have a function which takes in a parameter t and outputs a point {x, y}, which together describe a continuous and smooth curve. This function is fairly nasty and built of InterpolatingFunctions, a typical example might look something like:

curve = {(0.001936 Cos[2.90772-InterpolatingFunction[Domain: {{0.,10.}}Output: scalar][t]])
/(1+1.020758 Cos[5.84877-InterpolatingFunction[Domain: {{0.,10.}}Output: scalar][t]]),
-((0.001936 Sin[2.90772-InterpolatingFunction[Domain: {{0.,10.}}Output: scalar]t]])
/(1+1.0207584 Cos[5.8487-InterpolatingFunction[Domain: {{0.,10.}}Output: scalar][t]]))}

[BTW: How do I get Mathematica output to format nicely here?] In case it's relevant, these draw out variations on conic sections and cycloids (e.g., the path of the Moon around the Sun), so I know it's very smooth and could even compute some bounds on its maximum curvature if that ends up being needed.

I also have a region defined by t which is considerably nicer, something like:

region = Region[Disk[{Cos[t], Sin[t]}, 1]]

What I need is the smallest value of t within some interval tMin to tMax (corresponding to the domain of the InterpolatingFunctions) where the curve[t] crosses the boundary of region[t], or to know if it doesn't cross at all. I could also start such that curve[tMin] is inside or outside region[tMin].

I tried NSolve with things like

NSolve[curve \[Element] RegionBoundary[region] && tMin < t < tMax, t]

This works if the curve is defined via normal functions, but NSolve refuses to work with InterpolatingFunctions. I also tried FindRoot with variations on

FindRoot[RegionDistance[region, curve], {t, tMin}]

Which does something with interpolating functions but has some serious problems of its own. Firstly, if there are no solutions, it sometimes spits out a seemingly arbitrary wrong answer (with a warning message), and I don't know how to filter those out. Secondly, it sometimes answers outside of [tMin, tMax], and I don't know if that means there are no solutions within the desired interval or if it just homed in on another. Thirdly, even if it does give a valid answer within the range, I don't know if I can conclude that this is the lowest valid t just because I put tMin as the starting point for the search (I presume not in general, but perhaps if my InterpolatingFunctions aren't too pathological? - see top of post) I am considering writing my solver that samples along curve at fixed intervals to find a coarse-grained solution and using that as a starting point for FindRoot, but that seems slow and only partially solves the problems mentioned.

Is there another solver in Mathematica I should use for this problem? Are there ways of patching the issues I listed with FindRoot?

EDIT: Runnable Example Here is concrete code to show what I mean, where I go out of the way to force curve to be an interpolating function. Not that, in general, I'd want the region to be a function of t, but I'm omitting this here to show it as a graph.

dataX = Table[{p, 1.2p}, {p, 0, 1, 0.01}];
dataY = Table[{p, p(2-p)Sin[10p]}, {p, 0, 1, 0.01}];
testCurve = {Interpolation[dataX, t], Interpolation[dataY, t]};
testRegion = Region[Disk[{0.95, 0}, 0.3]];

Show[
ParametricPlot[testCurve, {t, 0, 1}],
Graphics[{Opacity[0.5],Disk[{0.95, 0},0.3]}],
Graphics[{Red, PointSize->Large,Point[testCurve/.t->0]}],
Graphics[{Purple, PointSize->Large,Point[testCurve/.t->0.93]}]
]

Example Sketch

If tMin=0, I start at the red dot and want to get where the line first enters the shaded area; if tMin=0.93, I start at the purple dot and want when the blue line next leaves the shaded area.

$\endgroup$

3 Answers 3

5
$\begingroup$
dataX = Table[{p, 1.2p}, {p, 0, 1, 0.01}];
dataY = Table[{p, p(2-p)Sin[10p]}, {p, 0, 1, 0.01}];
testCurve = {Interpolation[dataX, t], Interpolation[dataY, t]};
testRegion = Region[Disk[{0.95, 0}, 0.3]];
  • Plot the SignedRegionDistance function can find all the four times.

( plot0 is the signed distance between the testCurve and testRegion.When the point of testCurve is outside the testRegion, the sign >0 and when the point of testCurve is inside the testRegion, the sign <0, so we only need to find the roots of such SignedRegionDistance.)

plot0 = Plot[SignedRegionDistance[testRegion, testCurve], {t, 0, 1}, 
   Mesh -> {{0}}, MeshFunctions -> {#2 &}, MeshStyle -> Red];
p = Cases[plot0[[1]], GraphicsComplex[{pts__}, __] :> pts, -1];
meshindexs = First@Cases[plot0[[1]], Point[i_] :> i, -1];
times = SortBy[p[[meshindexs]], First][[;; , 1]]

enter image description here

{0.604517, 0.65739, 0.915983, 0.964347}.

  • color the points.
colors = Join[{Red}, ColorData[97] /@ Range[Length@times - 1]];
Show[ParametricPlot[testCurve, {t, 0, 1}], 
 Graphics[{Opacity[0.5], Disk[{0.95, 0}, 0.3]}], 
 Graphics[{PointSize -> Large, 
   Thread[{colors, Table[Point@testCurve, {t, times}]}]}]]

enter image description here

  • We set MaxRecursion -> 4 or another iternation MaxRecursion -> 1 etc. to test the animation.
reg[center_] := Disk[center, 0.3];
findTimes[testCurve_, testRegion_] := Module[{plot0, p, meshindexs},
   plot0 = 
    Plot[SignedRegionDistance[testRegion, testCurve], {t, 0, 1}, 
     Mesh -> {{0}}, MeshFunctions -> {#2 &}, MeshStyle -> Red, 
     MaxRecursion -> 4];
   p = Cases[plot0[[1]], GraphicsComplex[{pts__}, __] :> pts, -1];
   meshindexs = First@Cases[plot0[[1]], Point[i_] :> i, -1];
   times = SortBy[p[[meshindexs]], First][[;; , 1]]];
Manipulate[
 Module[{times},
  times = findTimes[testCurve, reg[center]]; 
  Show[ParametricPlot[testCurve, {t, 0, 1}], 
   Graphics[{Opacity[0.5], reg[center]}], 
   Graphics[{Red, AbsolutePointSize[8], 
     Point[testCurve /. t -> times[[1]]], Blue, 
     Point[Table[testCurve, {t, Rest@times}]]}]]], {{center, {1, 0}}, 
  Locator}]

enter image description here

$\endgroup$
3
  • $\begingroup$ I have to say I don't understand the solution. Is plot0 actually part of finding the solutions, or just displaying them? Would that make it very slow? What does the mesh do? How would this work in the general case where testRegion is time dependent (and so can't be plotted)? Or if testCurve is periodic? $\endgroup$ Mar 6 at 22:58
  • 1
    $\begingroup$ @ScienceSnake plot0 is the signed distance between the testCurve and testRegion, when the point of testCurve outside the testRegion, the sign >0 and when the point of testCurve inside the testRegion, the sign <0, so we only need to find the roots of such SignedDistance. $\endgroup$
    – cvgmt
    Mar 6 at 23:03
  • $\begingroup$ Isn't this hugely inefficient? Solving an equation by plotting it and manipulating the underlying points seems like it must involve a lot more function evaluation and computation than solving it without resorting to graphical objects? Or am I misunderstanding how Mathematica works completely? $\endgroup$ Mar 6 at 23:17
2
$\begingroup$

You didn't provide runable curve-function, that's why I define my own (examplary, simple)

curve = Function[u, (1 - u)/2 {-1, -1} + (1 + u)/2  {1, 1}] 
region= Circle[{0, 0}, 1]

NMinimize evaluates the boundary crossing points

min1=NMinimize[{RegionDistance[reg, ip[u]], 0 > u }, u]   
(*{3.59472*10^-11, {u -> -0.707107}}*)
min2=NMinimize[{RegionDistance[reg, ip[u]], 0 < u }, u]   
(*{2.8084*10^-11, {u -> 0.707107}}*)

Show[{Graphics[{region, PointSize[Large], Red, 
Point[ip[u /. min1[[2]]]], Point[ip[u /. min2[[2]]]]}], 
ParametricPlot[ip[u], {u, -1, 1}]}]

enter image description here

addendum

To get all intersections for your testCurve try

data = Table[{{p}, {1.2 p, p (2 - p) Sin[10 p]} }, {p, 0, 1, 0.01}];
testRegion = Circle[{0.95, 0}, 0.3] 
testCurve = Interpolation[data ]
(u /.FindInstance[{RegionDistance[testRegion, testCurve[u]] == 0, 0 < u < 1}, u, 10])
sol =DeleteDuplicates[%, (#1 - #2)^2 < 10^-5 &]
(*{0.604516, 0.657398, 0.915953, 0.964371}*)

Smallest parameter gives the first intersection point

Show[{ParametricPlot[testCurve[t], {t, 0, 1}], 
Graphics[{ testRegion, Red, PointSize -> Large,Point[testCurve[sol[[1]]]], Blue, Point[testCurve[Rest[sol]]]}] }]

enter image description here

$\endgroup$
4
  • $\begingroup$ Thanks! That seems to do the trick. The difference between Mathematica's numerical methods never ceases to confuse me. However, is there any guarantee that NMinimize will spit out the first crossing? In the example I've added, the line crosses the region 4 times at t>0, can I be sure that NMinimize will give the first one? $\endgroup$ Mar 6 at 18:00
  • 1
    $\begingroup$ @ScienceSnake See my modified answer! $\endgroup$ Mar 6 at 21:25
  • $\begingroup$ Thanks, that works great! I guess I need to do a little bit of work to put in a reasonable number of solutions to find so that I'm sure not to miss the first one, but that's fine. My only issue is that, if there are no solutions, FindInstance keeps running for ages and printing tons of error messages. Is there a way to get it quit neatly and just output no solutions? $\endgroup$ Mar 6 at 22:51
  • 1
    $\begingroup$ @ScienceSnake Try TimeConstrained[] evaluation $\endgroup$ Mar 7 at 13:13
1
$\begingroup$

This is not really an answer, but this seemed like the best way of comparing the two actual solutions proposed

Firstly, I'll start with a more complicated test than the toy example in the question to be more thorough:

r=0.5
dataX = Table[{t, -Cos[1.234t]}, {t, 0, 100, 0.01}];
dataY = Table[{t,1.3Sin[1.234t]-0.3}, {t, 0, 100, 0.01}];
testCurve = {Interpolation[dataX, t], Interpolation[dataY, t]};
testRegion = Region[Disk[2{Cos[t],Sin[t]}, r]];
testBoundary = Circle[2{Cos[t],Sin[t]}, r];

Here the curve describes a particle going around in an ellipse, while the test region is moving in a circle outside of it and its only when the two motions coincide (which they do because they have different periods) that the particle crosses into the region. The SignedRegionDistance is messy to handle numerically with lots of local minima and true solutions that are sometimes very close and sometimes very far apart, as seen below.

Signed Distance

Based on the answers above, I came up with 3 candidate solutions:

cvgmt[region_, curve_, T_]:=Module[{plot0,p, sols},
plot0 = Plot[SignedRegionDistance[region, curve], {t, 0, T}, Mesh -> {{0}}, MeshFunctions -> {#2 &}];
p = Cases[plot0\[[1]], GraphicsComplex[{pts__}, __] :> pts, -1];
sols = Cases[plot0\[[1]], Point[i_] :> i, -1];
If[Length[sols]>0, 
SortBy[p\[[First@sols]], First]\[[1,1]],
T]
]


ulrichRegion[region_, curve_, T_]:= Module[
{sols=Check[FindInstance[{SignedRegionDistance[region, curve] ==0, 0<t<T}, t, 8, WorkingPrecision->10], {{t->T}}]},
Sort[t/.sols]\[[1]]
]


ulrichBoundary[boundary_, curve_, T_]:= Module[
{sols=Check[FindInstance[{RegionDistance[boundary, curve] ==0, 0<t<T},t, 8, WorkingPrecision->10], {{t->T}}]},
Sort[t/.sols]\[[1]]
]

Testing them on the example above gives:

{runtime, sol}=Timing[cvgmt[testRegion, testCurve, 100]]
SignedRegionDistance[testRegion, testCurve]/.t->sol

out = {0.501759,23.7918}
out = -0.00188925
{runtime, sol}=Timing[ulrichRegion[testRegion, testCurve, 100]]
SignedRegionDistance[testRegion, testCurve]/.t->sol

out = {1.0071,23.7903}
out = 2.66343*^-13
{runtime, sol}=Timing[ulrichBoundary[testBoundary, testCurve, 100]]
SignedRegionDistance[testRegion, testCurve]/.t->sol

out = {1.37135,23.7903}
out = 7.23437*^-11

So FindInstance is 2 to 3x slower than the plot method, but gives much more accurate solutions, even with a relatively low WorkingPrecision. Interestingly, using SignedBoundaryDistance with the region is notably faster than BoundaryDistance with the boundary of said region. Also, experimenting around shows that the later often requires telling FindInstance to search for a larger number of solutions before it correctly finds the first one.

Next I wanted to see what happened when there are no valid solutions. This is done in this problem by shrinking the testRegion by setting r=0.4. The output, in the same order as before, is now:

Out = {0.476073,100}

TemplateBox[{"FindInstance", "nsmet", "\"The methods available to FindInstance are insufficient to find the requested instances or prove they do not exist.\"", 2, 363, 25, 28874346979540420740, "Local"}, "MessageTemplate"]
Out = {4.8762,100}

TemplateBox[{"FindInstance", "nsmet", "\"The methods available to FindInstance are insufficient to find the requested instances or prove they do not exist.\"", 2, 365, 26, 28874346979540420740, "Local"}, "MessageTemplate"]
Out = {4.82453,100}

In all cases it returns the default value I wanted. But the FindInstance methods still output an error message despite the enclosing Check and are way, way slower. If there's a way to prevent this huge slowdown I'd consider the question solved thanks to insights from both answers.

$\endgroup$
8
  • $\begingroup$ This test is not right. Please note that when we want to test a moving testRegion, the center of testRegion must depend on another variable, for example s, and we need to defind testCurve by two variales t and s. $\endgroup$
    – cvgmt
    Mar 7 at 12:05
  • $\begingroup$ For example, we can test Clear[testRegion]; logspiralMap[t_] := With[{a = .05, b = .15}, a*E^(b*t) {Cos[t], Sin[t]}]; logspiral2d[s_] := ParametricPlot[logspiralMap[t], {t, 0, s}, PlotStyle -> Red, PlotRange -> All]; r = .5; testRegion[s_] := Disk[2 {Cos[s], Sin[s]}, r] ; Manipulate[ Graphics[{logspiral2d[s][[1]], testRegion[s]}, PlotRange -> 8], {s, .1, 40}, SaveDefinitions -> True] $\endgroup$
    – cvgmt
    Mar 7 at 13:37
  • $\begingroup$ No, the test as posted is correct. The curve is the set of points at all t outputted by an interpolating function. The region t is also a function of the same t. I'm interested in the lowest value of t (above some minimum) the output of the function lies within the region. Edited the OP to clarify that. You can think of it as the curve describing the trajectory of a particle, and I want to know when that particle enters a moving region. $\endgroup$ Mar 7 at 14:00
  • 1
    $\begingroup$ OK I see now. This question is only about the intersection of a moving disk and a moving point,not as my though that the intersection of the disk and the trail of curve.Both of the two case can be done by RegionDistance and SignedRegionDistance. $\endgroup$
    – cvgmt
    Mar 7 at 15:01
  • 1
    $\begingroup$ FindRoot to improve cvgmt. cvgmt1[region_, curve_, T_] := Module[{plot0, p, sols, data}, plot0 = Plot[SignedRegionDistance[region, curve], {t, 0, T}, Mesh -> {{0}}, MeshFunctions -> {#2 &}]; p = Cases[plot0[[1]], GraphicsComplex[{pts__}, __] :> pts, -1]; sols = Cases[plot0[[1]], Point[i_] :> i, -1]; data = If[Length[sols] > 0, SortBy[p[[First@sols]], First][[1, 1]], T]; FindRoot[ SignedRegionDistance[region, curve], {t, data}]]; {runtime, sol} = Timing[cvgmt1[testRegion, testCurve, 100]]; SignedRegionDistance[testRegion, testCurve] /. sol $\endgroup$
    – cvgmt
    Mar 8 at 21:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.