32
$\begingroup$

There are no issues trying to find the intersection points of two defined curves.

    f[x_] := 2 (x - 1) (x - 1.5) (x - .5) (x + .5) (x + 1) (x + 1.5);
    g[x_] := 0.4 x - 0.4;
    Solve[f[x] == g[x], x]
    (*{{x -> -1.39849}, {x -> -1.20949}, {x -> -0.331841}, {x -> 
      0.426865}, {x -> 1.}, {x -> 1.51295}}*)
    Manipulate[Plot[Evaluate[{f[x], g[x]}], {x, -1.5, 1.52}, 
               Epilog -> {Red, PointSize[Large], 
               Point[{#, f[#]} & /@ (x /. 
               Solve[f[x] == g[x], x])]}], {a, -0.8, .5, .2}]

enter image description here

When I try to perform the same but this time with an interpolation function I run into problems.

points = Table[{x, f[x]}, {x, -1.5, 1.5, .25}];
if = Interpolation[points];

I can obtain the first intersection by using the following.

Neither NSolve nor Solve will work with interpolation function.

With Solve this will happen

Solve[if[x] == g[x], x]

Solve::inex: Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help. >>

With NSolve, you get a similar error.

NSolve[if[x] == g[x], x, Reals]

(NSolve[InterpolatingFunction[{{-1.5,1.5}},<>][x]==-0.4+0.4 x,x,Reals])

FindRoot[if[x] == g[x], {x, -1.5}]
(*{x -> -1.36794}*)

Which makes this question different from marking points of intersections between two curves as the problem is specific when using InterpolationFunction

What is the best way to find all intersection points between the range (-1.5 and 1.5)?

Improve this question
$\endgroup$
1

6 Answers 6

Reset to default
27
$\begingroup$

[I gave a similar response some time ago either in StackOverflow or MSE but now I cannot find it.]

One way is to track the solution to the ODE that runs over the difference if[x]-g[x]. Use WhenEvent to record axis crossings. This will find all zeros that do not have multiplicity (that is, that cross transversally). Should also find any that are of odd multiplicity since they still cross the axis, albeit at slope of zero.

Reap[
  odsoln = NDSolve[{y'[x] == if'[x] - g'[x], 
      y[-1.5] == if[-1.5] - g[-1.5], WhenEvent[y[x] == 0, Sow[x]]}, 
     y[x], {x, -1.5, 1.5}];][[2, 1]]

(* Out[21]= {-1.36793657661, -1.19575651253, -0.330396218898, \
0.425433482282, 0.999999923787} *)

Plot[y[x] /. odsoln[[1]], {x, -1.5, 1.5}]

enter image description here

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ Nice .. Have to wonder why FindRoot doesn't have an option to find multiple roots since this is really all it would need to do.. $\endgroup$
    – george2079
    Feb 14, 2014 at 20:59
15
$\begingroup$

I take it that the essential question is how to solve an equation in which a term involves an InterpolatingFunction. If it is merely to plot the points, then I would use andre's method. Otherwise, I would use an approach like Daniel Lichtblau's with a small modification. The rest of this is essentially an extended comment to Daniel's answer and george2079's comment to it.

NDSolve is particularly well-suited to this task. The reasons are manifold. First, the process of integration runs over the whole interval and is likely to discover all the roots where a sign change occurs. Second, root-finding algorithms are invoked when an event in WhenEvent is detected. Further the initial condition for an event will be close to the root and the root-finding algorithms will converge relatively quickly, so quickly that in the OP's example the time it took NDSolve to find all of the roots was about the same as using FindRoot on each of them. So NDSolve is being used to track the function and invoke Finally, if the function is highly oscillatory, NDSolve will adapt its step size so that zero-crossings are not (likely to be) missed.

You get more accurate results if the equation to be solved is passed to WhenEvent instead of y[x] == 0. It seems an important principle, even if the interpolating function is only approximate. The symbol y[x] represents another layer of approximation; it seems clearer to use the actual equation whose solution is desired. It takes slightly less time if you ask that no solution be returned. In the example below I added an oscillatory component.

With[{eq = if[x] - g[x] + Sin[200 x]/10},
 roots = Reap[
    NDSolve[{y'[x] == D[eq, x], y[0] == (eq /. x -> 0), 
       WhenEvent[eq == 0, Sow[x]]},
      {}, {x, -1.5, 1.5}];][[2, 1]]
 ]
