**Multiple Nonlinear Equations using the Newton-Raphson Method**

For systems of equations the Newton-Raphson method is widely used, especially for the equations arising from solution of differential equations.

Write a Taylor expansion in several variables.

The Jacobian matrix is defined as

and the Newton-Raphson method is

Since the Jacobian depends on the iterate, it must be evaluated at each iteration. Sometimes the method converges even though the Jacobian is not reevaluated at each iteration. In this formulation of the method the right-hand side gradually (hopefully!) goes to zero. There are two alternate ways to write the Newton-Raphson method.

The first one of these is useful in theoretical developments (although an LU decomposition (link) is used rather than a matrix inverse (link)). The second one is such that upon convergence the right-hand side approaches a constant that doesn't change with further iterations.

**Newton-Raphson method using MATLAB**

Next let us apply the Newton-Raphson method to the system of two nonlinear equations solved above using optimization methods. The equations to solve are

and the Jacobian is

Prepare the following script (but without the ';' at the end of each line).

% Newton Raphson solution of two nonlinear algebraic equations

% set up the iteration

error1 = 1.e8;

xx(1) = 0; % initial guesses

xx(2) = 0;

iter=0;

maxiter=30.

% begin iteration

while error1>1.e-12

iter=iter+1;

x = xx(1);

y = xx(2);

% calculate the functions

f(1) = 10*x+3*y*y-3;

f(2) = x*x-exp(y) -2.;

% calculate the Jacobian

J(1,1) = 10;

J(1,2) = 6*y;

J(2,1) = 2*x;

J(2,2) = -exp(y);

% solve the linear equations

y = -J\f';

% move the solution, xx(k+1) - xx(k), to xx(k+1)

xx = xx + y';

% calculate norms

error1=sqrt(y(1)*y(1)+y(2)*y(2))

error(iter)=sqrt(f(1)*f(1)+f(2)*f(2));

ii(iter)=iter;

if (iter > itermax)

error1 = 0.;

s=sprintf('****Did not converge within %3.0f iterations.****',itermax);

disp(s)

end

% check if error1 < 1.e-12

end

x = xx(1);

y = xx(2);

f(1) = 10*x+3*y*y-3;

f(2) = x*x-exp(y) -2.;

% print results

f

xx

iter

% plot results

semilogy(ii,error)

xlabel('iteration number')

ylabel('norm of functions')

clear ii

clear error

The first three lines initialize three things. The initial conditions for x1 and x2 are set here, and the value of error is set here. The definition of error is

and we would like it to be zero. The M-file is structured so that the calculation in the while-end
loop will continue indefinitely until error is less than
10^{-12}. Thus the value of error is initially set to
a large value, 10^{8}, so that this condition is not satisfied before beginning. (Some languages test
the condition on entry, some on exit; to be safe assume both!)

Inside the while-end loop, the values of
x_{1} and x_{2} are moved to local variables, which
are easier to use in the statements following. Then the two functions are evaluated. Next the
elements of the Jacobian are evaluated. The Newton-Raphson method uses

J^{k} ( x^{k+1}
x^{k}) = f^{k}

If we let

y = x^{k+1 } x^{k}

then the Newton-Raphson method is

J^{k} y = f^{k}

Likewise the value of x^{k+1} is

x^{k+1} = x^{k} + y

Thus the next statements compute y and
x^{k}. Since we are not sure if this problem has a
solution,

we will test 'error', which is the square root of the sum of the squares of the functions, which
should be zero. If error < 10^{-12} then the loop ends; finally it computes f and prints it. Thus we find
the solution x and the function f, which is to be zero. Note that we have to recalculate f with the
latest value of x in order to find f for the last solution.

Before we can rely on the results, we must check the M-file carefully. First remove all the ';' so that the calculations are printed out as they are made. Change the initial values from [0,0] to [1.5, 2.5]. The functions should be 30.75 and -11.9325. Run the M-file (by issuing the command >nr and quickly press the apple-period. (Alternatively, insert the command 'pause' in the code where you would like to stop and look at the results.) Look at the results on the screen for the first time through the calculation. The values of f and the Jacobian are correct.

J_{11} =10, J_{12} = 15

J_{21} = 3, J_{22} = 12.1825

We already know that MATLAB solves the linear equations correctly, but watch that you have the minus sign! Now put back in the ';' and run the problem for as many iterations as it needs. (If it takes hundreds of iterations you probably made a mistake. Stop and check your work). The results are:

xx = 1.44555236880496 2.41215834803947

error = 8.953656560885539e-13

f =

1.0e-14 *

0.17763568394003 0

The solution has been found. The change in the solution from one iteration to another was less
than 10^{-12}, and the function itself is smaller than
10^{-14} for the first element and zero to machine
accuracy for the second element. The iterative behavior is shown in the figure. It took about 20
iterations before the error started significant reductions, but once that point was reached the error
decreased rapidly.

**Figure NR2. Value of Error in the Newton-Raphson Method for Example
II.**

The successive substitution method did not work for this example.

Take Home Message: Convergence may not occur, but once the iterate is close, convergence is fast.