PHYS 251 - Introduction to Computer Techniques in Physics

Numerical Analysis - Algebraic Equations

Numerical Solution Using a Generalized Newton's Method

Notation: Equation numbers are placed in parentheses, (Eq. 1), at the end of the sentence preceding the mathematical presentation of the equation. They should not be read as if they are part of the sentence. When reference is made to an equation, the word "equation" will be spelled out and not abbreviated. Mathematical and Greek symbols are not in all cases precise, since they may depend on the browser used.

Mathematical Functions

It is assumed throughout this discussion that the reader is familiar with the concept of mathematical functions, such as, f(x) = x² - 2x + 36, and f(x) = sin(x). If this is not the case, then the reader should consult an introductory calculus textbook.

System of Linear Algebraic Equations

It is natural to begin a consideration of numerical methods with algebraic and transcendental systems, since the arithmetical operations performed by modern computers are those of classical algebra.

The general system of n linear algebraic equations of nth order can be written as (Eq. 1)

$a$₁₁ x₁ + a₁₂ x₂ + a₁₃ x₃ + ... + a_1n x_n = b₁,
a₂₁ x₁ + a₂₂ x₂ + a₂₃ x₃ + ... + a_2n x_n = b₂,
a₃₁ x₁ + a₃₂ x₂ + a₃₃ x₃ + ... + a_3n x_n = b₃,
...
a_n1 x₁ + a_n2 x₂ + a_n3 x₃ + ... + a_nn x_n = b_n.

If we let A be an ( n x n ) matrix of the coefficients, a_nn, x be the column vector of the variables, x_n, and b be the column vector of the right-hand side, b_n, the system can be written in compact vector notation as (Eq. 2)

$$A x = b.

For theoretical discussions Equation 2 is more convenient even though Equation 1 and 2 are equivalent.

For computational purposes, one wants to know that our system of equations, Equation 1, has one and only one solution before investing time and effort searching for something that does not exist. The basic theorem says that in order for Equation 1 to have one and only one solution, the determinant of A, written det ( A ), must be different from zero (Eq. 3)

$det\; ($A ) = | A | not= 0.

Cramer's Rule for Solving a System of Linear Algebraic Equations

If one is dealing with systems of 2, 3, or 4 linear algebraic equations, a method, known as Cramer's rule, is a reasonable technique for finding a numerical solution. Let us demonstrate the technique through an example (Eq. 4)

$x$₁ + x₂ - x₃ = 2,
x₁ - x₂ + x₃ = 0,
- x₁ + x₂ + x₃ = 0,

whose solution is x₁ = 1, x₂ = 1, and x₃ = 0.

For this example the determinant of the coefficients a_nn is (Eq. 5a)

$$| A | = | ( a₁₁, a₁₂, a₁₃ ), ( a₂₁, a₂₂, a₂₃ ), ( a₃₁, a₃₂, a₃₃ ) |
= a₁₁ ( a₂₂ a₃₃ - a₂₃ a₃₂ ) + a₁₂ ( a₂₃ a₃₁ - a₂₁ a₃₃ ) + a₁₃ ( a₂₁ a₃₂ - a₂₂ a₃₁ ) = - 4.

Thus the solution for x₁ using Cramer's rule is (Eq. 5b)

$x$₁ = | A₁ | / | A | = ( -4 ) / ( -4 ) = 1,

where | A₁ | is the determinant of A in which the column of coefficients, a₁₁ = 1, a₂₁ = 1, and a₃₁ = -1, is replaced by the column of b's, b₁ = 2, b₂ = 0, and b₃ = 0. In a similar fashion x₂ is given by (Eq. 5c)

$x$₂ = | A₂ | / | A | = ( -4 ) / ( -4 ) = 1,

where | A₂ | is the determinant of A in which the column of coefficients, a₁₂ = 1, a₂₂ = -1, and a₃₂ = 1, is replaced by the column of b's, b₁ = 2, b₂ = 0, and b₃ = 0. Finally, the solution for x₃ is given by (Eq. 5d)

$x$₃ = | A₃ | / | A | = ( 0 ) / ( -4 ) = 0,

