 Re: New Root Seeking Algorithms - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: Re: New Root Seeking Algorithms (/thread-128291.html) Re: New Root Seeking Algorithms - hugh steers - 11-11-2007 hello namir, this is an interesting idea. what it seems you are doing here is effectively creating a polynomial approximation in the neighborhood of the root and using this for a convergent guess. presumably, 3 points would lead to a parabolic approximation. it occurred to me that your approach ought to be similar to the ridders method which also uses parabolic approximation, although this is not the case. running the ridders method gives far more function evals than you have here. also you are forced to bracket the root for ridders rather than supplying a simple guess. this leads to some points. i know that, at this stage, you are experimenting with the algorithm, but it will need some damage control. for example, if St is not better than any of the previous values, it will be discarded and you will be in a repeating loop. i was wondering if this could be protected against by, for example, always ensuring that the 3 values kept bracket the root when possible; ie if they ever bracket the root, dont ever through away a value that results in a non-bracketing. once you have a root bracket, you can always do a bisection or possible secant when your St is weak. sometimes, like in your last example, the guess does not lead to the closest root. this could mean, for example, that a systematic campaign of starting values, may miss a root. however, this is true of newton too. i translated your VBA into LUA. here is a transcript of the converted program. you will see that it is very similar to your original. i have delberately tried to keep it so. abs = math.abs exp = math.exp function F(X) return exp(X) - 3 * X ^ 2 --return exp(X) - exp(-3) --return 0.005 * (X+5) * (X+3) * (X+1) * (X-5) * (X-3) * (X-1) end function RootByProbingSteps(guess) local NUM_STEPS = 3 local SLOPE_SHIFT = 5 local STEP_FACTOR = 0.15 local XToler = 1e-8 local FxToler = 1e-8 local X local Xnew = {} local St = {} local Fx = {} X = guess R = 2 FxVal = F(X) NumCalls = 1 h = 0.01 * (1 + abs(X)) St = h * FxVal / (F(X + h) - FxVal) NumCalls = NumCalls + 1 St = St * (1 + STEP_FACTOR) St = St * (1 - STEP_FACTOR) for I = 1,3 do Xnew[I] = X - St[I] Fx[I] = F(Xnew[I]) NumCalls = NumCalls + 1 end repeat Sum = 0 for I = 1, NUM_STEPS do Prod = St[I] for J = 1, NUM_STEPS do if I != J then Prod = Prod * (0 - Fx[J]) / (Fx[I] - Fx[J]) end end Sum = Sum + Prod end St = Sum Xnew = X - St Fx = F(Xnew) NumCalls = NumCalls + 1 print(NumCalls, Xnew, Fx) -- remove worst fx() value for I = 1, NUM_STEPS do for J = I + 1, NUM_STEPS + 1 do if abs(Fx[I]) > abs(Fx[J]) then Buff = Fx[I] Fx[I] = Fx[J] Fx[J] = Buff Buff = St[I] St[I] = St[J] St[J] = Buff Buff = Xnew[I] Xnew[I] = Xnew[J] Xnew[J] = Buff end end end R = R + 1 until abs(Xnew - Xnew) <= XToler or abs(Fx) <= FxToler end i get the same results as you do except for the second example, where i get different numbers of NumCalls. More unfortunately. > RootByProbingSteps(-2) 6 -2.987219019322218338792071 0.00064041138381297847374604 7 -2.999613965317089589248981 0.1922324533141019322479e-4 8 -2.999999566911772811857812 0.0002156219784549971291e-4 9 -2.999999999999290070803826 0.03534529342597e-12 > RootByProbingSteps(-1) 6 -2.77277612899414810051559 0.01270121980745993109340363 7 -2.950229467075658006405727 0.00254062871917551029598691 8 -2.997235451795202188353435 0.0001378291804025967635871 9 -2.99999028549716645878965 0.0048365896598271835341e-4 10 -2.999999999560278729216809 0.002189243297610219e-8 > RootByProbingSteps(0) 6 -2.20194163657231040709631 0.06080115925608043659108856 7 -2.636856692870105576230211 0.02179886448380611134617641 8 -2.911035175923267599061736 0.00463229871509569129373929 9 -2.993687080339088281855584 0.00031529593458869318827063 10 -2.999939414762663026366751 0.0301645272870873645619e-4 11 -2.999999988926047158788972 0.055133965026063744e-8 > RootByProbingSteps(1) 6 -1.3813889457765245250755 0.2014422982797888952523052 7 -1.984256999208900514530413 0.08769565758395040579026306 8 -2.527063422843063025775631 0.03010622058300368778323235 9 -2.874397356686962420726117 0.00666308019176851773462704 10 -2.98646226635607850631352 0.00067858697161281661401285 11 -2.999770146137248562266027 0.1144506527526535822222e-4 12 -2.999999874202242755785782 0.626310193438372435e-8 > RootByProbingSteps(2) 6 -0.44806273902053333180207 0.5890775326673590685193263 7 -1.126892727991626695537758 0.2742515013318770758184717 8 -1.827612065763559220999504 0.1110100141795754388565973 9 -2.44017324764718846757097 0.0373586839896545465477001 10 -2.821342460717671200475147 0.00973890939136019973572929 11 -2.975445704868764523985792 0.00123761861602227210536844 12 -2.999323739049539909131346 0.3368043727374816436782e-4 13 -2.999999060537709997516538 0.0004677309523215980816e-4 14 -2.999999999994832780936982 0.25726068876829e-12 compare this to the first example, where both the answers, NumCalls and the error are the same. > RootByProbingSteps(7) 6 4.94935387865757248624523 67.59546577772976358203062 7 4.490838244780918934086408 28.69329836270818012735 8 4.100808982436306600194201 9.93921864919306164613123 9 3.849532676496498443596026 2.51440191702456352024564 10 3.751984513614966147802354 0.37338608067801946048483 11 3.733822283783513707463962 0.01443585375480962405536 12 3.733080996527429246405105 0.3819525654754575604e-4 13 3.733079028669344302216203 0.070901869376517e-8 > RootByProbingSteps(6) 6 4.285617914368213914777376 17.54786001602111130461905 7 3.991563636848715502088463 6.34173689463161285061912 8 3.805326804336186839048452 1.49839801545011555074212 9 3.740479979446648712053394 0.14462975514681120562327 10 3.733221647178180174021681 0.00276847076865353853075 11 3.733079124549602732329736 0.0186166474107468597e-4 12 3.733079028632949295861292 0.000262208555262e-8 > RootByProbingSteps(5) 6 3.850545108451683894747413 2.53859363395771158802991 7 3.757822511447820942221247 0.49131793540958246580583 8 3.734418165868706300372532 0.02602365624880021132889 9 3.733083682007478130284536 0.903184947158656975e-4 10 3.733079028835428155946392 0.39325672716101e-8 if i get time, i'll try some more tests. Re: New Root Seeking Algorithms - hugh steers - 11-11-2007 aha! sorry i was solving the wrong equation for case (2) should have been. function F(X) return exp(-X) - exp(-3) end i can now reproduce your results exactly.