Re: New Root Seeking Algorithms



#2

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.


#3

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.


Possibly Related Threads...
Thread Author Replies Views Last Post
  [HP Prime] Using n-root symbol and exponent problem uklo 7 1,792 11-11-2013, 01:39 AM
Last Post: Alberto Candel
  OT--TI-36X Algorithms Matt Agajanian 48 6,814 09-01-2013, 08:13 PM
Last Post: robert rozee
  Cubic root (-8) = 2 ? Gilles Carpentier 37 6,153 08-12-2013, 10:26 PM
Last Post: jep2276
  HELP WANTED ON ALGORITHMS Joerg Woerner 19 3,015 04-27-2013, 12:56 PM
Last Post: Eric Smith
  Square Root Simplifier for HP39gII Mic 4 1,299 03-11-2013, 08:25 AM
Last Post: Eddie W. Shore
  Seeking DIY4X update David Griffith 10 1,759 11-15-2012, 12:41 AM
Last Post: Eric Smith
  Cube root on standard calculator Thomas Klemm 22 4,040 11-09-2012, 06:11 AM
Last Post: Pierre
  Seeking new resources for TI NSpire Programming Namir 6 1,497 07-29-2012, 03:57 AM
Last Post: fhub
  ROOT bug? HP 48S/48G Eddie W. Shore 8 2,812 07-13-2012, 07:05 PM
Last Post: Eddie W. Shore
  OT: primitive mult/div algorithms Egan Ford 9 1,665 05-27-2012, 09:27 PM
Last Post: Egan Ford

Forum Jump: