(*********************************************************************** This file was generated automatically by the Mathematica front end. It contains Initialization cells from a Notebook file, which typically will have the same name as this file except ending in ".nb" instead of ".m". This file is intended to be loaded into the Mathematica kernel using the package loading commands Get or Needs. Doing so is equivalent to using the Evaluate Initialization Cells menu command in the front end. DO NOT EDIT THIS FILE. This entire file is regenerated automatically each time the parent Notebook file is saved in the Mathematica front end. Any changes you make to this file will be overwritten. ***********************************************************************) BeginPackage["Ersek`RootSearch`"]; RootSearch::usage= "RootSearch[lhs==rhs,{x,xmin,xmax}] tries to find all numerical solutions \ to the equation (lhs==rhs) with values of x between xmin and xmax."; RootSearchSamples::usage= "RootSearchSamples returns a list of points where the equation given to \ RootSearch was sampled during the last use of RootSearch. If \ RootSearch[lhs==rhs,{x,xmin,xmax}] was evaluated. Then RootSearchSamples \ returns {{x1,y1},{x2,y2}, ...} where yn is the value of (lhs-rhs) at x=xn. Points where (lhs-rhs) did not evaluate to a numeric value are not included \ in the list returned."; MaxBrentSteps::usage= "MaxBrentSteps is an option for RootSearch which specifies the maximum \ number of steps the algorithm should take when attempting to find a specific \ root using Brent's method."; MaxSecantSteps::usage= "MaxSecantSteps is an option for RootSearch which specifies the maximum \ number of steps the algorithm should take when attempting to find a specific \ root using a modified secant method."; RootTest::usage= "RootTest is an option for RootSearch which determines whether specific \ values are roots of an equation. Once the core algorithm for \ RootSearch[lhs==rhs,{x,xmin,xmax}] is finished a list of points \ {{x1,y1},{x2,y2}, ...} is effectively made where each xi is a likely root and \ yi is the value of (lhs-rhs) when (x=xi). Only the likely roots for which \ RootTest is True are returned. The RootTest setting should be a function of \ one or two arguments where xi is the first argument and yi is the second \ argument."; InitialSamples::usage= "InitialSamples is an option for RootSearch which specifies the number of \ nearly equally spaced samples that should be taken before root finding \ methods are used to converge on roots."; InitialPrecision::usage= "InitialPrecision is an option for RootSearch which specifies the \ precision of the values where a function is intially sampled. If the \ InitialPrecision setting is a positive value less than or equal to \ $MachinePrecision, then machine numbers are used. If the setting is larger \ than $MachinePrecision the setting specifies the level of arbitrary precision \ used."; (*GuardDigits::usage="A RootSearch option.";*) Begin["`Private`"]; Unprotect[RootSearch,InitialSamples,RootTest,MaxBrentSteps,MaxSecantSteps, PrecisionGoal,InitialPrecision]; RootSearch::iopnm= "The setting of the `1` option in RootSearch must be a NonNegative Integer."; RootSearch::prec= "The setting of the RootSearch option `1` must be a Positive number less than or equal to $MaxPrecision."; RootSearch::var= "RootSearch was given an an equation involving `1` when it expected to solve an equation in `2`."; RootSearch::vars= "RootSearch was given an an equation with variables (`1`) when it expected to solve an equation where the only variable is `2`."; RootSearch::smpls= "The setting of the InitialSamples Option in RootSearch must be an Integer between 2 and 2147483648."; RootSearch::pllim= "Range specification `1` is not in the form {x, xmin, xmax}. Neither xmin or xmax can have infinite magnitude."; RootSearch::eqf="`1` is not a well-formed equation."; RootSearch::optx="`1` are not recognized as RootSearch Options."; RootSearch::lowp="One or more RootSearch sample resulted in low precision."; RootSearch::numb="RootSearch took a number of inital samples and at each sample point the function sampled did not evaluate to a numeric value."; RootSearch::err="The Private RootSearch function `1` had an error in branch `2` when searching for a root near `3`. Please email a notebook containing this use of the RootSearch function to RootSearch author at (ted.ersek@navy.mil)."; RootSearch::numer="The Private RootSearch function `1` sampled the given equation at `2` and the result was not numeric."; Ulp2::fail="The private function Ulp2 in the RootSearch package was given two equal arguments which is not allowed. Please email a notebook containing this use of the RootSearch function to RootSearch author at (ted.ersek@navy.mil)."; RootSearch::pole="It seems the Private RootSearch function zBrent is converging on a pole near `1`."; Off[RootSearch::pole]; (*RootSearch::grdg="The setting for Options[RootSearch,GuardDigits] must be numeric and must not have the head Complex.";*) Options[RootSearch]={InitialSamples\[Rule]300, RootTest\[RuleDelayed] (Abs[#2]\[LessEqual]10^5(Abs[#1]+$MinMachineNumber)$MachineEpsilon&), MaxBrentSteps\[Rule]150,MaxSecantSteps\[Rule]150, PrecisionGoal\[RuleDelayed]$MachinePrecision, InitialPrecision\[RuleDelayed]$MachinePrecision}; RootSearch[ lhs_\[Equal]rhs_,{x_Symbol?(!NumericQ[#]&),xmin1_?NumericQ, xmax1_?NumericQ},opts___?OptionQ]/; Im[N[xmin1,17]]\[Equal]0&&Im[N[xmax1,17]]\[Equal]0&&Re[xmin1]$MaxPrecision,invalidQ=True; Message[RootSearch::prec,PrecisionGoal]]; If[(PositiveQ[InitPrec]===False)||InitPrec>$MaxPrecision,invalidQ=True; Message[RootSearch::prec,InitialPrecision]]; If[!invalidQ,{PrecGoal,InitPrec}={Max[N@PrecGoal,$MachinePrecision], Max[N@InitPrec,$MachinePrecision]}]; If[PrecGoal\[Equal]InitPrec\[Equal]$MachinePrecision, GrdDigts=0;{xmin,xmax}\[Equal]N[{xmin,xmax}],(*else*)GrdDigts= 10;{xmin,xmax}=N[{xmin,xmax},InitPrec]]; If[opts=!=Null, BadOpts=Complement[Part[Flatten[{opts}],All,1], Part[Options[RootSearch],All,1],SameTest\[Rule]Automatic]; If[Length[BadOpts]===1, Message[General::optx,First[BadOpts],RootSearch]]; If[Length[BadOpts]>1,Message[RootSearch::optx,ToString[BadOpts]]]]; Clear[SamplePoints]; If[(!TrueQ@IntegerQ@InitSmpls)||InitSmpls<2||InitSmpls>2147483647, invalidQ=True;Message[RootSearch::smpls]]; If[(!TrueQ@IntegerQ@MaxBrntStps)||Negative@MaxBrntStps,invalidQ=True; Message[RootSearch::iopnm,MaxBrentSteps]]; If[(!TrueQ@IntegerQ@MaxSecntStps)||Negative@MaxSecntStps,invalidQ=True; Message[RootSearch::iopnm,MaxSecantSteps]]; SymbolList= Union@Cases[{lhs, rhs}/.{_Rule\[Rule]Sequence[],_RuleDelayed\[Rule] Sequence[]},_Symbol?(!NumericQ[#]&& Context[#]=!="System`"&),{-1},Heads\[Rule]False]; If[Complement[SymbolList,{x}]=!={},invalidQ=True; If[Length[SymbolList]===1,Message[RootSearch::var,First@SymbolList,x], Message[RootSearch::vars,StringTake[ToString[SymbolList],{2,-2}],x]] ]; If[invalidQ,$Failed, (*else*) f=Evaluate[{SamplePoints[#],lhs-rhs/.x\[Rule]#}]&; f=ReplacePart[f,Set,{1,0}]; u=Function[ Evaluate[ withTemp[{setTemp[y,(lhs-rhs)/.x\[Rule]#], setTemp[yyp,((lhs-rhs)/D[lhs-rhs,x])/.x\[Rule]#]}, Hold[SamplePoints[#]=y, Which[TrueQ[y\[Equal]0],0,NumericQ[yyp],yyp,True, Indeterminate]]]]]; u=ReplacePart[ u/.{withTemp\[Rule]With,setTemp\[Rule]Set},{CompoundExpression},{{1, 2,0}},{{1}}]; If[TrueQ@FreeQ[u,Derivative,Heads\[Rule]True], roots= SearchForRootsU[f,u,xmin,xmax,MaxSecntStps,MaxBrntStps,InitSmpls, InitPrec,PrecGoal,GrdDigts], (*else*) roots=SearchForRootsF[f,xmin,xmax,MaxSecntStps,MaxBrntStps,InitSmpls, InitPrec,PrecGoal,GrdDigts] ]; Block[{$MinPrecision=-\[Infinity]}, If[MemberQ[ RootSearchSamples,_?((Precision[#]<3)&&(Abs[#]> 1000$MachineEpsilon)&),{-1}], Message[RootSearch::lowp]]]; If[roots==={},{}, (*else we move on*) roots=Select[Sort[roots],Test[#,f[#]]&]; If[PrecGoal\[NotEqual]$MachinePrecision, (*Bring Precision back down.*)roots= Map[If[TrueQ[(#\[Equal]0)&&(Accuracy[#]>PrecGoal)], SetAccuracy[#,PrecGoal],#]&,N[roots,PrecGoal], Heads\[Rule]False],(*else*)roots=N[roots] ]; Map[{x\[Rule]#}&,roots,Heads\[Rule]False] ] ] ] RootSearch[_,_,args__]/;!OptionQ[Flatten@{args}]&& Message[RootSearch::nonopt,First@Select[Flatten@{args},!OptionQ[#]&,1], 2,RootSearch[__]]:="Never get here" RootSearch[_,arg2_,___]/;(Head[arg2]=!=List||Length[arg2]=!=3)&& Message[RootSearch::pllim,arg2]:="Never get here" RootSearch[_,{symb_,_,_},___]/;(Head[symb]=!=Symbol||NumericQ[symb])&& Message[RootSearch::itraw,symb]:="Never get here" HoldPattern[ RootSearch[arg1_,___]/; Head[arg1]=!=Equal&&Message[RootSearch::eqf,arg1]]:="Never get here" HoldPattern[ RootSearch[_, arg2:{_,xmin_,xmax_},___]/;(Im[N[xmin,17]]\[NotEqual]0|| Im[N[xmax,17]]\[NotEqual]0||Re[xmax]\[LessEqual]Re[xmin]|| Not[NumericQ[xmin]]||Not[NumericQ[xmax]])&& Message[RootSearch::pllim,arg2]]:="Never get here" SearchForRootsU= Function[{f,u,xmin,xmax,MaxSecantIter,MaxBrentIter,InitialSamples, InitPrec,PrecGoal,GrdDigts}, Module[{xList,yList,ZeroPosn,CrossingPosn1,CrossingPosn2,CrossingVals1, CrossingVals2,MinPos1,MinPos2,soln1,soln2,soln3,soln4,soln5,xyPnts, xPnts,yPnts,xs,signs}, xList=InitialPoints[xmin,xmax,InitialSamples,InitPrec,PrecGoal]; yList= Developer`ToPackedArray[MakeReal[Map[f,xList,Heads\[Rule]False]]]; If[FreeQ[yList,_?NumericQ,{1}], Message[RootSearch::numb];{}, (*else we move on.*) ZeroPosn=Position[yList,_?(#\[Equal]0&),Heads\[Rule]False]; If[PrecGoal\[Equal]$MachinePrecision, soln1=Extract[xList,ZeroPosn], (*else*) soln1=SetPrecision2[Extract[xList,ZeroPosn],PrecGoal+GrdDigts]; zeroList=MakeReal[Map[f,soln1,Heads\[Rule]False]]; ZeroPosn=Position[zeroList,_?(#\[Equal]0&),Heads\[Rule]False]; soln1=Extract[soln1,ZeroPosn] ]; (*Find the zero crossings of f[x].*) CrossingPosn1= Flatten@Position[ Partition[yList,2,1],_List?(NegativeQ[First[#]*Last[#]]&),{1}, Heads\[Rule]False]; CrossingPosn1= Complement[CrossingPosn1,Flatten[ZeroPosn], SameTest\[Rule]Automatic]; CrossingPosn1=Transpose[{CrossingPosn1,CrossingPosn1+1}]; (*Next use CrossingPosition[] to find positions where both f[x], u[x] change sign.*) CrossingPosn2=CrossingPosition[xList,CrossingPosn1,u]; If[CrossingPosn2=!={}, CrossingVals2= Map[Take[xList,#]&,CrossingPosn2,Heads\[Rule]False]; If[PrecGoal\[NotEqual]$MachinePrecision, CrossingVals2=SetPrecision2[CrossingVals2,PrecGoal+GrdDigts]]; (*Apply zBrent on u[x] at the zero crossings of f[x],u[x].*) soln2=Apply[zBrent[u,##,MaxBrentIter,InitPrec,PrecGoal+GrdDigts]&, CrossingVals2,1,Heads\[Rule]False], (*else*)soln2={} ]; (*Prepare to apply zBrent on f[x] where f[x] has a zero crossing, but u[x] does not.*) CrossingPosn1= Complement[CrossingPosn1,CrossingPosn2,SameTest\[Rule]Automatic]; If[CrossingPosn1=!={}, CrossingVals1= Map[Take[xList,#]&,CrossingPosn1,Heads\[Rule]False]; If[PrecGoal\[NotEqual]$MachinePrecision, CrossingVals1=SetPrecision2[CrossingVals1,PrecGoal+GrdDigts]]; soln3=Apply[zBrent[f,##,MaxBrentIter,InitPrec,PrecGoal+GrdDigts]&, CrossingVals1,1,Heads\[Rule]False], (*else*)soln3={} ]; (*Find samples near a local min of Abs[f[x]].*) MinPos1=Position[Partition[yList,3,1],_?LocalExtremaQ,{1}, Heads\[Rule]False]+1; (*Prepare to Apply zBrent on u[ x] at samples near the local min of Abs[f[x]] where u[ x] changes sign.*) If[MinPos1=!={}, xs=Apply[Take[xList,{#-1,#+1}]&,MinPos1,{1},Heads\[Rule]False]; signs=Map[Sign2[MakeReal[u[#]]]&,xs,{-1},Heads\[Rule]False]; (*MinPos2 is the position of samples near local min of Abs[ f[x]] where u[ x] has opposite sign at the first and last samples.*) MinPos2=Extract[MinPos1, Position[ Replace[ signs,{{1,1,-1}\[Rule]0,{-1,1,1}\[Rule]0,{-1,-1, 1}\[Rule]0,{1,-1,-1}\[Rule]0},{-2}],_List,1]]; xMinVals= Transpose[{xs, Replace[ signs,{{1,1,-1}\[Rule]{2,3},{-1,1,1}\[Rule]{1,2},{-1,-1, 1}\[Rule]{2,3},{1,-1,-1}\[Rule]{1, 2},_\[Rule]0},{-2}]}]; (*xMinVals are samples near local min of Abs[f[x]] where u[ x] changes sign.*) xMinVals=Take@@@Cases[xMinVals,{{__},{__}}]; If[xMinVals=!={}, soln4=Apply[ zBrent[u,##,MaxBrentIter,InitPrec,PrecGoal+GrdDigts]&, xMinVals,1,Heads\[Rule]False], soln4={} ], (*else*)soln4={};MinPos2={} ]; (*Prepare to use GoldenSecant at samples near local min of Abs[ f[x]] where u[x] does not change sign.*) If[MinPos2=!={}, xPnts=Apply[Take[xList,{#-1,#+1}]&,MinPos2,{1}, Heads\[Rule]False]; If[PrecGoal\[NotEqual]$MachinePrecision, xPnts=SetPrecision2[xPnts,PrecGoal+GrdDigts]; yPnts=Map[f,xPnts,{-1},Heads\[Rule]False], (*else*) yPnts=Apply[Take[yList,{#-1,#+1}]&,MinPos2,{1}, Heads\[Rule]False] ]; xyPnts=Replace[Transpose[{xPnts,yPnts}],List\[Rule]Sequence,{3}, Heads\[Rule]True]; soln5=Apply[ GoldenSecant[f,##,MaxSecantIter,MaxBrentIter,InitPrec, PrecGoal+GrdDigts]&,xyPnts,1,Heads\[Rule]False], (*else*)soln5={} ]; Flatten[{soln1,soln2,soln3,soln4,soln5}] ]]]; SearchForRootsF= Function[{f,xmin,xmax,MaxSecantIter,MaxBrentIter,InitialSamples,InitPrec, PrecGoal,GrdDigts}, Module[{xList,yList,ZeroPosn,CrossingPosn1,CrossingVals1,MinPos1,soln1, soln2,soln3,xPnts,yPnts,xyPnts,zeroList}, xList=InitialPoints[xmin,xmax,InitialSamples,InitPrec,PrecGoal]; yList= Developer`ToPackedArray[MakeReal[Map[f,xList,Heads\[Rule]False]]]; If[FreeQ[yList,_?NumericQ,{1}],Message[RootSearch::numb];{}, (*else we move on.*) ZeroPosn=Position[yList,_?(#\[Equal]0&),Heads\[Rule]False]; If[PrecGoal\[Equal]$MachinePrecision, soln1=Extract[xList,ZeroPosn], (*else*) soln1=SetPrecision2[Extract[xList,ZeroPosn],PrecGoal+GrdDigts]; zeroList=MakeReal[Map[f,soln1,Heads\[Rule]False]]; ZeroPosn=Position[zeroList,_?(#\[Equal]0&),Heads\[Rule]False]; soln1=Extract[soln1,ZeroPosn] ]; (*Find the zero crossings of f[x].*) CrossingPosn1= Flatten@Position[ Partition[yList,2,1],_List?(NegativeQ[First[#]*Last[#]]&),{1}, Heads\[Rule]False]; CrossingPosn1= Complement[CrossingPosn1,ZeroPosn,SameTest\[Rule]Automatic]; CrossingPosn1=Transpose[{CrossingPosn1,CrossingPosn1+1}]; (*Find samples near a local min of Abs[f[x]].*) MinPos1=Flatten@ Position[Partition[yList,3,1],_?LocalExtremaQ,{1}, Heads\[Rule]False]+1; If[CrossingPosn1=!={}, (*Apply zBrent on f[x] at the zero crossings of f[x].*) CrossingVals1= Map[Take[xList,#]&,CrossingPosn1,Heads\[Rule]False]; If[PrecGoal\[NotEqual]$MachinePrecision, CrossingVals1=SetPrecision2[CrossingVals1,PrecGoal+GrdDigts]]; soln2=Apply[zBrent[f,##,MaxBrentIter,InitPrec,PrecGoal+GrdDigts]&, CrossingVals1,1,Heads\[Rule]False], (*else*)soln2={} ]; (*Prepare to Apply GoldenSecant on f[x] where Abs[ f[x]] has a local min.*) If[MinPos1=!={}, xPnts=Map[Take[xList,{#-1,#+1}]&,MinPos1,Heads\[Rule]False]; If[PrecGoal\[NotEqual]$MachinePrecision, xPnts=SetPrecision2[xPnts,PrecGoal+GrdDigts]; yPnts=MakeReal[Map[f[#]&,xPnts,{-1},Heads\[Rule]False]], (*else*) yPnts=Map[Take[yList,{#-1,#+1}]&,MinPos1,Heads\[Rule]False] ]; xyPnts=Replace[Transpose[{xPnts,yPnts}],List\[Rule]Sequence,{3}, Heads\[Rule]True]; soln3=Apply[ GoldenSecant[f,##,MaxSecantIter,MaxBrentIter, PrecGoal+GrdDigts]&,xyPnts,1,Heads\[Rule]False], (*else*)soln3={} ]; Flatten[{soln1,soln2,soln3}] ]]]; GoldenSecant= Function[{f,xa,xb,xc,fa,fb,fc,MaxSecantIter,MaxBrentIter,PrecGoal}, Module[{x1,x2,x3,f1,f2,f3,newX,newF,dx,GiantStepCount,InfiniteStep, Passed,iter,MaxStep,tiny,RootFound, PrevDirection},({newX,newF,GiantStepCount,InfiniteStep}={xb,fb,0, False}; If[PositiveQ[ Abs[newX-xc]-Abs[newX-xa]],{x2,f2,x3,f3}={xa,fa,xc,fc},{x2,f2, x3,f3}={xc,fc,xa,fa}, Message[RootSearch::err,"GoldenSecant","1",newX]]; PrevDirection=Sign2[x2-x3]; Passed=Catch[ Do[If[f2-newF\[Equal]0, If[PositiveQ[Abs[newX-x2]-Abs[x3-newX]],newX=(newX+x2)/2, newX=(newX+x3)/2,(*Neither True or False*)Message[ RootSearch::err,"GoldenSecant","2",newX];Throw[False]]; newF=MakeReal[f[newX]], False,(*Neither True or False*)Message[RootSearch::err, "GoldenSecant","3",newX];Throw[False]]; If[PositiveQ[Abs[f2]-Abs[newF]], If[NonPositiveQ[Abs[newX-x2]-Abs[newX-x3]], {x1,f1,x2,f2}={x2,f2,newX,newF}, (*else*){x1,f1,x2,f2,x3,f3}={x3,f3,newX,newF,x2,f2}, (*Neither True or False*)Message[RootSearch::err, "GoldenSecant","4",newX];Throw[False] ], (*else*) If[NonPositiveQ[Abs[x1-x2]-Abs[newX-x2]], {x3,f3}={newX,newF}, (*else*){x3,f3,x1,f1}={x1,f1,newX,newF}, (*Neither True or False*)Message[RootSearch::err, "GoldenSecant","5",newX];Throw[False] ], (*Neither True or False*)Message[RootSearch::err, "GoldenSecant","6",newX];Throw[False] ]; tiny=Ulp2[x2,x3]; If[NonPositiveQ[Abs[x3-x2]-2*tiny], (*converged on a root*)RootFound=x2;Throw[True], False, (*Neither True or False*) Message[RootSearch::err,"GoldenSecant","7",newX] ]; newX=SecantEstimate[x1,f1,x2,f2,x3,tiny,PrecGoal, Unevaluated@PrevDirection,Unevaluated@GiantStepCount]; If[newX-x3\[Equal]0, newX=x3+Sign2[x3-x2]*Ulp2[x3,x2], False, (*Neither True or False*)Message[RootSearch::err, "GoldenSecant","8",newX] ]; newF=MakeReal[f[newX]]; Which[ (*test1*) Not@NumericQ[newF]&&(newF=!= Infinity)&&(newF=!=-Infinity), Message[RootSearch::numer,"GoldenSecant",newX]; Throw[False], (*test2*)TrueQ[newF\[Equal]0], (*found an exact root*)RootFound=newX;Throw[True], (*test3*)NegativeQ[newF*f2], RootFound={zBrent[f,x2,newX,MaxBrentIter,PrecGoal], zBrent[f,newX,x3,MaxBrentIter,PrecGoal]}; If[RootFound==={{},{}},Throw[False],(*else*)Throw[True]] ], (**Which does nothing if no test passes.**){MaxSecantIter} ]]; If[TrueQ@Passed,RootFound,{}] )]]; (********The last two arguments of SecantEstimate must be unevaluated since the values they refer to are changed by SecantEstimate.********) SecantEstimate= Function[{x1,f1,x2,f2,x3,tiny,FuncPrec,PrevDirection,GiantStepCount}, Module[{u1,u2,dx,InfiniteStep,xw=Abs[x3-x2],direction}, If[f2-f1\[Equal]0, dx=xw/2, (*else I add and then subtract x2 when computing dx to avoid \ aproblem with roundoff error. The problem surfaces when looking for a root of Abs[ Zeta[1/2+I y]] near (y=25.0109).*) dx=Abs[Plus@@{x2+f2*(x2-x1)/(f2-f1),-x2}], (*neither True or False*)Message[RootSearch::err,"SecantEstimate", "1",x2]; Throw[False] ]; direction=Sign2[x3-x2]; If[dx>100*xw, If[++GiantStepCount\[Equal]4,(*Assume we are not converging on a \ root.*)Throw[False],(*else*)dx=xw(2-GoldenRatio)], (*else*) GiantStepCount=0; If[NonNegativeQ[dx-xw], dx=xw(2-GoldenRatio), (*else*)If[direction\[NotEqual]PrevDirection, dx=Min[dx,xw(2-GoldenRatio)]], (*neither True or False*)Message[RootSearch::err,"SecantEstimate", "2",x2];Throw[False] ], (*neither True or False*)Message[RootSearch::err,"SecantEstimate", "3",x2];Throw[False] ]; dx=Max[tiny,dx]; PrevDirection=direction; If[FuncPrec\[Equal]$MachinePrecision, x2+dx*direction, (*else*)Block[{$MinPrecision=-Infinity}, x2+SetAccuracy[dx*direction, Accuracy[x2]+15]],(*neither True or False*)Message[ RootSearch::err,"SecantEstimate","4",x2]; Throw[False] ] ]]; SetAttributes[MakeReal,Listable]; MakeReal[x_Real]:=x; MakeReal[z_Complex]:=Abs[z]; MakeReal[x_?(Element[#,Reals]&)]:=x; MakeReal[z_?NumericQ]:=Abs[z]; MakeReal[x_DirectedInfinity]:=If[x===-\[Infinity],-\[Infinity],\[Infinity]]; MakeReal[_]=Indeterminate; SetPrecision2= Function[{expr,prec}, Replace[expr,{x_?(#\[Equal]0&)\[RuleDelayed]SetAccuracy[x,prec], x_?NumericQ\[RuleDelayed]SetPrecision[x,prec]},{-1}]]; If[FreeQ[Options[Precision],Verbatim[Rule][Round,_]], Precision2=(Block[{$MinPrecision=-\[Infinity]},Precision[#]]&), Precision2=(Block[{$MinPrecision=-\[Infinity]}, Precision[#,Round\[Rule]False]]&) ]; NegativeQ=(With[{y=Negative[#]}, If[y,True,False, If[InexactNumberQ[#],Negative[SetPrecision[#,15.]],y]]]&); PositiveQ=(With[{y=Positive[#]}, If[y,True,False, If[InexactNumberQ[#],Positive[SetPrecision[#,15.]],y]]]&); NonNegativeQ=(With[{y=NonNegative[#]}, If[y,True,False, If[InexactNumberQ[#],NonNegative[SetPrecision[#,15.]],y]]]&); NonPositiveQ=(With[{y=NonPositive[#]}, If[y,True,False, If[InexactNumberQ[#],NonPositive[SetPrecision[#,15.]],y]]]&); Sign2[x_?InexactNumberQ]:= If[Precision[x]\[GreaterEqual]3.0,Sign[x],Sign@SetPrecision[x,15.0]]; Sign2[DirectedInfinity[z_]]:=Sign2[z]; Sign2[x_]:=Sign[x] LocalExtremaQ=((Part[#,2]\[NotEqual]0)&&PositiveQ[First[#]/Part[#,2]-1]&& PositiveQ[Last[#]/Part[#,2]-1]&); (********* In CrossingPosition[] PosnList is a list of positions where corresponding values from xList give opposite sign for f[x].The expression CrossingPosition[..] returns is a list of positions from PosnList where corresponding values in xList give a change in sign for u[x] as well. *********) CrossingPosition=Function[{xList,PosnList,u}, Module[{lst,CrossingVals}, CrossingVals=Map[Take[xList,#]&,PosnList,Heads\[Rule]False]; lst= Thread[{Apply[{MakeReal@u[#1],MakeReal@u[#2]}&,CrossingVals,1, Heads\[Rule]False],PosnList}]; lst=Select[lst,NegativeQ[Part[#,1,1]*Part[#,1,2]]&]; Part[lst,All,-1] ]]; (********* Ulp2[x1,x2] returns the distance between x1 and another number very close to x1 and towards x2.The following cases are considered.(1) If x1 is a machine number Ulp2 returns the distance to the next machine number towards x2.(2) Otherwise Ulp2 returns the smallest positive epsilon that will give(x1+ epsilon)- x1 not equal to zero.When x1 is a machine number Ulp2 returns a machine number,otherwise Ulp2 returns a number with five more digits of accuracy than x1. *********) Ulp2=Function[{x1,x2}, Module[{eps,sgn=Sign2[x2-x1],tol}, Which[ (*test1*)sgn\[Equal]0, Message[Ulp2::fail];Throw[False], (*test2*)MachineNumberQ[x1], tol=$MachineEpsilon*Max[Abs[x1],$MinMachineNumber]/2; While[(eps=Plus@@{x1+sgn*tol,-x1})\[Equal]0.0,tol*=2]; Abs[eps], (*test3*)True, Block[{$MinPrecision=-Infinity}, tol=SetPrecision[16*10^-Accuracy[x1],20]; While[(eps=Plus@@{x1+sgn*tol,-x1};Precision[eps]\[LessEqual]0), tol*=2]; SetPrecision[Abs[eps],20] ] ]]]; InitialPoints= Function[{xmin,xmax,n,InitPrec,PrecGoal}, If[InitPrec\[Equal]$MachinePrecision, Module[{step=Re@N[xmax-xmin]/(n-1),pnts}, With[{noise=step/50,state=$RandomState,xmin2=Re@N@xmin, xmax2=Re@N@xmax}, SeedRandom[123456]; pnts=Range[xmin2+step,xmax2-step+Abs[xmax2]*$MachineEpsilon, step]; pnts=pnts+ Table[Random[Real,{-noise,noise}],Evaluate@{Length@pnts}]; pnts=Join[{xmin2,xmax2},pnts]; If[(xmin2<03)&&PoleTest[b,fb,c,fc,newB,newF], Message[RootSearch::pole,newB];Throw[False]], (*else*) biggerSample=0 ]; {b,fb}={newB,newF}; If[fb\[Equal]0,(*Found an exact root*)Throw[True]], {maxIter}]]; If[TrueQ[Passed],b,{}] )]]; BrentInterpolation=Function[{a,fa,b,fb,c,fc,tol1,FuncPrec,d,e}, (*Arguments e and d must be passed before evaluating since the \ variables they refer to are changed.*) Module[{xm=(c-b)/2,p,q,r,s}, If[NonNegativeQ[Abs[e]-tol1]&&PositiveQ[Abs[fa]-Abs[fb]], s=fb/fa; If[a-c\[Equal]0, (*take a secant step*) {p,q}={2*xm*s,1-s}, (*else attempt inverse quadratic interpolation*) {q,r}={fa/fc,fb/fc}; {p,q}={s(2*xm*q(q-r)-(b-a)(r-1)),(q-1)(r-1)(s-1)}, (*neither True or False*)Message[RootSearch::err, "BrentInterpolation","1",b]; Throw[False] ]; If[PositiveQ[p], q=-q, p=-p, (*neither True or False*)Message[RootSearch::err, "BrentInterpolation","2",b]; Throw[False] ]; If[NegativeQ[2 p-Min[Abs[e*q],3*xm*q-Abs[tol1*q]]], (*accept the interpolation*){d,e}={p/q,d}, (*else take a bisection step*){d,e}={xm,d}, (*neither True or False*)Message[RootSearch::err, "BrentInterpolation","3",b];Throw[False] ], (*else take a bisection step*){d,e}={xm,d}, (*neither True or False*)Message[RootSearch::err, "BrentInterpolation","4",b];Throw[False] ]; If[NegativeQ[Abs[d]-tol1], d=tol1*Sign2[xm], False, (*neither True or False*)Message[RootSearch::err, "BrentInterpolation","5",b]; Throw[False] ]; (*The value returned will normally be the next b.*) If[FuncPrec\[Equal]$MachinePrecision, b+d, (*else*)Block[{$MinPrecision=-Infinity}, b+SetAccuracy[d,Accuracy[b]+5]], (*neither True or False*)Message[RootSearch::err, "BrentInterpolation","6",b];Throw[False] ] ]]; RootSearchSamples:= With[{samples= Transpose[{Part[DownValues[SamplePoints],All,1,1,1], Part[DownValues[SamplePoints],All,2]}]}, Cases[samples,{_?NumericQ,_?NumericQ}]] Protect[RootSearch,InitialSamples,RootTest,MaxBrentSteps,MaxSecantSteps, PrecisionGoal,InitialPrecision]; End[]; EndPackage[];