-1

I am trying to code a Newton's Method with back stepping code that I wrote in Matlab to Python, but am having some trouble with the Python syntax. Matlab takes about 5 iterations, but my Python code is looping up to the max iteration of 1000 and having a domain error as the back stepping mechanism does not work (tries to calculate a negative log). I have not used Python in a while so I am most likely confusing some syntax of some sort.

Here is the Matlab code that works correctly:

x = 10;                                               %defines x
f = @(x) log(x);                                        %defines objective function
df = @(x) 1/x;                                          %defines first derivative
tol = .00001;                                           %defines our tolerance level
maxit = 1000;                                           %defines maximum iteration steps
maxsteps = 200;                                         %defines maximum backsteps
for i=1:maxit                                           %starts loop
    fval = f(x);                                        %value of function at f(x)
    fjac = df(x);                                       %value of jacobian at f(x)
    fnorm = norm(fval);                                 %calculates norm value at fval
    if fnorm<tol, return, end                           %if fnorm less than tol, end
    x
    d = -(fjac\fval);                                   %forms second part of iteration rule
    d
    fnormold = inf;                                     %sets arbitrary fnormold
    for backstep=1:maxsteps
        fvalnew = f(x+d);                               %calculates f(x+d)
        fnormnew = norm(fvalnew);                       %calculates norm of fvalnew
        if fnormnew<fnorm, break, end                   %implements 1st backstepping rule
        if fnormold<fnormnew, d=2*d; break, end         %implements 2nd backstepping rule
        fnormold = fnormnew;                            %updates fnormold
        d=d/2;
    end
    x=x+d;
end
disp(x)

Here is the Python code:

from math import log

x = 10

def f(x):
    f = x* log(x)
    return f

def df(x):
    df = 1/x
    return df

tol = .00001
maxit = 1000
maxsteps = 200
maxsteps = 200

for i in range(1, maxit):
    fval = f(x)
    fjac = df(x)
    fnorm = abs(fval)
    if fnorm < tol:
        print x
    d = -(fjac/fval)
    fnormold = float('Inf')
    for backstep in range(1, maxsteps):
        fvalnew = f(x+d)
        fnormnew = abs(fvalnew)
        if fnormnew < fnorm:
            break
        if fnormold < fnormnew:
            d= 2*d
            break
        fnormold = fnormnew
        d = d/2
    x = x+d
print x
  • Could you provide the full error traceback? – jonrsharpe Jun 24 '14 at 21:53
  • You need x=10.0 not 10, to avoid integer operations; your f function should be log(x), not `x*log(x)`; your expression `d = -fjac/fval` should be `=-fval/fjac`; and the Matlab/octave code is working in complex numbers, which your python code doesn't (and so it crashes the first time it's given a negative `x+d`. You can fix that last by either just returning the real component of `log(x)` for negative x, or you can use, eg, [`cmath`](https://docs.python.org/3/library/cmath.html). – Jonathan Dursi Jun 24 '14 at 22:03
  • thanks jonathan, this worked for me. – genghistong Jun 24 '14 at 22:15

1 Answers1

0
  • 1/x in df can be 0 in most cases since in python 2.x divison over integers returns integer
  • range based for is one index too less
Luka Rahne
  • 9,760
  • 3
  • 30
  • 55