where | A₃₂ | is the determinant of A in which the column of coefficients, a₁₃ = -1, a₂₃ = 1, and a₃₃ = 1, is replaced by the column of b's, b₁ = 2, b₂ = 0, and b₃ = 0.

Example Using Cramer's Rule for Solving a System of Linear Algebraic Equations

For an example of using the Excel function MDETERM to do this problem follow the link Cramer's Rule. As should be fairly evident, Cramer's rule becomes cumbersome when the system consists of more than about four equations.

Newton's Method for Solving One Equation in One Unknown

In addition to Cramer's rule there are other methods for the solution of systems of linear algebraic equations, such as Gaussian Elimination and LU Decomposition. They are called direct methods in as much as, barring roundoff error accumulation, they lead to a numerical solution in a prescriptive fashion with a fixed number of steps. However, in general, they do not apply to the solution of nonlinear algebraic and transcendental systems of equations. The generalized Newton's method is one that can be applied to both linear and nonlinear systems. However, before considering the generalized Newton's method for systems of equations, let us consider the classical Newton's method, often referred to as a Newton-Raphson technique, applied to one equation in one unknown.

Let us begin with the equation

$y\; =\; f(x)\; =\; 0,$

where f(x) is some function, such as the polynomial (Eq. 6)

$y\; =\; f(x)\; =\; x3+\; (31/2)\; x2-\; 2\; x\; -\; 2\; (31/2)\; =\; 0.$

The problem is one of finding the value or values of x such that y = 0, that is we are searching for the zeros of y = f(x). Figure 1 shows the function f(x) graphed for 0 < x < 2.5.

Note that the function changes sign between the values of x = 1 and x = 1.5. Thus the value of x that makes the function zero is approximately x = 1.4. The Newton's method is an iterative scheme that utilizes the tangent to the curve at some point, say x = a. Regardless of whether x = a > x_solution or x = a < x_solution, the tangent line crosses, in general, the x-axis at some point, say x₀. Since, in general, f( x₀ ) not= 0, x₀ is only a zero approximation to the root and we wish to develop a means of generating a new approximation, call it x₁.

To do so, let us note that the derivative of the function f ( x = a ), denoted by f '(a), is tangent to the curve at the point x = a. Hence, we can use the definition of the derivative as the limit of the divided difference (Eq. 7)

$(\; dy\; /dx\; )$_x=a = f '(x) = [ f(a) - f(x₀) ] / ( a - x₀ ),

in order to obtain an expression for the value of x₀ by setting f ( x₀ ) = 0 and re-arranging Equation 7, which yields (Eq. 8)

$x$₀ = a - f(a) / f '(x).

Similarly for any iteration, Equation 8 can be generalized to produce an iterative relation (Eq. 9)

$x$_n+1 = x_n - f ( x_n ) / f '( x _n ) ,

which provides an expression for calculating an improved estimate of the root in terms of the preceding estimate, the function evaluated for the preceding estimate, and the derivative of the function evaluated for the preceding estimate.

As for issues of convergence, we find that

1. Newton's method converges rapidly, but
2. it does not always converge.

An alternative route to Equation 9 is through a Taylor expansion

$f(x)\; =\; f\; (\; x$_n ) + ( x - x_n ) f '( x_n ) + [ ( x - x_n )² / 2 ] f ''( x_n ) + [ ( x - x_n )³ / 6 ] f '''( x_n ) + ...

in which we keep the first two terms only, set f(x) = 0, and re-label x as x_n+1. The result is again Equation 9.

In order to obtain an estimate of the error made in each iteration, let us subtract the i+1 iteration from the actual value of the root, x, and change the iteration index to i (Eq. 10)

$(\; x\; -\; x$_i+1 ) = ( x - x_i ) + f ( x_i ) / f '( x _i ).

The difference ( x - x_i+1 ) is the error in the i+1 iteration, call it e_i+1, while ( x - x_i ) is the error in the i iteration, call it e_i, or

$$e_i+1 = e_i + f ( x_i ) / f '( x _i ) ,

