(10C) Mortgage Loan Interest Rate

[Albert Chan]

PHP Code:

def nsolve_i(n, r, eps=1e-8, verbal=True):
    i = r - 1.0/(r*n*n)         # guess
    prev = r
    while 1:
        fi = i/(1.-(1.+i)**-n)  # newton's method
        di = (fi - r) / (fi - n*fi*(fi-i)/(1+i)) * i
        if verbal: print '%.15g\t%+.5e' % (i, -di)
        curr = abs(di)
        if curr >= prev: break  # bad di
        i = i - di
        if curr < eps: break
        prev = curr
    return i 


Code:

>>> nsolve_i(n=36, r=86.67/2500, eps=0)
0.0124110212798083      +9.35677e-05
0.0125045889822448      -4.29793e-08
0.012504546002918       -9.07992e-15
0.012504546002909       -9.44595e-17
0.0125045460029089      +3.54223e-16
0.012504546002908864

Here is a way to avoid Newton's non-convergence due to rounding error.
Stop when adjustment size is bigger than previously.

The last bad adjustment is thrown away, and return previous iteration.
I learn this trick from Kahan's cubic equation solver paper

Edit: above for illustration only.
With eps = 1e-8, we already reached IEEE double limitation. (i = 0.012504546002909)