%
%   NLDEMO1_LESSON6.M  Solve nonlinear demo using the linearized iteration method
%
%   This file computes the solution for a 3rd order nonlinear system via
%   the linearized iteration (successive substitution) method.  In this method, 
%   the coefficient matrix is a function of the solution vector.  However, with 
%   a given guess for x, A(x) can be determined, and the solution of a linear 
%   system, A(x) * x = b, can be obtained.  This process is continued until 
%   convergence (or we hit the max number of iterations).
%
%   Note that NLDEMO2 and NLDEMO3... solve this same problem using different methods.
%
%   File prepared by J. R. White, UMass-Lowell (last update:  Nov. 2017)
%
      
      clear all;   close all;   nfig = 0;
%
%   set reasonable initial guess 
      x2 = input(' input a guess for x2: ');  
      if x2 == 3, x2 = 3.001; end   % don't allow x2 = 3 as guess
      x1 = (-7+2*x2*x2)/3;          % satisfies eqn. 1
      x3 = -5/(3-x2)  ;             % satisfies eqn. 3
      xold = [x1 x2 x3]';     
%   
%   select eqn form (there were two obvious forms for A(x))
      form = menu('Select Equation Form', ...
                  '        Form #1        ', ...
                  '        Form #2        ');
%
%   start iteration loop
      itmax = 30;   it = 0;   tol = 1e-5;   emax = 1;   n = length(xold);
      fprintf(1,'\n Intermediate edit for NLDemo1 (Form #%1i) \n',form);
      while emax > tol && it <= itmax
        it = it+1;  
        x = xold;
%   
%   set up A and b matrices given current guess and solve for new solution vector
        switch form
          case 1
            A = [  3    -2*x(2)     0;          % coefficient matrix (form #1) 
                 -x(2)     3       -2;
                   0      -x(3)     3];
          case 2
            A = [  3    -2*x(2)     0;          % coefficient matrix (form #2)
                   0     3-x(1)    -2;
                   0       0       3-x(2)];
        end
        b = [-7 9 -5]';                         % RHS vector (constant for this case)
        xnew = A\b;                             % find solution vector
%
%   calc & edit error (intermediate results)
        emax = max(abs((xnew-xold)./xnew));
        fprintf(1,'  it = %3d      max error = %8.3e \n',it,emax);
        fprintf(1,'      xnew         xold     \n');
        for j = 1:n
          fprintf(1,' %10.5f   %10.5f    \n',xnew(j),xold(j));
        end
%
        xold = xnew;  % use current estimate as guess for next iteration
      end
%
%   print final max relative error and iteration count
      fprintf(1,'\n Number of iterations to convergence = %3d\n',it);
      fprintf(1,' Max relative error at convergence =  %8.3e\n',emax);
      if it >= itmax 
        fprintf(1,' ***** WARNING  --  Hit max number of iterations!!!  *****\n');
      end
%
%   end of file