Using Newton Raphson and Gauss Seidel for nonlinear systems

Question

I need help figuring out how to incorporate Newton Raphson and Gauss Seidel methods to solve a nonlinear system of equations in Matlab. This is my logic:

1. Linearize the system by finding the derivative of the equations and put into matrix. Transfer the system so we have [matrix of partial derivatives][H] = [original system], where H is what i am solving for. H is the stepsize.
2. First plug in my intital guesses for x into the first and last matrix, and then every time I solve for H, i'll add that stepsize to the previous guess of x to create a new x. I will then plug in that updated x, find another h, and continue that process.
3. Each iteration has an error between the old and new guess. When the error is smaller than the tolerance, the values are reported.

I have before coded the Newton Raphson and Gauss Seidel methods before but separately. My main confusion comes from how to accurately set this up together and with two different variables. First, I would like to know if my logic is right, and a pseduocode for this set up would be very helpful.

I have tried to code this, and I believe that the matrices are correct, and I think that the Gauss Seidel function is correct, but I am extremely confused on how to update the guess properly. Right now, it outputs an iteration of 0 once I got rid of the errors. I was told to do a guess of [0,0] for the linear equation and [2,2] for the nonlinear equation. Basically, I am not sure how to incorporate Newton Raphson into it. I have coded the NR method before but only with 1 variable so I do not understand how to transfer this logic into this code. Here is my seperate Newton Raphson code:

Answer 1

You are on the right track, but there are a few things we could improve in your code:

1. First of all, notice that both the Newton-Raphson method and the Gauss-Seidel method are numerical methods to solve system of equations. Avoid at all costs the use of symbolic variables. Not only it is much slower, it cannot solve all numerical problems.
2. Use and abuse MATLAB syntax to deal with matrices and vectors.

Your problem is to couple the Newton-Raphson and the Gauss-Seidel solvers. I will not solve the entirety of the problem, as you said it appears your Gauss-Seidel is working.

I will also slightly change some of the notation you provided so you can fully see the MATLAB easy to use syntax.

---

Nonlinear system of equations:

The provided system of nonlinear equations can be written in vector form as:

function fx = F(x)
    fx = zeros(2,1); % pre-allocate to create a column-vector
    fx(1) = x(1)-x(2)+1;
    fx(2) = x(1)^2+x(2)^2-4;
end

This creates the column-vector function F as a function of the variables x1 (x) and x2 (y).

The Newton-Raphson algorithm requires the evaluation of the Jacobian. You need to calculate these derivatives. The functions are simple, so you can do it manually. Create a function to calculate the Jacobian as a function of the variables.

function J = jacob(x)
    J = zeros(2,2);
    J(1,1) = 1; % df/dx
    J(1,2) = -1; % df/dy
    J(2,1) = 2*x(1); % dg/dx
    J(2,2) = 2*x(2); % dgdy
end

Your logic for the Newton-Raphson is correct. Let's use it to code our Newton-Raphson:

function x = newton(fun,jacobian,x0)
    % This function uses the Newton-Raphson method to find a root for the
    % nonlinear system of equations F, with jacobian J, given an initial
    % estimate x0.

    % Here we take advantage of MATLAB syntax to easily implement our
    % routine. fun and jacobian are function handles to the equations and 
    % its derivatives. So it is computationally much more efficient to use
    % it instead of working with symbolic variables.

    % Following the steps you provided and the Newton-Raphson for one variable:
    TOL = 0.00001;
    MAXITER = 100;

    err = 1;
    iter = 0;
    x = x0;
    while err>TOL && iter<MAXITER % repeat until error is smaller than tolerance
        % step 1: linearize the system
        fx = fun(x);
        J = jacobian(x);
        % notice that the system is already in matrix form.

        % ********************************************* %
        % step 1.2: solve for stepsize H              * %
        H = -J\fx; % SOLVE LINEAR SYSTEM OF EQUATIONS * %
        % ********************************************* %

        % step 2: add stepsize to the previous guess of x
        x = x+H;

        % step 3: calculate error
        err = norm(H);
    end
end

The above function uses the Newton-Raphson method to solve for the zero of the system of nonlinear equations. Let's check if it is working on our main program:

clc; clear; close all force;

x = newton(@F,@jacob,[2;2]) % implemented Newton-Raphson
x0 = fsolve(@F,[2;2]) % Matlab built-in function

norm(x-x0) % compare both solutions

ans =

   8.9902e-10

The difference between our function and MATLAB's built-in function is very small. It appears the function works.

Notice the step 1.2 highlighted above. That is the important step. In this step, I am using the MATLAB backlash operator to solve the linear system Ax=b. The following statements have the same functionality (solve a system of linear equations):

x = A\B
x = mldivide(A,B)

Provided that you have to use the Gauss-Seidel method to solve the linear system of equations, I will leave that modifications for you to do. In step 1.2, you will change the backlash operator for the Gauss-Seidel function:

H = -J\fx;
H = gaussSeidel(-J,fx);

And properly code your gaussSeidel function (transform your code in a function):

function x = gaussSeidel(A,b)
    % Solves the system of linear equations Ax=b by the Gauss-Seidel
    % method.
    %
    % ...
    %
end

That should do the trick.