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.expfunction 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)
endfunction RootByProbingSteps(guess)
local NUM_STEPS = 3
local SLOPE_SHIFT = 5
local STEP_FACTOR = 0.15local XToler = 1e-8
local FxToler = 1e-8local X
local Xnew = {}
local St = {}
local Fx = {}X = guess
R = 2
FxVal = F(X)
NumCalls = 1h = 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
endrepeat
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
endSt[4] = Sum
Xnew[4] = X - St[4]
Fx[4] = F(Xnew[4])
NumCalls = NumCalls + 1print(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] = BuffBuff = St[I]
St[I] = St[J]
St[J] = BuffBuff = 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.