gausseid.m

% Illustration of Gauss-Seidel iteration with relaxation

% The goal is to determine a value of x1 and x2 that simultaneously
% satisfy the following equations...
%
%   f1(x1,x2) = x1^2 + x2^2 + exp(x1) - 7.7183*x1 = 0
%   f2(x1,x2) = x2 + exp(x2) + x1^3 - 10.389      = 0

% we will take an x1 out of f1(x1,x2) and an x2 out of f2(x1,x2)
% thus we will have
%
%   x1 = g1(x1,x2) = (x1^2 + x2^2 + exp(x1))/7.7183 and
%   x2 = g2(x1,x2) = 10.389 - exp(x2) - x1^3

% in file 'g1.m', we would have:
%
%   function y = g1(x1,x2)
%   y = (x1^2 + x2^2 + exp(x1))/7.7183;

% in file 'g2.m', we would have
%
%   function y = g2(x1,x2)
%   y = 10.389 - exp(x2) - x1^3;


fprintf('\nSolution of Non-Linear Equations using Gauss-Seidel iteration\n');

      x1 = input('Enter x1 estimate:                   ');
      x2 = input('Enter x2 extimate:                   ');
       w = input('Enter relaxation factor:             ');
maxerror = input('Enter max relative error (percent):  ');
   maxit = input('Enter max # of iterations:           ');

  count = 0;
  error = 1;

  fprintf('\niteration      x1           x2      error\n');
  while (error > maxerror) & (count < maxit)
    count = count + 1;

    x1new = g1(x1,x2);
    error1 = abs((x1new - x1)/x1new*100); % determine error
    x1 = x1 + w * (x1new - x1);
    
    x2new = g2(x1,x2);
    error2 = abs((x2new - x2)/x2new*100); % Note: this could all be vectorized!
    x2 = x2 + w * (x2new - x2);

    error = max(error1,error2);           % find maximum error
  
   fprintf('%6g    %10g %10g %10g\n',count,x1,x2,error);

  end