# Find all roots of an interpolating function (solution to a differential equation)

I'm trying to find all the roots of the solution to a differential equation. Using NSolve or Reduce I don't get the roots, so I'm using an iterative method which I found in physicsforums.com. This method solves my problem, but you have to choose the increment and thus in some cases it might give headaches. I wonder if there is any more general approach.

Following is a sample code:

data = NDSolve[{1.09 x''[t] - 0.05 x'[t] + 1.1759 Sin[x[t]] == 0,
x == Pi/3, x' == 0}, x, {t, 0, 50}]


{{x->InterpolatingFunction[{{0.,50.}},<>]}}

First attempt with no success:

sol = NSolve[x'[t] == 0 /. data , t]


NSolve::ifun: Inverse functions are being used by NSolve, so some solutions may not be found; use Reduce for complete solution information. >>

{{t->InverseFunction[InterpolatingFunction[{{0.,50.}},<>],1,1][0.]}}

Second attempt with no success:

sol = Reduce[x'[t] == 0 /. data , t]


Reduce::inex: Reduce 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 Reduce require exact input, providing Reduce with an exact version of the system may help. >>

Reduce[{InterpolatingFunction[{{0.,50.}},<>][t]==0},t]

Third attempt, works fine, but manually choosing dt could cause problems with some equations:

dt = 0.1;
tmin = 0.;
tmax = 50.;

Union[Table[t /. FindRoot[x'[t] == 0 /. data ,{t, tInit, tmin, tmax}],
{tInit, tmin + dt,tmax - dt, dt}], SameTest->(Abs[#1 - #2] < 10^-2&)]


{0., 3.26812, 6.58301, 9.95657, 13.4054, 16.9538, 20.6403, 24.533, 28.7857, 34.2571}

Is there any more elegant method to find all roots in a range?

• Dec 16, 2012 at 23:23

I'm slightly surprised nobody's already mentioned the event location capabilities of Mathematica, as it's the most compact way to find the roots of an interpolating function that came from NDSolve[]. I don't have Mathematica on this machine I'm writing in, but I'd do something like this:

Reap[NDSolve[{1.09 x''[t] - 0.05 x'[t] + 1.1759 Sin[x[t]] == 0, x == Pi/3, x' == 0},
x, {t, 0, 50}, Method -> {"EventLocator",
"Event" -> x[t], "EventAction" :> Sow[t]}]]


which should yield the interpolating function and the list of roots.

At least, that's how its done in version 8. Version 9 has the WhenEvent[] function, which can be used like so:

Reap[NDSolve[{1.09 x''[t] - 0.05 x'[t] + 1.1759 Sin[x[t]] == 0, x == Pi/3, x' == 0,
WhenEvent[x[t] == 0, Sow[t]]}, x, {t, 0, 50}]]

• I did mention this method in my answer above by referring to Daniel's answer here, of which this answer is a variation.
– Jens
Dec 17, 2012 at 4:37
• Jens: Well, since the Brent algorithm used by the event locator can take options, you can still easily leverage the abilities of FindRoot[]: Method -> {"EventLocator", "Event" -> x[t], "EventAction" :> Sow[t], "EventLocationMethod" -> {"Brent", AccuracyGoal -> Infinity, "SolutionApproximation" -> "CubicHermiteInterpolation"}}. There should be an equivalent way of tweaking such options if you're taking the WhenEvent[] route. Dec 17, 2012 at 4:43
• (I'm quite sure a question like this has been asked before on this site, but I suck at searching for dupes.) Dec 17, 2012 at 4:44

Here is a function aptly named findAllRoots that is based on idea published a long time ago in the Mathematica Journal, I believe. I'll try to find the article again, but in the meantime here is my version of it. I added the option handling, in particular the possibility to show the plot of the function while finding its roots. The idea to extract data from a Plot (thereby leveraging its analysis of the curve) and to then use Split in this non-standard way is the main thing that I remember from the article.

Clear[findAllRoots]
SyntaxInformation[
findAllRoots] = {"LocalVariables" -> {"Plot", {2, 2}},
"ArgumentsPattern" -> {_, _, OptionsPattern[]}};
SetAttributes[findAllRoots, HoldAll];

Options[findAllRoots] =
Join[{"ShowPlot" -> False, PlotRange -> All},
FilterRules[Options[Plot], Except[PlotRange]]];

findAllRoots[fn_, {l_, lmin_, lmax_}, opts : OptionsPattern[]] :=
Module[
{pl, p, x, localFunction, brackets},
localFunction = ReleaseHold[Hold[fn] /. HoldPattern[l] :> x];
If[
lmin != lmax,
pl = Plot[localFunction, {x, lmin, lmax},
Evaluate@
FilterRules[Join[{opts}, Options[findAllRoots]], Options[Plot]]
];
p = Cases[pl, Line[{x__}] :> x, Infinity];
If[OptionValue["ShowPlot"],
Print[Show[pl, PlotLabel -> "Finding roots for this function",
ImageSize -> 200, BaseStyle -> {FontSize -> 8}]]],
p = {}
];
brackets = Map[
First,
Select[
(* This Split trick pretends that two points on
the curve are "equal" if the function
values have _opposite _ sign. Pairs of such
sign-changes form the brackets for the subsequent
FindRoot *)
Split[p, Sign[Last[#2]] == -Sign[Last[#1]] &],
Length[#1] == 2 &
],
{2}
];
x /. Apply[FindRoot[localFunction == 0, {x, ##1}] &,
brackets, {1}] /. x -> {}
]


You can specify PlotPoints as an option to findRoots, in order to increase the resolution of the initial search for zero crossings which are then used as brackets for FindRoot. This way, if you initially don't find all roots with the default settings, you can improve the completeness of the result at the expense of longer execution time.

End note

Here follows the interpolation function from your question, and the application of the above function:

data = NDSolve[{1.09 x''[t] - 0.05 x'[t] + 1.1759 Sin[x[t]] == 0,
x == Pi/3, x' == 0}, x, {t, 0, 50}];

f[t_] = (x /. First[data])[t];

findAllRoots[f[t], {t, 0, 50}, "ShowPlot" -> True] {1.59975,4.88823,8.22838,11.6339,15.1242,18.7282,22.4934,26.5081,30.9882,37.4065}

I've used this successfully before, and I've also compared it to Daniel's solution in this answer. I prefer the above method because his NDSolve approach ended up being less accurate in my applications.

Edit: IntervalRoots

One of the standard add-on packages that disappeared during the upgrade from Mathematica version 5.2 to 6 was NumericalMath. If you had a notebook that started with

<<NumericalMathIntervalRoots


then you could call a function like

IntervalNewton[Sin[t], t, Interval[{10, 20}], .1]


Interval[{12.5661,12.5664},{15.6503,15.7438},{18.8495,18.8506}]

Now if you want to try replacing the root bracketing in my function findAllRoots by this bracketing routine, you have to do get this package manually at MathSource. I thought I'd mention this for completeness. However, these old bracketing functions apparently can't handle InterpolatingFunction properly, so this is purely a historical footnote.

• Very nice, thanks a lot for your effort! Jun 24, 2015 at 15:39
• Do you mind if I include your function in a package? (for finding eigenvalues of 1D boundary value problems) Mar 29, 2018 at 11:10
• @KraZug Go ahead and include it - I'm glad it proved useful.
– Jens
Mar 29, 2018 at 15:58
• Hello, Jens. Thank you for making this! If a function F is complex-valued, is there a better way than checking both Re@F and Im@F with FindAllRoots, given that Abs@F looks incompatible with your algorithm? I ask this because I want to avoid doing the same calculation twice as F might be a heavy calculation. Jul 13, 2019 at 19:43
• @Jens In the End Note code, if you define g[t_]=f[t]-Pi/3, then 0 should be a root, but when I run your code it doesn't show 0 as a root. May 10, 2020 at 13:03

Using RootSearch.m (mentioned in the above link given by MrWizard)

http://library.wolfram.com/infocenter/Demos/4482/

download the .m file, put it in my current folder and did:

SetDirectory[NotebookDirectory[]];
Get["RootSearch.m"];

eq = 1.09 x''[t] - 0.05 x'[t] + 1.1759 Sin[x[t]] == 0;
ic = {x == Pi/3, x' == 0};
sol = First@NDSolve[Flatten[{eq, ic}], x[t], {t, 0, 50}];
ErsekRootSearchRootSearch[Evaluate[x[t] /. sol] == 0, {t, 0, 50}]

(*  {{t -> 1.59975}, {t -> 4.88823}, {t -> 8.22838}, {t -> 11.6339},
{t -> 15.1242}, {t -> 18.7282}, {t -> 22.4934}, {t -> 26.5081},
{t ->    30.9882}, {t -> 37.4065}}  *)


ps. on console I see this warning message

$MinPrecision::precset: Cannot set$MinPrecision to -\[Infinity]; value
must be a non-negative number or Infinity. >>


But the package has options. So these things migthbe possible to configure. see the notebook for more information • I found this as a possible cure for the MinPrecision error: set it to +Infinty.
– Jens
Dec 17, 2012 at 2:27
• I am truly overwhelmed by the responses and their quality. Thanks. It has helped me a lot and I hope they will also help more people. Not only because of the answers but studying the submitted codes I'm discovering new possibilities of Mathematica. Thanks. Dec 17, 2012 at 22:39

You can extract values for the roots directly from the the Plot by using the MeshFunctions option.

f = x /. First@NDSolve[
{1.09 x''[t] - 0.05 x'[t] + 1.1759 Sin[x[t]] == 0, x == Pi/3, x' == 0},
x, {t, 0, 50}
];


The default plot options give reasonable approximations to the roots:

plot = Plot[f[t], {t, 0, 50},
Mesh -> {{0}}, MeshFunctions -> {#2 &},
MeshStyle -> PointSize[Medium], PlotRange -> All
] Sort@Cases[Normal@plot, Point[{x_, y_}] -> x, Infinity]

(*{1.59974, 4.88815, 8.22838, 11.6334, 15.1239, 18.7279, 22.4934, 26.5078, 30.9879, 37.4064}*)


You can get better approximations by increasing the PlotPoints:

plot = Plot[f[t], {t, 0, 50},
Mesh -> {{0}}, MeshFunctions -> {#2 &},
MeshStyle -> PointSize[Medium],
PlotRange -> All, PlotPoints -> 2000
];
Sort@Cases[Normal@plot, Point[{x_, y_}] -> x, Infinity]

(*{1.59975, 4.88822, 8.22838, 11.6339, 15.1242, 18.7282, 22.4934, 26.5081, 30.9882, 37.4065}*)

• and if the accuracy is not sufficient, the found roots can be used as successive initial conditions in FindRoot. Mar 4, 2017 at 23:47