equation solving - Implement the Bisection algorithm elegantly and easily

Description:

Rencently, I have finished my course Numerical Analysis, so I'd like to implement many algorithm that I have learned from that course.By this practice, I hope that I can improve my programming skill and understand the knowledge of numerical analysis deeply.

Firstly,I implement the bisection to search the root of nonlinear equation

My trial

shrinkInterval[func_, {a_?NumericQ, b_?NumericQ}] /; a < b :=
 Module[{interval},
  interval = {a, (a + b)/2, b};

  Last@
   Flatten[
    Select[
     Thread@
     (Partition[#, 2, 1] & /@
      {Sign /@ (func /@ interval), interval}), 
     Times @@ (First@#) < 0 &], 1]
 ]
(*=====================================================*)
biSection[func_, {a0_?NumericQ, b0_?NumericQ}, \[Xi]_?NumericQ] :=

 NestWhile[
  shrinkInterval[func, #] &, {a0, b0}, Abs@(Subtract @@ #) > \[Xi] &]

Testing

biSection[#^3 + 4 #^2 - 10 &, {1, 1.5`20}, .5*10^-15]

{1.3652300134140964438, 1.3652300134140968879}

My question:

1, In my code, I think the function shrinkInterval is a little complex. I'd like to know more simple method to rewrite it?

2, Is there any good form/strategy to deal with the below condition that I need to write NumericQ many times?

myfun[{x_?NumericQ, y_?NumericQ,z_NumericQ},{a_?NumericQ, b_?NumericQ},
  {c_?NumericQ, d_?NumericQ}]...

Answer

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] &