Skip to main content

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

Comments

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...