Code Continuously Runs - Newton Raphson

Hello.

My program seems to run some sometimes, but other times it continuously says it is busy and produces no result until I manually stop the program. I have tested my program with two non-linear functions with two variables and it works fine, but as soon as I call a different .m file supplied to us by our prof, it seems to get stuck sometimes.

When I stop the program I get this error:

Operation terminated by user during test (line 70)

In NRPartA (line 52)

      f2 = a(dx2);

I was just wondering if someone could explain and help me solve this problem.

Here is my code. If you guys want to see the test file please let me know. It is just a file with 32 power flow equations and a bunch of constant. It looks really hairy...:

if true

 %Define a function that will read in an m file containing our set of
%nonlinear functions, an initial set of values to use as a guess for our
%functions, a tolerance value so the program knows when to quit, and an
%iteration value for each initial value or next guess.
function F = NRPartA(a, initial, condition, step) 

%Defines the number of functions found in our called .m file.
n = length(initial);
%Initial values for the y matrix.
yinitial = a(initial);
%Creates a matrix for our starting point values.
startpoint = zeros(1,n);
%Initial Newton Raphson guess values.
dy = zeros(1,n);
%Initial set-up of the Jacobian Matrix.
Jacobian = zeros(n,n);

func = 1;

while func

%Find new value of dy.
for i = 1:n
  dy(i) = startpoint(i) - yinitial(i);
end

%Exit loop when tolerance value is met.
if abs(dy) < condition
  func = 0;
  break;
end

%Producing the Jacobian Matrix.
for i = 1:n
  %Filling in the empty matrix.
  for j = 1:n
      x1 = initial(j);

        xp = x1 + step;
        xm = x1 - step;

        dx1 = initial;
        dx1(j) = xp;
        dx2 = initial;
        dx2(j) = xm;

        f1 = a(dx1);
        f2 = a(dx2);

        guess1 = f1(i); 
        guess2 = f2(i);

        %Find each element of the Jacobian Matrix. 
        J = (guess1 - guess2)/(2 * step);
        Jacobian(i,j) = J; 
    end
end

%Sparse the Jacobian Matrix to make the program faster.
SJacobian = sparse(Jacobian);

%Start LU factorization process.
[L,U,P] = lu(SJacobian);

b = dy;
b = dy * P;
c = zeros(n,1);

%Forward substitution.
for i = 1:n
  c(i)= (b(i) - L(i, :) * c)/L(i,i); 
end

new = zeros(n,1);

%Backward substitution.
for i = n:-1:1
  new(i) = (c(i) - U(i,:) * new)/U(i,i);
  initial(i) = initial(i) + new(i);
end
 yinitial = a(initial);

end

%Output values after convergence.
F = initial;

end

end

Answer 1

[James Tursa]

I haven't looked at your algorithm in detail, but just the convergence code seems suspicious to me:

%Find new value of dy.
for i = 1:n
  dy(i) = startpoint(i) - yinitial(i);
end

%Exit loop when tolerance value is met.
if abs(dy) < condition
  func = 0;
  break;
end

I would have expected the tolerance check to be norm(dy) < condition since dy is a vector. abs(dy) doesn't make sense in this context. Also, I would expect startpoint to change each iteration but you never change it. Did you mean something like this instead?

%Find new value of dy.
for i = 1:n
  dy(i) = startpoint(i) - yinitial(i);
end
startpoint = yinitial;

%Exit loop when tolerance value is met.
if norm(dy) < condition
  func = 0;
  break;
end

[Kody Haugli]

I have tried your changes, but it still does the same thing. I actually just noticed that it switches where the error is found. Sometimes it is at f1 = a(dx1); and sometimes it is at f2 = a(dx2);

[James Tursa]

Then I would back up a step and start plotting your norm(dy) values to see if your algorithm even looks like it is converging. E.g., pre-allocate a normdy vector and save the iterative norm(dy) values in it and plot them up.

The place at which the code breaks out when you type Ctrl-C is irrelevant.