The Museum of HP Calculators

HP Forum Archive 18

[ Return to Index | Top of Index ]

Re: New Root Seeking Algorithms
Message #1 Posted by hugh steers on 11 Nov 2007, 5:17 p.m.

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[4] 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[4] 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[1] = h * FxVal / (F(X + h) - FxVal)

NumCalls = NumCalls + 1 St[2] = St[1] * (1 + STEP_FACTOR) St[3] = St[1] * (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[4] = Sum Xnew[4] = X - St[4] Fx[4] = F(Xnew[4]) NumCalls = NumCalls + 1

print(NumCalls, Xnew[4], Fx[4])

-- 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[1] - Xnew[2]) <= XToler or abs(Fx[1]) <= 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
Message #2 Posted by hugh steers on 11 Nov 2007, 5:21 p.m.,
in response to message #1 by hugh steers

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.


[ Return to Index | Top of Index ]

Go back to the main exhibit hall