To install click the Add extension button. That's it.
The source code for the WIKI 2 extension is being checked by specialists of the Mozilla Foundation, Google, and Apple. You could also do it yourself at any point in time.
How to transfigure the Wikipedia
Would you like Wikipedia to always look as professional and up-to-date? We have created a browser extension. It will enhance any encyclopedic page you visit with the magic of the WIKI 2 technology.
Try it — you can delete it anytime.
Install in 5 seconds
Yep, but later
4,5
Kelly Slayton
Congratulations on this excellent venture… what a great idea!
Alexander Grigorievskiy
I use WIKI 2 every day and almost forgot how the original Wikipedia looks like.
Let be a square system of n linear equations, where:
When and are known, and is unknown, we can use the Gauss–Seidel method to approximate . The vector denotes our initial guess for (often for ). We denote as the k-th approximation or iteration of , and is the next (or k+1) iteration of .
Matrix-based formula
The solution is obtained iteratively via
where the matrix is decomposed into a lower triangular component , and a strictly upper triangular component such that .[4] More specifically, the decomposition of into and is given by:
Why the matrix-based formula works
The system of linear equations may be rewritten as:
The Gauss–Seidel method now solves the left hand side of this expression for , using previous value for on the right hand side. Analytically, this may be written as:
Element-based formula
However, by taking advantage of the triangular form of , the elements of can be computed sequentially for each row using forward substitution:[5]
Notice that the formula uses two summations per iteration which can be expressed as one summation that uses the most recently calculated iteration of . The procedure is generally continued until the changes made by an iteration are below some tolerance, such as a sufficiently small residual.
Discussion
The element-wise formula for the Gauss–Seidel method is similar to that of the Jacobi method.
The computation of uses the elements of that have already been computed, and only the elements of that have not been computed in the (k+1)-th iteration. This means that, unlike the Jacobi method, only one storage vector is required as elements can be overwritten as they are computed, which can be advantageous for very large problems.
However, unlike the Jacobi method, the computations for each element are generally much harder to implement in parallel, since they can have a very long critical path, and are thus most feasible for sparse matrices. Furthermore, the values at each iteration are dependent on the order of the original equations.
The Gauss–Seidel method sometimes converges even if these conditions are not satisfied.
Golub and Van Loan give a theorem for an algorithm that splits into two parts. Suppose is nonsingular. Let be the spectral radius of . Then the iterates defined by converge to for any starting vector if is nonsingular and .[8]
Algorithm
Since elements can be overwritten as they are computed in this algorithm, only one storage vector is needed, and vector indexing is omitted. The algorithm goes as follows:
algorithm Gauss–Seidel method isinputs:A, boutput:φChoose an initial guess φ to the solutionrepeat until convergence
forifrom 1 untilndoσ ← 0forjfrom 1 untilndoifj ≠ ithenσ ← σ + aijφjend ifend (j-loop)
φi ← (bi − σ) / aiiend (i-loop)
check if convergence is reached
end (repeat)
Examples
An example for the matrix version
A linear system shown as is given by:
We want to use the equation
in the form
where:
We must decompose into the sum of a lower triangular component and a strict upper triangular component :
The inverse of is:
Now we can find:
Now we have and and we can use them to obtain the vectors iteratively.
First of all, we have to choose : we can only guess. The better the guess, the quicker the algorithm will perform.
We choose a starting point:
We can then calculate:
As expected, the algorithm converges to the exact solution:
In fact, the matrix A is strictly diagonally dominant (but not positive definite).
Another example for the matrix version
Another linear system shown as is given by:
We want to use the equation
in the form
where:
We must decompose into the sum of a lower triangular component and a strict upper triangular component :
The inverse of is:
Now we can find:
Now we have and and we can use them to obtain the vectors iteratively.
First of all, we have to choose : we can only guess. The better the guess, the quicker will perform the algorithm.
We suppose:
We can then calculate:
If we test for convergence we'll find that the algorithm diverges. In fact, the matrix A is neither diagonally dominant nor positive definite.
Then, convergence to the exact solution
is not guaranteed and, in this case, will not occur.
An example for the equation version
Suppose given k equations where xn are vectors of these equations and starting point x0.
From the first equation solve for x1 in terms of For the next equations substitute the previous values of xs.
To make it clear consider an example.
Solving for and gives:
Suppose we choose (0, 0, 0, 0) as the initial approximation, then the first approximate solution is given by
Using the approximations obtained, the iterative procedure is repeated until the desired accuracy has been reached. The following are the approximated solutions after four iterations.
0.6
2.32727
−0.987273
0.878864
1.03018
2.03694
−1.01446
0.984341
1.00659
2.00356
−1.00253
0.998351
1.00086
2.0003
−1.00031
0.99985
The exact solution of the system is (1, 2, −1, 1).
An example using Python and NumPy
The following numerical procedure simply iterates to produce the solution vector.
importnumpyasnpITERATION_LIMIT=1000# initialize the matrixA=np.array([[10.0,-1.0,2.0,0.0],[-1.0,11.0,-1.0,3.0],[2.0,-1.0,10.0,-1.0],[0.0,3.0,-1.0,8.0],])# initialize the RHS vectorb=np.array([6.0,25.0,-11.0,15.0])print("System of equations:")foriinrange(A.shape[0]):row=[f"{A[i,j]:3g}*x{j+1}"forjinrange(A.shape[1])]print("[{0}] = [{1:3g}]".format(" + ".join(row),b[i]))x=np.zeros_like(b,np.float_)forit_countinrange(1,ITERATION_LIMIT):x_new=np.zeros_like(x,dtype=np.float_)print(f"Iteration {it_count}: {x}")foriinrange(A.shape[0]):s1=np.dot(A[i,:i],x_new[:i])s2=np.dot(A[i,i+1:],x[i+1:])x_new[i]=(b[i]-s1-s2)/A[i,i]ifnp.allclose(x,x_new,rtol=1e-8):breakx=x_newprint(f"Solution: {x}")error=np.dot(A,x)-bprint(f"Error: {error}")