1 Introduction
Let
and
a2,1×1 +a2,2×2 +···+a2,nxn = b2 . = .
an,1×1 +an,2×2 +···+an,nxn = bn a1,1 a1,2 ··· a1,n
a2,1 a2,2 ··· a2,n A=. . .. .
AM6007 Assignment 4
Write a console app (program) in C# which solves a system of linear equations using the iterative Gauss-Jacobi method.
2 Gauss-Jacobi Method
Consider a linear system of n equations with n unknowns: a1,1×1 +a1,2×2 +···+a1,nxn = b1
The system can be written in matrix/vector notation as: Ax = b
The Gauss-Jacobi method splits A as follows: a1,10···000
···00 ··· 0 0
.. .−. . . .
−a1,n . .. ..
0 a2,2 ··· 0 −a2,1 0 A=…..−..
−a1,2 · · ·
0 ··· −a2,n
. . . . . .
0 0 ··· an,n −an,1 −an,2
. . . ··· 0 0 0 ··· 0
….
an,1 an,2 b1
· · · an,n x1
b2 b=.,
x2 x= .
. . bn xn
1
That is:
where
A=D−L−U
1. D consists of the diagonal elements of A with zeroes everywhere else,
2. L consist of the negative of the lower diagonal elements of A with zeroes on the diagonal and above,
3. U consist of the negative of the upper diagonal elements of A with zeroes on the diagonal and below.
Hence our system becomes:
(D − L − U)x = b
Dx = (L+U)x+b
The inverse of D exists if none of the diagonal elements of A are zero so that: x = D−1(L + U)x + D−1b
where:
10···0 a1,1
01···0 D−1 = a2,2 .
. . .. . ….
1
an,n
xk+1 = D−1(L + U)xk + D−1b
where xk is the current estimate of the solution and xk+1 is the next estimate. Usually
taking x0 to be all zeroes is a good enough initial guess. Furthermore let: T = D−1(L + U)
c = D−1b
and notice that neither T nor c change during the calculation and can be calculated
once at the begining. The iteration simplifies to: xk+1 = Txk + c
The stopping condition for ending the iteration is usually related to how close sub- sequent estimates are, such as:
|xk+1 − xk| < ε |xk |
2
00··· The Gauss-Jacobi iterative scheme then follows:
where ε is a tolerance (say 10−7).
Thinking ahead to coding the solution, note that calculating T involves a matrix
multiplication, calculating c involves multiplying a matrix by a vector. The iteration in- volves multiplying a matrix by a vector and vector addition. Implementing the stopping condition involves vector subtraction and calculating the norm (or absolute value) of a vector.
3 Implementation in C#
Details of the requirements:
1. Read and understand the previous section.
2. Create a class named Matrix which encapsulates an n by n square matrix. 3. Create the following methods:
(a) Store the matrix data in a suitably sized 2D array.
(b) A default constuctor which sizes the matrix to 3 by 3 and sets all values to zero.
(c) A constructor public Matrix(int size);
Make sure size > 1 otherwise throw an exception. Size the matrix accordingly and set all values to zero.
(d) public static Matrix operator*(Matrix lhs, Matrix rhs)
Implements matrix multiplication and checks that the matrices are of the same size otherwise throw an exception.
(e) public static Vector operator*(Matrix lhs, Vector rhs)
Implements multiplication of a matrix by a vector and checks that the matrix and vector are of the same size otherwise throw an exception
(f) public double this[int row, int col]
Implements an indexer function to allow users to access or set elements of the matrix. Check for valid indices and throw an exception if they are invalid.
(g) public override string ToString()
Returns a string representation of the matrix for use in Console.WriteLine.
4. Create a class named Vector which encapsulates a vector of size n. 5. Create the following methods:
(a) Store the vector data in a suitably sized 1D array.
(b) A default constuctor which sizes the vector to 3 and sets all values to zero.
(c) A constructor public Vector(int size);
Make sure size > 1 otherwise throw an exception. Size the vector accordingly and set all values to zero.
3
(d) static public Vector operator +(Vector a, Vector b)
Implements vector addition and checks that the vectors are of the same size otherwise throw an exception.
(e) static public Vector operator -(Vector a, Vector b)
Implements vector subtraction and checks that the vectors are of the same size otherwise throw an exception
(f) public double this[int index]
Implements an indexer function to allow users to access or set elements of the vector. Check for valid index and throw an exception if it is invalid.
(g) public double Norm()
Implements the norm (abolute value) of a vector.
(h) public override string ToString()
Returns a string representation of the vector for use in Console.WriteLine.
6. Create a class named LinSolve which implements the Gauss-Jacobi method
7. Create the following methods/data:
(a) Provide a data member maxiters (= 100) to prevent infinite loops.
(b) public Vector Solve(Matrix A, Vector b)
Solves the matrix equation Ax = b and returns the resulting vector of the solution.
8. Your code should run against the code shown in Fig. 1. See also the output in
Fig. 2.
9. Submit your answer by putting all code into a single file with an extension .cs. Check that your code runs and works before submitting. There will be marks allocated for tidy, readable and commented code.
4
Figure 1: Test your implementation against this code.
Figure 2: Sample output.
5