1. Bisection algorithm
The algorithm itself is fairly straightforward and "fast" in some sense: the number of iterations is roughly Log2
of the ratio of the initial interval length and the desired accuracy. My point is that the time spent by Flatten
, Select
, Thread
, etc. in your function is fairly small. The significant time-waster is reevaluating the function at the end points of the interval at each iteration. These were evaluated in the previous step and discarded. So you end up with three times as many evaluations of func
as necessary. For a function like
func = #^3 - 4 #^2 - 10 &
which is really fast to evaluate, it's not much of a big deal. But for a complicated function like
func[t_?NumericQ] := 1 + NIntegrate[Sin[x^2] - x, {x, 0, t}]
it makes a big difference.
Since you apparently want something Mathematica-like, as opposed to a procedural loop, here's a take on it. Note, as belisarius commented, we really ought to check whether the exact root has accidentally been found. That adds some case-checking. I'll discuss pattern-checking the arguments later. You may note I defined a realQ
to replace NumericQ
. This may be over-fastidious, but complex numbers will not work.
To address the problem of repeated function evaluations, I added the signs to the data being nested. So instead of just an interval {a, b}
being passed, I now have
{{a, sgna}, {b, sgnb}}
being passed. I used DeleteDuplicatesBy
to delete the old element with the same sign of func
at the midpoint. Since the midpoint comes first, the old one will be deleted. The order of a
and b
does not matter, so splitInterval
ignores it. The top function biSection2
orders the final output.
splitInterval[func_, v : {{a_, sgna_}, {b_, sgnb_}}] :=
With[{sgn = Sign[func[(a + b)/2]]},
If[sgn == 0,
{{(a + b)/2, sgn}, {(a + b)/2, sgn}},
DeleteDuplicatesBy[Join[{{(a + b)/2, sgn}}, v], Last]] (* see Pre V10 note at bottom *)
]
ClearAll[biSection2];
realQ = Quiet@Check[BooleanQ[# < 0], False] &;
biSection2[func_, {a0_?realQ, b0_?realQ}, ξ_?realQ] :=
With[{sgna = Sign[func[a0]], sgnb = Sign[func[b0]]},
Which[
sgna == 0
, {a0, a0},
sgnb == 0
, {b0, b0},
True
, Sort[
First /@
NestWhile[splitInterval[func, #] &, {{a0, sgna}, {b0, sgnb}},
Abs@(Subtract @@ #[[All, 1]]) > ξ &]]
] /; sgna*sgnb <= 0
]
Example, and comparison with OP's biSection
:
Needs["GeneralUtilities`"];
func[t_?NumericQ] := 1 + NIntegrate[Sin[x^2] - x, {x, 0, t}];
(op = biSection[func, {1, 2.`20}, 10^-14]) // AccurateTiming
(me = biSection2[func, {1, 2.`20}, 10^-14]) // AccurateTiming
op - me
(*
0.530108
0.186535
{0.*10^-20, 0.*10^-20}
*)
2. Pattern checking
There are several questions that deal with this issue on the site. Some are listed below. Some deal with related issues such as generating error messages as well as with argument checking. Some deal with the use of Condition
(/;
), which appears at the end of the definition of biSection2
. The ones marked with an asterisk address nearly the same question as the OP's question 2.
Using a PatternTest versus a Condition for pattern matching
_?NumericQ equivalent for lists
How to program a F::argx message?
*More elegant way
*Quick way to use conditioned patterns when defining multi-argument function?
*How shortening argument test when declaring functions?
Let me just make a few remarks about the choices represented above.
If shrinkInterval
or splitInterval
is to be used only internally by biSection
, one might put the responsibility of checking the arguments on biSection
. If the argument checks were significantly time-consuming and shrinkInterval
were iterated many times, it might make sense for the sake of efficiency to do this. This is the approach in the code above, even though ?NumericQ
is very fast and the number of iterations not very great.
The Condition
/; sgna*sgnb <= 0
inside With
in the definition of biSection2
checks the condition after sgna
and sgnb
have been computed and before the Which
statement is evaluated. If the condition is not true, then the evaluation of bisection2
halts and returns the unevaluated expression. See the documentation and the linked questions for more information. Checking the arguments is safer, and should be considered required if the function is meant to be standalone and available to users; however, it is not uncommon to have a user-level function that checks the arguments and then calls an internal function.
Three fun ways to handle myfun
gleaned from the linked questions:
I usually copy-paste-paste-paste for this, but I rarely have more than four or five arguments that need it. Instead of myfun[{x_?NumericQ, y_?NumericQ, z_?NumericQ}, {a_?NumericQ, b_?NumericQ}, {c_?NumericQ, d_?NumericQ}]
, you could...
ClearAll[myfun];
myfun[{x_, y_, z_}, {a_, b_}, {c_, d_}] /;
VectorQ[{x, y, z, a, b, c, d}, NumericQ] :=
Total /@ {{x, y, z}, {a, b}, {c, d}}
ClearAll[myfun, allNumericQ];
SetAttributes[allNumericQ, HoldAllComplete];
allNumericQ[f_[args__]] := VectorQ[Flatten[{args}], NumericQ]; (* caveat Flatten *)
myfun[{x_, y_, z_}, {a_, b_}, {c_, d_}]?allNumericQ :=
Total /@ {{x, y, z}, {a, b}, {c, d}}
ClearAll[myfun];
myfun[{x_, y_, z_}, {a_, b_}, {c_, d_}] /; TrueQ[$inmyfun] :=
Total /@ {{x, y, z}, {a, b}, {c, d}};
myfun[args___] := Block[{$inmyfun = True},
myfun[args]] /; ! TrueQ[$inmyfun] && VectorQ[Flatten[{args}], NumericQ];
All the above work as follows:
myfun[{1, 2, 3}, {4, Sqrt[5]}, {6., 7}]
(* {6, 4 + Sqrt[5], 13.} *)
myfun[{1, x, 3}, {4, Sqrt[5]}, {6., 7}]
(* myfun[{1, x, 3}, {4, Sqrt[5]}, {6., 7}] *)
Caveat: Beware using allNumericQ
if the args
might contain a large amount of data.
Pre V10 users
For users of V9 and earlier, you can use these replacements for DeleteDuplicatesBy
and BooleanQ
. The code for DeleteDuplicatesBy
below is basically Mr.Wizard's myDeDupeBy
in his answer to his question, DeleteDuplicatesBy is not performing as I'd hoped. Am I missing something? Perhaps everyone should be using it.
DeleteDuplicatesBy[x_, f_] := GatherBy[x, f][[All, 1]]
BooleanQ = MatchQ[#, True | False] &