which characterizes the error at one interation to the next in terms of the function and its derivative. However, we wold like to specify the error in terms of the tangent and its rate of change, f ''(x). We can do this by using the first three terms in the Taylor expansion, which yields the result that (Eq. 11)

$$e_i+1 = - e_i² [ f ''( x_i ) / f '( x _i ) ] .

Thus for the tangent and its rate of change being reasonably constant near a root, the convergence according to Equation 11 is quadratic. But, note that if one makes an initial bad guess for the starting value x = a, so that e₀ is greater than one, then the method diverges. Hence graphing the function in order to ascertain the approximate value of the root is almost a necessity in order converge the method.

Example of Newton's Method for One Equation in One Unknown

As an example let us consider the function given in Equation 6 above, whose solution is x = 2^1/2 = 1.41421 to five decimals. The solution using the Newton's method is given in the following document Newton's Method.

Generalized Newton's Method for One Equation in One Unknown

In general, the criterion for terminating the iterations in Newton's method is that the following expression has been satisfied (Eq. 12)

$|\; x$_i+1 - x_i | < e ,

where e is a small quantity that establishes the accuracy of the calculation and is dependent on the requirements of the application. It should be obvious that depending on the nature of the function f ( x ), such as being in the neighborhood of a local maximum or minimum where the tangent becomes nearly parallel to the x axis, the convergence can require a large number of iterations to achieve any reasonable accuracy.

For those functions requiring a large number of iterations, one can attempt to construct a line rather than the tangent at x₀ = a for which x_solution < x₊ < x₁. Thus x₊ is better approximation to the root than is x₁. The equation of such a line can be written as (Eq. 13)

$y\; =\; f\; (\; x$₀ ) + t f '( x₀ ) ( x - x₀ ) ,

where t is a number close to unity. Setting y = 0 and x = x₁ = x₊, Equation 13 yields

$x$₁ = x₀ - [ f ( x₀ ) / t f '( x₀ ) ] ,

which can be generalized to produce a recursion relation for x. However, for notational purpose, let us set w = 1 / t with 0 < w < 2 so that the recursion relation becomes (Eq. 14)

$x$_i+1 = x_i - w [ f ( x_i ) / f '( x_i ) ] .

Equation 14 is known as a generalized Newton's method and w is called an over-relaxation parameter. Of course we recover the usual Newton's method when we set w = 1.

Generalized Newton's Method for a System of Two Algebraic Equations in Two Unknowns

Suppose one wishes to solve a system of two equations in two unknowns, such as (Eq. 15)

$y$₁ = f₁( x₁, x₂ ) = 0 ,
y₂ = f₂( x₁, x₂ ) = 0 .

A generalization of Equation 14 for iterative relations for x₁ and x₂ in Equation 15 will require us to change notation, since we are identifying the various independent variables x with a subscript. Thus, in order to identify the number of iterations, we will use the notation, x_i^(k), for the kth iteration of the ith variable. Also, since we are dealing with functions that are dependent on more than one independent variable, we must take the partial derivative with respect to each of the variables rather than the derivative df/dx = f '(x). The symbol we will use for the partial derivative is df₁( x₁, x₂ )/dx₁. Using this notation, the generalization of the generalized Newton's method for two functions in two unknowns is (Eq. 16)

$x$₁^(k+1) = x₁^(k) - w [ f₁ ( x₁^(k), x₂^(k) ) / df₁( x₁^(k), x₂^(k) )/dx₁ ] ,
x₂^(k+1) = x₂^(k) - w [ f₂ ( x₁^(k+1), x₂^(k) ) / df₂( x₁^(k+1), x₂^(k) )/dx₂ ].

It is reasonably obvious as to what would be the generalization of Equation 16 to a system of three algebraic equations in three unknowns. The point to note is that the new estimate for a variable is used in subsequent calculations involving that variable as is evident in Equation 16. That is to say, the k+1st estimate of x₃ would use the k+1st estimates for x₁ and x₂, but the kth estimate for x₃.

Example of Generalized Newton's Method for Two Equations in Two Unknowns

As an example let us consider the functions given below