mysolve help

I need to write a function MySolve that will solve the nonlinear system of equations in the form

0=f(x)

where I have f:R^n -> R^n is an arbitrary function with n-dimensions input and n-dimensional output. The function should implement the newton iteration.

so far i have

 function [x,converged]=MySolve(f,xold,tol,maxit)

 % f:function defining the non-linear equations. f should accept a single
 % arguement x and return a single value y
 % x0: initial guess for Newtons iteration (n-dimentional vector)
 %tol stop the iteration and set converged to true denoted by 1 if
 %correction ek and residual rk both have norms less than tol.
 % x: solution of nonlinear system (if newton iteration has converged)
 % converged: flag indicating of newton iteration has converged

 %maxit maximum number of iterations to be tried. 

 x=xold; 
 h=1e-10;

 % run a loop from 1 to maxit

         for k=0:maxit   

         %Need to call in MyJacobian    

         J=Myjacobian(f,x,h);

         % Newton iteration

         x= xold-(J\f(xold));

            if(max(abs(x-xold)))<tol && (max(abs(x-xold)))<tol

            converged=1; % Newtons iteration has converged

            else

            converged=0; %Newtons iteration hasn't converged

            end

 xold=x;

        end

 end

my jacobian is

 function df=Myjacobian(f,x,h)
 % f: function to be differentiated
 %x: point where jacobian is taken
 % h: parameter for finite differences
 % outputs df: m*n matrix, jacobian of f in x.

 n=length(x); % defines number of rows 
 fx=f(x);    
m=length(fx); % defines number of colums
df=zeros(n,m); % matrix of zeros 

    for i=1:n; 

        % runs a loop that takes two values x1 and x2 and places it
        % into my empty matrix df the process then produces 2 matircies of 
        % df1 and df2

       x(i)=x(i)+h;
        x1= f(x);
        x(i)=x(i)-(2*h);
        x2= f(x);
    df(:,i)=(x1-x2)/(2*h);
    x(i)=x(i)+h;

    end

end

so yeah a start but a few things missing if anyone could help it would be greatly appreciated thanks

Answer 1

[Geoff]

Does it work? What specific answer do you require? I can offer a couple of things just from a quick glance...

Your comment says you iterate over 1:maxit, and the code says 0:maxit.

When you set the converged flag, you still iterate up to maxit because you're not checking it. You should do this:

converged = 1;
break;

Also, it's good practise to set converged to zero BEFORE the loop, and DON'T set it to zero again if not converged. Only change it if you DO converge.

In that respect, you might choose to use a while-loop instead of a for-loop. There are some who argue that the prolific use of break is bad coding practise. But in practise it's quite valid, often more efficient, and incredibly useful. Nevertheless, here is the alternative (drawback is you have to remember to increment 'k', so this becomes less desirable if you practise the use of the continue keyword)

k = 1;
converged = 0;

while ~converged && k <= maxit
  % etc
  k = k + 1;
end

I notice your if-statement tests exactly the same thing twice. It's no wonder you didn't notice because it's not very legible. Let me lay out the function more clearly using your exact code, with my suggestions:

converged = 0;

for k = 1:maxit

    J = Myjacobian(f,x,h);
    x = xold - (J \ f(xold));

    if max(abs(x-xold)) < tol
      converged = 1;
      break;
    end

end

Notice the use of a bit of whitespace. I also removed the parentheses enclosing the max() call... There is no need for them and it makes your code less clear.

Just a point of note.... I'm not going to think too hard here, but shouldn't you be storing xold again after your tests? ie at the end of the loop:

xold = x;

Then where you initialise converged, you should do:

xold = NaN;

That takes advantage of a handy MatLab side-effect to avoid extra logic for the first iteration, or duplicating code.

[ edit: oh, yes... and don't pass xold as a function parameter - that's meaningless ]

That's my lot. Think about consistent and practical use of whitespace...

[Joe Howes]

thanks will take a look at it and in term of answers i am looking for i am using rosenbrocks banana to find the minimum of f(x)=((1-x1)^2)+100(x2-x1^2)^2

[Joe Howes]

oh and just figured where one of mymistakes were thank you one was meant to be the correction norm max(abs(x-xold)); and the other was meant to be the residual max(abs(f(xold))); where part of the problem was

[Geoff]

Well, that's probably the subject of a new question, and you still haven't said what that actual question might be. It seems to be that in this question, you've asked for comments on your code with absolutely no mention of that function. I think you'll find my suggestions useful in terms of making your code work correctly and also making it easier to read.

[Geoff]

Hehe, yes, you would have seen that more easily if you hadn't had so many parentheses and so little whitespace =)

[Joe Howes]

ok the code i am using to check that it is working is

%% Test Newton iteration using Rosenbrock's banana

% we find the minimum of

obj=@(x)((1-x(1))^2+100*(x(2)-x(1)^2)^2);

% by solving grad obj(x)=0 with a Newton iteration:

h=1e-6;

f=@(x)(Myjacobian(obj,x,h)'); %f is grad obj

x0=[0;0]; %x0 is the initial guess

[x,conv]=MySolve(f,x0);

% the solution should be [1;1] after 5 (or so) iterations

x

which should give me an answer that converges