(*
  {-0.304168, -0.310282, -0.329946, -0.349218, -0.355714, -1.16778, 
   -1.17454, -1.19405, -1.21169, -1.22167, -1.34002, -1.34738, -1.3668, 
   -1.38682, -1.391, 0.410548, 0.423902, 0.43804, 0.991196, 1.00719, 1.01613}
*)

Plot[if[x] - g[x] + Sin[200 x]/10, {t, -1.5, 1.5}, PlotPoints -> 250, 
 Mesh -> {roots}, MeshFunctions -> {# &}, 
 MeshStyle -> {PointSize[Medium], Red}]

Mathematica graphics

Using the actual equation instead of y[x] produces more accurate results:

With[{eq = if[x] - g[x] + Sin[200 x]/10},
 roots2 = Reap[
   NDSolve[{y'[x] == D[eq, x], y[0] == (eq /. x-> 0), 
      WhenEvent[y[x] == 0, Sow[x]]},
     {}, {x, -1.5, 1.5}];][[2, 1]]
 ];

Through[{Max, Mean, Min}[Abs[if[#] - g[#] + Sin[200 #]/10 & /@ Sort@roots]]]
Through[{Max, Mean, Min}[Abs[if[#] - g[#] + Sin[200 #]/10 & /@ Sort@roots2]]]
(*
  {2.55121*10^-10, 2.31578*10^-11, 9.71445*10^-17}
  {3.15835*10^-7,  1.6696*10^-7,   7.10046*10^-8}
*)
Improve this answer
$\endgroup$
10
$\begingroup$

Here is a method that works with your example :

First, one plots the two functions with a special mesh :

f[x_] := 2 (x - 1) (x - 1.5) (x - .5) (x + .5) (x + 1) (x + 1.5);
g[x_] := 0.4 x - 0.4;
points = Table[{x, f[x]}, {x, -1.5, 1.5, .25}];
if = Interpolation[points];

graphic = 
 Plot[Evaluate[{if[x], g[x]}], {x, -1.5, 1.52}, Mesh -> {{0.}}, 
  MeshFunctions -> {if[#] - g[#] &}, MeshStyle -> PointSize[Large]]

enter image description here

Then one extracts the points from the graphic. There are duplicata that must be removed :

pointsList = 
 Cases[Normal[graphic], Point[x_] :> x, {1, Infinity}]  //
 Gather[#, Abs[#1[[1]] - #2[[1]]] < 0.0001 &] & //
 First /@ # &

{{0.425434, -0.230426}, {-1.19576, -0.877972}, {1., -0.000093329}, {-0.330396, -0.532158}, {-1.36794, -0.947173}, {1.51695, 0.206782}}

Be carefull with this method:
The documentation of MeshFunctions says that the MeshFunctions should normally be chosen to be continuous monotonic functions. This is not the case here.

Improve this answer
$\endgroup$
2
  • $\begingroup$ About the MeshFunctions docummentation, in the docummentation of Mathematica 13.2 I find a single statement related to this: "The m_i should normally be chosen to be continuous monotonic functions." reference.wolfram.com/language/ref/MeshFunctions.html. In the same page there is an example that does not use a monotonic function. It is unclear to me if there is indeed a limitation on the use of MeshFunctions. The statement above is quite unclear. $\endgroup$
    – Davi
    Aug 16 at 12:24
  • $\begingroup$ This solution worked for me. Sine in my case the y-axis value is known (I am dealing with the intersection of a numerical function with a horizontal line), here is a variation to the last part: pointsList = Cases[Normal[plotTest], Point[x_] :> x, {1, Infinity}][[All,1]] // DeleteDuplicates[ #, Abs[#1 - #2] < 0.001 & ] & $\endgroup$
    – Davi
    Aug 16 at 12:33
7
$\begingroup$

The RootSearch package http://library.wolfram.com/infocenter/MathSource/4482/ by Ted Ersek returns all the roots over a specified range as a list of replacement rules.

A RootsInRange function was presented in "Finding Roots in an Interval" in The Mathematica Journal 7(2), 1998 (http://www.mathematica-journal.com/issue/v7i2/ updated code below):

f[x_] := 2 (x - 1) (x - 1.5) (x - .5) (x + .5) (x + 1) (x + 1.5);
g[x_] := 0.4 x - 0.4;
points = Table[{x, f[x]}, {x, -1.5, 1.5, .25}];
if = Interpolation[points];

RootsInRange[{f1_, f2_}, {t_, tmin_, tmax_}, opts___] := 
  Module[{p, r, pts, x, f = Function[t, Evaluate[f1 - f2]]}, 
    p = Plot[f[t], {t, tmin, tmax}, 
    Evaluate[Sequence @@ FilterRules[{opts}, Options@Plot]]]; 
    pts = Cases[First[p], Line[{x__}] -> x, Infinity]; 
    r = Abs[Subtract @@ Last@PlotRange[p]]; 
    pts = Map[First, 
      Select[Split[pts, 
        Sign[Last[#2]] == -Sign[Last[#1]] && 
          Abs[Last[#2] - Last[#1]] < r/2 &], Length[#1] == 2 &], {2}]; 
    pts = (FindRoot[f[t] == 0, {t, Sequence @@ ##}, 
    Evaluate[
      Sequence @@ FilterRules[{opts}, Options@FindRoot]]] &) /@ pts;
    Print[Plot[{f1, f2}, {t, tmin, tmax}, 
      Epilog -> {Green, AbsolutePointSize[4], 
      Point /@ ({t, f1} /. pts)}, 
    Evaluate[Sequence @@ FilterRules[{opts}, Options@Plot]]]]; pts]

RootsInRange[f1_ == f2_, {t_, tmin_, tmax_}, opts___] := 
  RootsInRange[{f1, f2}, {t, tmin, tmax}, opts]

For example

RootsInRange[if[x] == g[x], {x, -1.5, 1.5}]

RootsInRange[if[x] == g[x] - Sin[200 x]/10, {x, -1.5, 1.5}, PlotPoints -> 20]
Improve this answer
$\endgroup$
3
  • $\begingroup$ When I apply the RootsInRange function on a very oscillated input (seismic signal) the roots are well found for a tiny range af frequencies. Unfortunately for a large spectrum these function (i) does not find all the roots or (ii) puts forward a coordinates problem. Consequently I use the solution using NDSolve, even if th computation time is much more important. $\endgroup$
    – Antoine
    Nov 27, 2014 at 12:22
  • $\begingroup$ Can you show with a concrete example how that is related to the question? Add some code maybe? $\endgroup$
    – Öskå
    Nov 27, 2014 at 13:07
  • $\begingroup$ Did you try increasing the number of PlotPoints or MaxRecursion options when dealing with the seismic signal? RootsInRange is not fully adaptive (you could try Ted Ersek's package instead). However, since you have the data, why not just subtract of the values of the intersecting function at the data ordinates and look for sign changes in successive data? This will be much more efficient that Interpolating a highly oscillatory function ... $\endgroup$
    – TheDoctor
    Nov 29, 2014 at 0:55
4
$\begingroup$

Another variation.. Here I'm using Plot's automatic incrementation to trace out a sufficiently smooth set of points, then feeding zero crossings as initial points into FindRoot..

  x /. FindRoot[if[x] - g[x], {x, #[[1, 1]], #[[2, 1]]}] & /@ 
           Select[ Partition[
                 Sort@Last@Last@Reap[
                   Plot[if[x] - g[x], {x, -1.5, 1.5},
                        EvaluationMonitor :> Sow[{x, if[x] - g[x]}]]],
                              2,1] ,
                             #[[2, 2]] #[[1, 2]] <= 0 & ]

{-1.36794, -1.19576, -0.330396, 0.425433, 1.}

Improve this answer
$\endgroup$
1
  • $\begingroup$ The RootSearch package <library.wolfram.com/infocenter/MathSource/4482> by Ted Ersek returns all the roots over a specified range as a list of replacement rules This package .A RootsInRange function was presented in "Finding Roots in an Interval" in The Mathematica Journal 7(2), 1998: should be sufficient: $\endgroup$
    – TheDoctor
    Feb 19, 2014 at 7:24
2
$\begingroup$

It seems this problem has been sufficiently answered, but I'd like to mention Stan Wagon's FindAllCrossings2D[].

Improve this answer
$\endgroup$

Your Answer

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

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