martes, 27 de julio de 2010

Gaussian Elimination

miércoles, 21 de julio de 2010

LU DECOMPOSITION

A procedure for decomposing an N×N matrix A into a product of a lower triangular matrix L and an upper triangular matrix U,
 LU=A.

LU decomposition is implemented in Mathematica as LUDecomposition[m].
Written explicitly for a 3×3 matrix, the decomposition is

 [l_(11) 0 0; l_(21) l_(22) 0; l_(31) l_(32) l_(33)][u_(11) u_(12) u_(13); 0 u_(22) u_(23); 0 0 u_(33)]=[a_(11) a_(12) a_(13); a_(21) a_(22) a_(23); a_(31) a_(32) a_(33)]
 [l_(11)u_(11) l_(11)u_(12) l_(11)u_(13); l_(21)u_(11) l_(21)u_(12)+l_(22)u_(22) l_(21)u_(13)+l_(22)u_(23); l_(31)u_(11) l_(31)u_(12)+l_(32)u_(22) l_(31)u_(13)+l_(32)u_(23)+l_(33)u_(33)]=[a_(11) a_(12) a_(13); a_(21) a_(22) a_(23); a_(31) a_(32) a_(33)].

This gives three types of equations

 i<j    l_(i1)u_(1j)+l_(i2)u_(2j)+...+l_(ii)u_(ij)=a_(ij)

 i=j    l_(i1)u_(1j)+l_(i2)u_(2j)+...+l_(ii)u_(jj)=a_(ij)

 i>j    l_(i1)u_(1j)+l_(i2)u_(2j)+...+l_(ij)u_(jj)=a_(ij).

This gives N^2 equations for N^2+N unknowns (the decomposition is not unique), and can be solved using Crout's method. To solve the matrix equation

 Ax=(LU)x=L(Ux)=b,

first solve Ly=b for y. This can be done by forward substitution

for i=2, ..., N. Then solve Ux=y for x. This can be done by back substitution


for i=N-1, ..., 1




martes, 20 de julio de 2010

SYSTEMS OF LINEAR EQUATIONS: SOLVING BY SUBSTITUTION

