06-23-2022, 04:36 PM
(06-11-2022 05:34 PM)Albert Chan Wrote: [ -> ]Rational Halley's modified slope = f' * (1 - (f*f'')/(f')^2/2) = f' - (f''/2)*(f/f')
With slope adjustment, (f''/2)*(f/f'), getting bigger in size relative to f',
it may turn slope to the wrong direction! Or, worse ... slope of 0!
Below example, with i0 = pmt/fv = 1e-9, 1st iteration get worse (turned negative!)
True rate = 5.32155425915678; Rate iteration should not go backward.
lua> n,pv,pmt,fv = 10, -100, 10, 1e10
lua> loan_rate2(n,pv,pmt,fv)() + 0
-0.08736585216654838
Checking f and its derivative at i=i0, Halley's slope had wrong sign.
lua> f0, f1, f2 = 999999995.5, -4500000102.007542, 143515080310.56577
lua> f1 - (f2/2)*(f0/f1) -- bad slope
11446119499.270123
Technically, with guess rate from the edge, this should not happen.
The reason is a badly estimated f'', over-estimated almost 9 times!
With more precision, this is true f0, f1, f2, at i=i0
lua> f0, f1, f2 = 999999995.5, -4500000038.5, 16499999810.25
lua> f1 - (f2/2)*(f0/f1) -- true halley's slope
-2666666750.1851845
lua> 1e-9 - f0/_ -- this should be i1
0.37499998756770886
Instead of special cased at *exactly* i=0, we should widen it.
Below code treat eps=i*i, such that 1+eps==1, as special.
From above examples, slope is less affected by catastrophic cancellations.
But, I fix it anyway for Newton's method too (= Halley's method code, with f''=0)
For naming consistency, Newton's = iter_i(), Halley's = iter_i2()
Code:
function iter_i2(n, pv, pmt, fv, i)
pmt, fv = pmt or -1, fv or 0
i = i or edge_i(n, pv, pmt, fv)
return function()
local eps, f = i*i
if 1+eps==1 then
local a, b = (pv+fv)/n, pv-fv
local fp, fpp = (a+b)*0.5, (n*n-1)*a/6
f, fp, fpp = (a+pmt)+fp*i, fp+fpp*i, fpp*(1-1.5*i)
eps = -f*fp / (fp*fp - .5*f*fpp)
else
local x = i/expm1(n*log1p(i))
local k, y = (pv+fv)*x, n*x-1
local num, den = y+(n-1)*i, i+eps
f = k + pv*i + pmt
x = k*num - pv*den -- f' * -den
y = k*(num+i+1)*(num+y) -- f'' * den^2
eps = f*x*den / (x*x-.5*f*y)
end
i = i + eps -- Halley's method
return i, eps, f
end
end
function iter_i(n, pv, pmt, fv, i)
pmt, fv = pmt or -1, fv or 0
i = i or edge_i(n, pv, pmt, fv)
return function()
local eps, f = i*i
if 1+eps==1 then
local a, b = (pv+fv)/n, pv-fv
local fp, fpp = (a+b)*0.5, (n*n-1)*a/6
f, fp = (a+pmt)+fp*i, fp+fpp*i
eps = -f / fp
else
local x = i/expm1(n*log1p(i))
local k, y = (pv+fv)*x, n*x-1
local num, den = y+(n-1)*i, i+eps
f = k + pv*i + pmt
eps = f / (k*num/den - pv)
end
i = i + eps -- Newton's method
return i, eps, f
end
end
function edge_i(n,pv,pmt,fv)
local a, b = pmt/fv, -pmt/pv
if abs(a)>abs(b) and b>-1 or b~=b then a,b = b,a end
return a, b -- a = small edge or NaN
end
lua> n,pv,pmt,fv = 10, -100, 10, 1e10
lua> iter_i(n,pv,pmt,fv)() + 0 -- Newton's method
0.22222222032098768
lua> iter_i2(n,pv,pmt,fv)() + 0 -- Halley's method converge faster
0.3749999875677088