19
$\begingroup$

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}]...
$\endgroup$
2
  • 2
    $\begingroup$ Your algorithm seems to fail if you "accidentally" find an exact root in the bisection process $\endgroup$ Commented Dec 28, 2014 at 14:29
  • $\begingroup$ Something related. Just to note a common problem with bisection methods and how to fix it. $\endgroup$ Commented Jun 1, 2015 at 1:32

5 Answers 5

20
$\begingroup$

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] &
$\endgroup$
2
  • $\begingroup$ Dear Michael E2,Thanks for your detailed answer sincerely:). And I have understood your idea about implementation of Bisection.Owing to that I didn't have MMA 10, so I cannot run the function BooleanQ, DeleteDuplicatesBy and Needs["GeneralUtilities"],So according to your idea, I replaced the DeleteDuplicatesBy[Join[{{(a + b)/2, sgn}}, v], Last] to Delete[Insert[v, {(a + b)/2, sgn}, 2], sgn]]. Do you think it is efficient? $\endgroup$
    – xyz
    Commented Dec 29, 2014 at 2:19
  • 1
    $\begingroup$ @ShutaoTang You're welcome. I added my recommendations for replacements for DeleteDuplicatesBy and BooleanQ. I did not test your Delete version, but I do not think it will work if the function is increasing near the root. (Instead of GeneralUtilities`AccurateTiming, use AbsoluteTiming.) $\endgroup$
    – Michael E2
    Commented Dec 29, 2014 at 2:43
9
$\begingroup$

In the following well known implementation of biSection, the use of shrinkInterval is not needed. Of course your \[Xi] has to be positive. My condition on the input is a little bit stronger and makes the use of NumericQ superfluous.

biSection[ func_, {a_, b_}, eps_?Positive] /; func[a]func[b]<0 :=
   Module[ {c, d, e},
    {c,d}= If[func[a]>0, {b,a}, {a,b}];
    While[d-c>eps,
      With[{e=(c+d)/2},
        If[func[e]<0, c=e, d=e]]
     ];
     {c,d}
   ]

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

(* {1.3652300134140964438,1.3652300134140968879} *)

UPDATE

After almost one and a half year, when someone upvoted this answer, I had another look at this solution and found that it is obviously incorrect when d is less than c. This can trivially be repaired by replacing While[d-c>eps with While[Abs[d-c]>eps, or by using the following variant:

biSection[func_,{a_,b_},eps_?Positive]/;(a<b && func[a]func[b]<0):=
  Module[{c=a,d=b}, With[{test=If[func[a]>0, Greater, Less]},
    While[d-c>eps,With[{e=(c+d)/2},If[test[func[e],0],c=e,d=e]]];
    {c,d}]]

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

(* {1.3652300134140964438,1.3652300134140968879} *)
(* {1.3652300134140964438,1.3652300134140968879} *)
$\endgroup$
0
5
$\begingroup$

I know this question is quite old, but here's my attempt. It's as fast as Michael's, but much simpler:

Bisection[f_, int_, tol_, niter_] := Block[
  {m = tol + 1, prev, ym, yl = f[Last@int]},
  NestWhile[(
     prev = m;
     m = Total@#/2;
     ym = f[m];
     If[ym*yl > 0, yl = ym; {First@#, m}, {m, Last@#}]
     ) &,
   int,  ym != 0 && Abs[m - prev] > tol &, 2, niter]
  ]

Testing it with

func[t_?NumericQ] := 1 + NIntegrate[Sin[x^2] - x, {x, 0, t}];

Bisection[func, {1, 2.`20}, 10^-14, 1000] // Timing
{0.109375, {1.9252809180739163253, 1.9252809180739234307}}
$\endgroup$
0
$\begingroup$
Attributes@bisection = {HoldAll};
bisection::nochangeofsign = "sign does not change, try other bounds";
Options@bisection = {AccuracyGoal -> $MachinePrecision/2, 
   Rational -> False};
Clear@bisection; 
bisection[expr_, {x_, x1_, x2_}, OptionsPattern[]] := 
 With[{f1 = Function[x, #] &@expr},
  With[{s = 
     Which[LessEqual @@ #, 1, GreaterEqual @@ #, -1, True, 
        Message@bisection::nochangeofsign] &@{f1@x1, 0, f1@x2}},
   If[Not@OptionValue@Rational, 1., 1] First@
     NestWhile[
      If[s f1@#2 > 0, {First@#, #2}, {#2, Last@#}] &[#, 
        Total@#/2] &, {x1, x2}, 
      Abs[First@# - Last@#] > 
        Power[10, -OptionValue@AccuracyGoal] &]]]

I think it is not good to use constants or variables in the core of the program: my experience is that it is slower (and it is also ugly). Please someone do the Timing and try to improve.

As the result is Rational, it is a good idea to keep Rational in Options just in case the user needs Rational approximation.

Usage:

bisection[Cos@x - x, {x, 0, 1.5}]

bisection[-Cos@x + .8, {x, 0, 1}, AccuracyGoal -> 10, 
 Rational -> True]

This is a funny example:

bisection[Tan[x/4] - 1, {x, 3, 4}, AccuracyGoal -> 20, 
 Rational -> True]
$\endgroup$
1
  • $\begingroup$ I think that the NumericQ question is solved by the HoldAll. $\endgroup$ Commented Apr 22, 2020 at 9:09
0
$\begingroup$

My program is not bad but has design flaws: the Rational option is useless and the interface will fail in certain cases. So I need to update:

Attributes@bisection = {HoldAll};
bisection::nochangeofsign = "sign does not change, try other bounds";
Options@bisection = {AccuracyGoal -> $MachinePrecision/2, 
   Sow -> False};
Clear@bisection; 
bisection[expr_, {x_, x1_, x2_}, OptionsPattern[]] := 
 Catch@If[OptionValue@Sow, {First@#, First@Last@#} &, First]@
   Reap@With[{sign = 
       Which[LessEqual @@ #, 1, GreaterEqual @@ #, -1, True, 
          Message@bisection::nochangeofsign] &@(*Throw@*){ReleaseHold[
          Hold@expr /. x -> x1], 0, 
         ReleaseHold[Hold@expr /. x -> x2]}},
     If[OptionValue@Sow, #, First@#] &@
      With[{hook = If[OptionValue@Sow, Sow, Identity]}, 
       NestWhile[
        If[sign ReleaseHold[Hold@expr /. x -> hook@#2] > 
             0, {First@#, #2}, {#2, Last@#}] &[#, Total@#/2] &, {x1, 
         x2}, Abs[First@# - Last@#] > 
          Power[10, -OptionValue@AccuracyGoal] &]]]

Usage:

bisection[-Cos@x + x, {x, 0, 1.5}, Sow -> True]

This example is fun:

bisection[x*x - 2, {x, 1, 2}, AccuracyGoal -> 20]

Comparing this version with the former does not teach anything on numerical analysis but it teaches a bit on the Wolfram language.

$\endgroup$

Your Answer

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

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