The method of solving "by substitution" works by solving one of the equations (you choose which one) for one of the variables (you choose which one), and then plugging this back into the other equation, "substituting" for the chosen variable and solving for the other. Then you back-solve for the first variable. Here is how it works. (I'll use the same systems as were in a previous page.) Solve the following system by substitution.

    The idea here is to solve one of the equations for one of the variables, and plug this into the other equation. It does not matter which equation or which variable you pick. There is no right or wrong choice; the answer will be the same, regardless. But — some choices may be better than others.

    For instance, in this case, can you see that it would probably be simplest to solve the second equation for "y =", since there is already a y floating around loose in the middle there? I could solve the first equation for either variable, but I'd get fractions, and solving the second equation forx would also give me fractions. It wouldn't be "wrong" to make a different choice, but it would probably be more difficult. Being lazy, I'll solve the second equation for y:


    Now I'll plug this in ("substitute it") for "y" in the first equation, and solve for x:

    Now I can plug this x-value back into either equation, and solve for y. But since I already have an expression for "y =", it will be simplest to just plug into this:

        Then the solution is (xy) = (5, 4).


    Warning: If I had substituted my "–4x + 24" expression into the same equation as I'd used to solve for "y =", I would have gotten a true, but useless, statement:

    Twenty-four does equal twenty-four, but who cares? So when using substitution, make sure you substitute into the other equation, or you'll just be wasting your time.

GAUSSIAN ELIMINATION

Solving three-variable, three-equation linear systems is more difficult, at least initially, than solving the two-variable systems, because the computations involved are more messy. You will need to be very neat in your working, and you should plan to use lots of scratch paper. The method for solving these systems is an extension of the two-variable solving-by-addition method, so make sure you know this method well and can use it consistently correctly.

Though the method of solution is based on addition/elimination, trying to do actual addition tends to get very messy, so there is a systematized method for solving the three-or-more-variables systems. This method is called "Gaussian elimination" (with the equations ending up in what is called "row-echelon form").
Let's start simple, and work our way up to messier examples

Solve the following system of equations

It's fairly easy to see how to proceed in this case. I'll just back-substitute the z-value from the third equation into the second equation, solve the result for y, and then plug z and y into the first equation and solve the result for x.

    Then the solution is (x, y, z) = (–1, 2, 3).
The reason this system was easy to solve is that the system was "triangular"; this refers to the equations having the form of a triangle, because of the lower equations containing only the later variables.



Gaussian elimination

GAUSS-JORDAN ELIMINATION

Gauss-Jordan Elimination is a variant of Gaussian Elimination. Again, we are transforming the coefficient matrix into another matrix that is much easier to solve, and the system represented by the new augmented matrix has the same solution set as the original system of linear equations. In Gauss-Jordan Elimination, the goal is to transform the coefficient matrix into a diagonal matrix, and the zeros are introduced into the matrix one column at a time. We work to eliminate the elements both above and below the diagonal element of a given column in one pass through the matrix.
The general procedure for Gauss-Jordan Elimination can be summarized in the following steps:
  1. Write the augmented matrix for the system of linear equations.
  2. Use elementary row operations on the augmented matrix [A|b] to transform A into diagonal form. If a zero is located on the diagonal, switch the rows until a nonzero is in that place. If you are unable to do so, stop; the system has either infinite or no solutions.
  3. By dividing the diagonal element and the right-hand-side element in each row by the diagonal element in that row, make each diagonal element equal to one.
Since the matrix is representing the coefficients of the given variables in the system, the augmentation now represents the values of each of those variables. The solution to the system can now be found by inspection and no additional work is required. Consider the following example:

It is now obvious, by inspection, that the solution to this linear system is x=3, y=1, and z=2. Again, by solution, it is meant the x, y, and z required to satisfy all the equations simultaneously.


Elimination Gauss-Jordan

THE SECANT METHOD

A potential problem in implementing the Newton-Raphson method is the evaluation of the derivative. Although this is not inconvenient for polynomials and many other functions, there are certain functions whose derivatives may be extremely difficult or inconvenient to evaluate. For these cases, the derivative can be approximated by a backward finite divided difference, as in
This approximation can be substituted to yield the following it equation:
Equation is the formula for the secant method. Notice that the approach requires two
Initial estimates of x. However, because f(x) is not required to change signs between
the estimates, it is not classified as a bracketing method.




THE NEWTON·RAPHSON METHOD

Perhaps the most widely used of all root-locating formulas is the Newton-Raphson equation.
If the initial guess at the root is Xi. a tangent can be extended from the point [Xi,f(Xi)]. The point where this tangent crosses the x  axis usually represents an improved estimate of the root

The Newton-Raphson method can be derived on the basis of this geometrical interpretation (an alternative method based on the Taylor series). As in the first derivative al x is equivalent to the slope:

This can be rearranged to yield



CONVERGENCE

Notice that the true percent relative error for each iteration of Example is roughly proportional,(by a factor of about 0.5 to 0.6) to the error from the previous iteration. This property, called linear convergence, is characteristic of fixed-point iteration. Aside from the "rate" of convergence, we must comment at this point about the "possibility" of convergence. The concepts of convergence and divergence can be depicted graphically. Recall that in , we graphed a function to visualize its structure and behavior. Such an approach is employed for the function

f(x) =(e^(-x) – x).

 An alternative graphical approach is to separate the equation into two component parts, as in

 f1(x) =f2(x)

Then the two equations

 Y1 = f1(x)                and

y2= f2(x)

can be plotted separately .

The x values corresponding to the intersections of these functions represent the roots of f(x) = O.

The two-curve method can now be used to illustrate the convergence and divergence of tixed-point iteration. First, can be fe-expressed as a pair of equations y1 = x and Y2 = g(x). These two equations can then be plotted separately. As was the case with and the r00tS of f(x) = 0 cOrrespond to the abscissa value at the intersection of the two curves. The function Y1 = x and four different shapes for Y2 = g(x) are plotted in For the first case, the initial guess of Xo is used to determine the corresponding point on the Y2 curve [xo. g(xo)].

 The point (X1, X1) is located by moving left horizontally to the Y I curve. These movements are equivalent to the first iteration in the fixed-point method:

 X I = g(xo)

 Thus, in both the equation and in the plot, a starting value of xo is used to-obtain an estimate of X1., The next iteration consists of moving to [Xl, g(xl)] and then to (X2, X2), This iteration is equivalent to the equation X2 = g(x1) The solution is convergent because the estimates of x move closer root with each iteration, where the iterations diverge from the root. Notice that Convergence seems to occur only when the absolute value of the slope of Y2=g(x) is less than the slope of Y1 =x, that is, when g’(x) < 1 . Box provides a theoretical derivation of this result.

X2 =g(X1)

SIMPLE FIXED•POINT ITERATION

As mentioned above, open methods employ a formula to predict the root. Such a formula can be developed for simple fixed-poil1t iteration (or, as it is also called, one-point iteration or successive substitution) by rearranging the function f(x) = 0 so that x is or side of the equation: x=g(x) This transformation can be accomplished either by algebraic manipulation or by simply adding x to both sides of the original equation.

For example, x^2-2x+3=0 Can be simply manipulated to yield .

 x=(x^2+3)/2

 Whereas sin x=0 could be put into the form of equation by adding x to both sides to yield X=sin x +x The utility of Equation is that it provides a formula to predict a new value of x as a function of an old value of x.

 Thus, given an initial guess at the root Xi, can be used to compute a new estimate Xi+l as expressed by the iterative formula x_(i+1)=g(x_i) As with other iterative formulas in this book, the approximate error for this equation can be determined using the error estimator :

 ε=(x_(i+1)-x_i)/x_(i+1)


OPEN METHODS

For the bracketing methods in the previous chapter, the root is located within an interval prescribed by a lower and an upper bound. Repeated application of these methods always results in doser estimates of the true value of the root.

Such methods are said to be convergent because they move closer to the truth as the computation progresses. For the contrast, the open methods described in this chapter are based on formulas that require only a single starting value of x or two starting values that do not necessarily bracket the root. As such, they sometimes diverge or move away from the true root as the computation progresses.

 However, when the open methods converge, they usually do so much more quickly lhan the bracketing methods. We will begin ours discussion of open techniques with a simple version that is useful for illustrating their general form and also for demonstrating the concept of convergence.

Pitfalls of the False-Position Method

Although the false-position method would seem to always be the bracketing method of Preference, there are cases where it performs poorly. In fact, as in the following example, there are certain cases where bisection yields superior results.

THE FALSE-POSITION METHOD

Although bisection is a perfectly valid technique for determining roots, its “brute-force” approach is relatively inefficient. False position is an alternative based on a graphical insight. 

A shortcoming of the bisection method is that, in dividing the interval from Xl to Xu into equal halves, no account is taken of the magnitudes of f(Xl) and f(xu). For example, if f(Xl) is much closer to zero than f(xu), it is likely that the root is closer to xl than to Xu. An alternative method that exploits this graphical insight is to join f(Xl) and f(Xu) by a straight line. 

The intersection of this line with the x axis represents an improved estimate of the root. The fact that the replacement of the curve by a straight line gives a “false position” of the root is the origin of the name, method of false position, or in Latin. Regula falsi. It is also called the linear interpolation method.

THE BISECTION METHOD

In general, if f(x) is real and continuous in the interval from Xl to XU and f(xl) and f(xU) have opposite signs, that is, . f(xl)*f(xU ) < 0 then there is at least one real root between Xl and Xu. Incremental search methods capitalize on this observation by locating un interval where the function changes sign.

 Then the location of the sign change (and consequently, the root) is identified more precisely by dividing the interval into a number of subintervals Each of these subinterval s is searched to locate the sign change. The process is repeated and the root estimate refined by dividing the subintervals into finer increments. 

The bisection method, which is alternatively called binary chopping, interval halving, or Bolzano’s method, is one type of incremental search method in which the interval is always divided in half. If a function changes sign over an interval, the function value at the midpoint is evaluated.

 The location of the root is then determined as lying at the midpoint of the subinterval within which the sign change occurs. The process is repeated to obtain refined estimates.

GRAPHICAL METHODS

A simple method for obtaining an estimate of the root of the equation f(x) = 0 is to make a plot of the function and observe where it crosses the x axis. This point, which represents the x value for which f(x) = 0, provides a rough approximation of the root.

BRACKETING METHODS

This chapter on roots of equations deals with methods that exploit the. fact that a fuction typically changes sign in the vicinity of a root. These techniques are called bracketing methods because two initial guesses for the root are required. As the name implies, these guesses must "bracket," or be on either side of. the foal.

 The particular methods described herein employ different strategies to systematically reduce the width of the bracket and hence, home in on the correct answer. As a prelude to these techniques, we will briefly discuss graphical methods for depicting functions and their roots.

 Beyond their utility for providing rough guesses, graphical techniques are also useful for visualizing the properties of the functions and the behavior of the various numerical methods.

sábado, 15 de mayo de 2010

ROOTS OF EQUATIONS


The purpose of calculating the roots of an equation is to determine the values of x for which holds:

f (x) = 0 (28)

The determination of the roots of an equation is one of the oldest problems in mathematics and there have been many efforts in this regard. Its importance is that if we can determine the roots of an equation we can also determine the maximum and minimum, eigenvalues of matrices, solving systems of linear differential equations, etc ...


The determination of the solutions of equation (28) can be a very difficult problem. If f (x) is a polynomial function of grade 1 or 2, know simple expressions that allow us to determine its roots. For polynomials of degree 3 or 4 is necessary to use complex and laborious methods. However, if f (x) is of degree greater than four is either not polynomial, there is no formula known to help identify the zeros of the equation (except in very special cases).

There are a number of rules that can help determine the roots of an equation:

• Bolzano's theorem, which states that if a continuous function, f (x) takes on the ends of the interval [a, b] values of opposite sign, then the function accepts at least one root in that interval.
• In the case where f (x) is an algebraic function (polynomial) of degree n and real coefficients, we can say that will have n real roots or complex.
• The most important property to verify the rational roots of an algebraic equation states that if p / q is a rational root of the equation with integer coefficients:

then the denominator q divides the coefficient a and the numerator p divides the constant term a0.

Example: We intend to calculate the rational roots of the equation:
3x3 3x2 - x - 1 = 0

First, you make a change of variable x = y / 3:

and then multiply by 32:

3y2 y3-3y = -90

with candidates as a result of the polynomial are:

Substituting into the equation, we obtain that the only real root is y = -3, that is to say, (which is also the only rational root of the equation). Logically, this method is not as effective, so we can serve only as guidelines.

Most of the methods used to calculate the roots of an equation are iterative and are based on models of successive approximations. These methods work as follows: from a first approximation to the value of the root, we determine a better approximation by applying a particular rule of calculation and so on until it is determined the value of the root with the desired degree of approximation.