Using FindRoot with an interpolating function

Question

This question bears resemblance to a few other questions on mathematica.SE about finding points of intersection of crossing curves. I know that the guidebook of numerics has an entry about the whole curve crossing thingy (can't seem to find the link right now).

However, my question is a little different. Yes, crossing curves are involved.

I have two curves that cross each other at two points:

curve1 = 3 x^2 + 3 x;
curve2 = 1.8 x ^2 + 2;
Plot[
 {curve1, curve2},
 {x, -5, 5},
 PlotRange -> All
 ]

Curves crossing

With Roots[...] I can find the points at which these curves cross each other, so:

Roots[curve1 == curve2, x]

x==-3.04699||x==0.546988

So this is nice and happy! Now, if I were to get data out of the individual plots, interpolate this data and fold it into and InterpolatingFunction, I am unable to use FindRoot[...] to do the same as Root[...]

pic1 = Plot[curve1, {x, -5, 5}];
Data1 = Cases[Normal[pic1], Line[Data1_] :> Data1, Infinity];
intplC1 = Data1 // Flatten // Interpolation
pic1 = Plot[curve2, {x, -5, 5}];
Data2 = Cases[Normal[pic1], Line[Data2_] :> Data2, Infinity];
intplC2 = Data2 // Flatten // Interpolation


FindRoot[intplC1 == intplC2, {x, 0.2}]

FindRoot::nlnum: The function value {InterpolatingFunction[{{1.,528.}},{4,7,0,{528},{4},0,0,0,0,Automatic},{{<<1>>}},{Developer`PackedArrayForm,{<<1>>},{-5.,60.,-4.99693,59.9172,<<43>>,3.80528,-1.7292,3.78277,<<478>>}},{Automatic}]-<<1>>} is not a list of numbers with dimensions {1} at {x} = {0.2}. >>

So my question(s) are:

I am thinking I am not using Cases[...] correctly here despite the fact that Data1//ListLinePlot and Data2//ListLinePlot seem to plot fine enough.
How can I use FindRoot[...] on my InterpolatingFunctions to find the multiple roots in this case?
For this situation, am I right to assume that Roots[...] is well and sufficient?

Michael E2 · Accepted Answer · 2013-04-19T20:43:34.967

8

Try

pic1 = Plot[curve1, {x, -5, 5}];
Data1 = Cases[Normal[pic1], Line[Data1_] :> Data1, Infinity];
intplC1 = Flatten[Data1, 1] // Interpolation
pic1 = Plot[curve2, {x, -5, 5}];
Data2 = Cases[Normal[pic1], Line[Data2_] :> Data2, Infinity];
intplC2 = Flatten[Data2, 1] // Interpolation

FindRoot[intplC1[x] == intplC2[x], {x, 0.2}]

Flatten in your original code, completely flattens the coordinates. It becomes a list

{ x1, y1, x2, y2, ... }

Further, the argument x must be supplied to the interpolating function.

Addendum

If you have interpolation from data, then the following could be used to find most intersections.

SeedRandom[2];
xBase = Sort[RandomReal[{0, 10}, 10]];
data1 = Table[{xBase[[i]], RandomReal[{0, 1}]}, {i, 10}];
intplC1 = data1 // Interpolation;
data2 = Table[{xBase[[i]], RandomReal[{0, 1}]}, {i, 10}];
intplC2 = data2 // Interpolation;

x0 = Pick[MovingAverage[Data1[[All, 1]], 2], 
       Negative /@ Times @@@ Partition[
         Subtract @@@ Transpose@{data1[[All, 2]], data2[[All, 2]]}, 2, 1]]

FindRoot[intplC1[x] == intplC2[x], {x, x0}]

(* {1.37302, 2.29548, 5.5938, 5.92218, 7.6267} *)

(* {x -> {1.13464, 2.72027, 5.38934, 5.9991, 7.65357}} *)

Plot[{intplC1[x], intplC2[x]}, Evaluate@{x, Sequence @@ intplC1[[1, 1]]}]

Plot of interpolating functions

Notes:

1) You have to pass multiple initial points inside its own list (x0 = {1.13464, ...}).

2) Interpolations can have multiple intersections between interpolated points, and the method above ignores that possibility. In that case, starting points could be added manually.

edited Apr 19 '13 at 20:43

answered Apr 19 '13 at 12:36

Michael E2

235,386
17
334
747

1

I was just going to answer about the same :( – swish Apr 19 '13 at 12:41
I still have an error message that suggests that the "function value is not a list of numbers" at x=0.2 – dearN Apr 19 '13 at 12:42
1

In addition, you can supply a list of starting values to FindRoot, like {x, {0.2, -3.}} and booth roots will be found. – BoLe Apr 19 '13 at 12:45
@MichaelE2 I get what you mean now when you say x must be provided to the interpolation. – dearN Apr 19 '13 at 12:51
@BoLe What if I don't know what the guess values should be and how many there should be? – dearN Apr 19 '13 at 12:51
@drN Try it now. Cut and paste had missed a letter. – Michael E2 Apr 19 '13 at 12:53
@MichaelE2 I modified my code. I used Interpolation[Flatten[curve1,1],x]. This worked. I forgot that I needed to tell it to use x :P – dearN Apr 19 '13 at 12:54
@BoLe And to answer my own question, I just used the minimum and maximum x limits as my guess values. I think that does the trick. – dearN Apr 19 '13 at 12:54
@drN It's a question of approaching I guess. At least in case of a clean rootfinding off a plot, I see no problem. – BoLe Apr 19 '13 at 12:59
What if I were to have more than one crossing and I still wanted to use FindRoot[...]. I noticed that FindRoot[...] doesn't allow for more than 2 guess values. I can use an {x,xmin,xstart,xend} but that didn't quite help. – dearN Apr 19 '13 at 14:34

score 4 · Answer 2 · answered Apr 19 '13 at 12:55

4

Why Roots in the first place? Also, why not define curves as pure functions?

f = 3 #^2 + 3 # &;
g = 1.8 #^2 + 2 &;

FindRoot[f@x - g@x, {x, {-3, 1}}]

(* {x -> {-3.04699, 0.546988}} *)

Supply a list of initial values and you can get more different solutions at the same time.

Coordinates of intersections:

{x, f@x} /. FindRoot[f@x - g@x, {x, {-3, 1}}]

(* {{-3.04699, 0.546988}, {18.7114, 2.53855}} *)

answered Apr 19 '13 at 12:55

BoLe

5,819
15
33

Very elegant! Nice – Sos Apr 19 '13 at 13:31

score 2 · Answer 3 · answered Apr 04 '19 at 02:19

You can use NDSolve to find the intersection points of two interpolating functions, similar to the accepted answer to How to find all the local minima/maxima in a range. Using Michael's example:

SeedRandom[2];
xBase = Sort[RandomReal[{0, 10}, 10]];
data1 = Table[{xBase[[i]], RandomReal[{0, 1}]}, {i, 10}];
intplC1 = data1 // Interpolation;
data2 = Table[{xBase[[i]], RandomReal[{0, 1}]}, {i, 10}];
intplC2 = data2 // Interpolation;
Plot[{intplC1[t], intplC2[t]}, {t, 1.1, 7.7}]

Finding the intersection points using NDSolve:

ints = Reap[
    NDSolve[
        {
        y'[t] == intplC1'[t] - intplC2'[t],
        y[2] == intplC1[2] - intplC2[2],
        WhenEvent[y[t]==0, Sow[t]]
        },
        y,
        {t, 1.1, 7.7}
    ]
][[2, 1]]

{1.13464, 2.72027, 5.38934, 5.9991, 7.65357}

Visualization:

Plot[
    {intplC1[t], intplC2[t]},
    {t, 1.1, 7.7},
    Epilog -> {PointSize[Large], Red, Point[Thread[{ints, intplC1[ints]}]]}
]

Using FindRoot with an interpolating function

3 Answers3

Linked