Computational Methods

Some problems have no tidy algebraic answer, and many that do have one are faster to solve by a fixed recipe. Computational methods are those recipes: finite sequences of arithmetic steps that a person with a calculator — or a machine — can carry out without thinking. This unit gives you three of them. Two solve a system of linear equations, one by systematic cancelling (Gauss elimination) and one by repeated guessing that improves itself (Gauss–Seidel); the third, the simplex method, finds the best value of a linear quantity subject to linear restrictions. The unit closes with the question every computed answer must answer: how wrong is it?

15.1 System of linear equations: Gauss Elimination Method

The augmented matrix

Write the system

as a table of its numbers only. The coefficient matrix A has entries a_{ij}; attaching the right-hand column b to it gives the augmented matrix, written [A \,|\, b]. Every step of the method is a change to the rows of this table, and the vertical bar just reminds you which column is the right-hand side.

The three elementary row operations

Only three moves are allowed, and none of them changes the solution set:

1. swap two rows (R_i \leftrightarrow R_j);
2. multiply a row by a non-zero constant (R_i \rightarrow kR_i);
3. add a multiple of one row to another (R_i \rightarrow R_i - mR_j).

The third is the workhorse. To kill the entry in column k of row i, use the multiplier

where row k is the pivot row and a_{kk} the pivot.

When the pivot is zero: partial pivoting

If a pivot comes out 0 you cannot divide by it. Swap the pivot row with a lower row whose entry in that column is non-zero — operation (1), so nothing is lost. Choosing the row with the largest entry in that column is called partial pivoting; it also keeps rounding errors small, which is why every computer program does it.

A zero pivot is not "no solution"

A zero in the pivot position only means swap rows. You may conclude something about the system only when a whole row of the coefficient part becomes zero: - 0\;0\;0\,|\,k with k \ne 0 — impossible equation, the system is inconsistent (no solution); - 0\;0\;0\,|\,0 — a redundant equation, so there are infinitely many solutions.

15.2 System of linear equations: Gauss–Seidel Method

Elimination is direct: a fixed number of steps and you are finished. The Gauss–Seidel method is iterative: you start with a guess and improve it over and over, stopping when the numbers stop changing to the accuracy you want.

Setting up the iteration

Solve each equation for the unknown on its own diagonal:

Start from x^{(0)} = y^{(0)} = z^{(0)} = 0 unless the question says otherwise.

Will it converge?

If the coefficient matrix is diagonally dominant, Gauss–Seidel converges from any starting guess. If it is not, the iteration may still converge — but it may also fly apart, and then no amount of patience helps. If the rows of the given system are not dominant, reorder the equations before you start.

The reason iteration works is easy to see with two unknowns. Solving for x moves you horizontally onto the first line; solving for y moves you vertically onto the second. The path is a staircase that closes in on the intersection.

15.3 Linear programming problems (LPP): simplex method (two variables)

A linear programming problem asks for the largest (or smallest) value of a linear expression when the variables obey linear inequalities.

The feasible region of a two-variable problem is a convex polygon, and the optimum always sits at a corner of it. That is the whole idea: instead of testing infinitely many points, test the corners.

Slack variables and standard form

An inequality cannot be put in a table, so turn each \le constraint into an equation by adding a non-negative slack variable that soaks up the unused capacity:

The tableau and the pivot rules

The starting tableau lists the constraint equations, then the objective row. The basic variables are the slacks (value d_i) and the non-basic variables x, y are zero — that is the corner (0,0).

Each tableau is one corner of the feasible region, and each pivot is a walk along one edge to a better corner. The simplex method never goes backwards, which is why it stops.

15.4 How accurate is a computed answer?

Gauss–Seidel and the simplex method both stop somewhere. A computed answer is therefore only as good as the error attached to it, and every numerical method comes with an error estimate. Two more retained methods make the point sharply; Unit 17 §17.3 revises the integration rules in exam form.

Newton–Raphson: doubling the correct digits each step

To solve f(x) = 0, take a guess x_0, follow the tangent to the curve down to the x-axis, and use where it lands as the next guess. The tangent at (x_n, f(x_n)) meets y = 0 at

Order of a quadrature rule

For \int_a^b f(x)\,dx with n equal strips of width h = \dfrac{b-a}{n} and ordinates y_0, y_1, \ldots, y_n:

The trapezium rule replaces the curve by chords, so its error is proportional to h^2; Simpson's rule replaces it by parabolas through three ordinates at a time, so its error is proportional to h^4. Halving h divides the trapezium error by 4 but the Simpson error by 16.

Chapter summary

• Gauss elimination. Form [A\,|\,b]; use R_i \rightarrow R_i - m_{ik}R_k with m_{ik} = a_{ik}/a_{kk} to make it upper triangular; back substitute. Swap rows if a pivot is 0 (partial pivoting).
• A row 0\;0\;0\,|\,k with k \ne 0 means no solution; 0\;0\;0\,|\,0 means infinitely many.
• Gauss–Seidel. x^{(k+1)} = \frac{1}{a_{11}}(b_1 - a_{12}y^{(k)} - a_{13}z^{(k)}), and so on, always using the newest value available. Converges for certain when |a_{ii}| > \sum_{j\ne i}|a_{ij}| in every row — reorder the equations if it is not so.
• Direct methods (elimination) finish in a set number of steps; iterative methods (Gauss–Seidel) are stopped when successive iterates agree to the required accuracy.
• LPP. Maximise Z = c_1x + c_2y over a convex feasible region; the optimum lies at a corner.
• Simplex. Add slacks, write Z - c_1x - c_2y = 0. Pivot column = most negative Z entry; pivot row = smallest positive ratio RHS \div column entry; pivot to make the column (1,0,\ldots,0). Stop when no Z entry is negative.
• Newton–Raphson. x_{n+1} = x_n - f(x_n)/f'(x_n); correct digits roughly double each step.
• Quadrature. Trapezium \frac{h}{2}[y_0+y_n+2(\text{rest})] with |E| \le \frac{(b-a)h^2}{12}\max|f''|; Simpson \frac{h}{3}[y_0+y_n+4(\text{odd})+2(\text{even})], n even, with |E| \le \frac{(b-a)h^4}{180}\max|f^{(4)}|.

Practice questions

Group A — Multiple choice (1 mark each)

1. In Gauss elimination the augmented matrix is reduced to [1]aa diagonal matrixban upper triangular matrixcan identity matrixda symmetric matrix
2. The Gauss–Seidel method is certain to converge if the coefficient matrix is [1]asymmetricbsingularcdiagonally dominantdtriangular
3. In a simplex tableau the entering variable is chosen from the objective row's [1]alargest positive entrybmost negative entrycsmallest ratiodlast column
4. A slack variable added to 3x + 4y \le 12 satisfies [1]as \le 0bs = 0cs \ge 0ds = 12
5. If a row of the reduced augmented matrix reads 0\;0\;0\,|\,5, the system has [1]aa unique solutionbno solutioncinfinitely many solutionsdexactly two solutions
6. Simpson's rule can be applied only when the number of strips n is [1]aoddbevenca multiple of 3dprime
7. For the Newton–Raphson method, x_{n+1} equals [1]ax_n + \frac{f(x_n)}{f'(x_n)}bx_n - \frac{f'(x_n)}{f(x_n)}cx_n - \frac{f(x_n)}{f'(x_n)}d\frac{f(x_n)}{f'(x_n)}

Answers to Group A

1. (b) — forward elimination puts zeros below the diagonal so the last row has one unknown.

2. (c) — diagonal dominance guarantees convergence from any starting guess.

3. (b) — the most negative coefficient promises the fastest increase in Z.

4. (c) — slack measures unused capacity, which cannot be negative.

5. (b) — it states 0 = 5, an impossibility, so the system is inconsistent.

6. (b) — the ordinates are taken three at a time (two strips per parabola).

7. (c) — the tangent at x_n meets the x-axis at x_n - f(x_n)/f'(x_n).

Group B — Short answer (5 marks each)

1. Solve by Gauss elimination: \;x + 2y + 3z = 14, \;2x + 5y + 2z = 18, \;3x + y + 5z = 20. [5]
2. State the three elementary row operations and explain, with the system 2x + 4y = 6, \;3x + y = -1, why a row interchange can be necessary during Gauss elimination. Solve the system. [5]
3. Carry out three iterations of the Gauss–Seidel method on \;8x - 3y + 2z = 20, \;4x + 11y - z = 33, \;6x + 3y + 12z = 35, starting from (0,0,0). Give your answers to four decimal places. [5]
4. Test whether \;x + 5y = 7, \;4x + y = 6 is suitable for Gauss–Seidel as written. Rearrange if necessary and perform three iterations. [5]
5. Maximise Z = 5x + 7y subject to x + y \le 4, 3x + 8y \le 24, x, y \ge 0 by the simplex method. [5]
6. Using the Newton–Raphson method with x_0 = 1, find the root of x^3 + x - 3 = 0 correct to four decimal places. [5]

Answers to Group B

1. Augmented matrix rows (1,2,3\,|\,14), (2,5,2\,|\,18), (3,1,5\,|\,20).

m_{21} = 2: R_2 - 2R_1 = (0,1,-4\,|\,-10). m_{31} = 3: R_3 - 3R_1 = (0,-5,-4\,|\,-22). m_{32} = -5: R_3 + 5R_2 = (0,0,-24\,|\,-72).

Back substitution: -24z = -72 \Rightarrow z = 3; y - 4(3) = -10 \Rightarrow y = 2; x + 2(2) + 3(3) = 14 \Rightarrow x = 1.

x = 1,\quad y = 2,\quad z = 3

Check in equation 2: 2(1) + 5(2) + 2(3) = 2 + 10 + 6 = 18. Correct.

2. The operations are (i) R_i \leftrightarrow R_j, (ii) R_i \rightarrow kR_i with k \ne 0, (iii) R_i \rightarrow R_i - mR_j. None of them changes the solution set.

A swap is needed whenever the intended pivot a_{kk} is zero, because the multiplier m_{ik} = a_{ik}/a_{kk} would require division by zero. (Swapping to the largest available entry — partial pivoting — also limits rounding error.)

Here the pivot 2 is non-zero, so no swap is needed: m_{21} = \frac32 gives R_2 - \frac32 R_1 = (0, -5\,|\,-10), so y = 2 and then 2x + 8 = 6 \Rightarrow x = -1. So x = -1, y = 2.

3. The rows are dominant (8 > 5, 11 > 5, 12 > 9). The formulas are x = \frac{20 + 3y - 2z}{8}, y = \frac{33 - 4x + z}{11}, z = \frac{35 - 6x - 3y}{12}.

k	x^{(k)}	y^{(k)}	z^{(k)}
1	2.5000	2.0909	1.1439
2	2.9981	2.0138	0.9142
3	3.0266	1.9825	0.9077

Row 1 in detail: x = \frac{20+0-0}{8} = 2.5; y = \frac{33-4(2.5)+0}{11} = \frac{23}{11} = 2.0909; z = \frac{35-6(2.5)-3(2.0909)}{12} = \frac{13.7273}{12} = 1.1439.

Continuing gives x \to 3.0168, y \to 1.9859, z \to 0.9118; after three iterations every value is already correct to one decimal place.

4. Row 1 has |1| < |5|, so the system is not diagonally dominant as written and Gauss–Seidel is unsafe. Interchange the equations:

4x + y = 6, \qquad x + 5y = 7

Now 4 > 1 and 5 > 1: both rows dominant. With x = \frac{6-y}{4}, y = \frac{7-x}{5} from (0,0):

k	x^{(k)}	y^{(k)}
1	1.500000	1.100000
2	1.225000	1.155000
3	1.211250	1.157750

The exact solution is x = \frac{23}{19} = 1.210526, y = \frac{22}{19} = 1.157895.

5. Slacks: x+y+s_1 = 4, 3x+8y+s_2 = 24, and Z - 5x - 7y = 0.

basic	x	y	s_1	s_2	RHS	ratio
s_1	1	1	1	0	4	4
s_2	3	\mathbf{8}	0	1	24	3 ←
Z	-5	-7	0	0	0

y enters, s_2 leaves, pivot 8. After R_2 \div 8, R_1 - R_2, Z + 7R_2:

basic	x	y	s_1	s_2	RHS	ratio
s_1	\mathbf{\frac58}	0	1	-\frac18	1	\frac85 ←
y	\frac38	1	0	\frac18	3	8
Z	-\frac{19}{8}	0	0	\frac78	21

x enters, s_1 leaves, pivot \frac58. After R_1 \times \frac85, R_2 - \frac38 R_1, Z + \frac{19}{8}R_1:

basic	x	y	s_1	s_2	RHS
x	1	0	\frac85	-\frac15	\frac85
y	0	1	-\frac35	\frac15	\frac{12}{5}
Z	0	0	\frac{19}{5}	\frac25	\frac{124}{5}

x = \frac85 = 1.6,\qquad y = \frac{12}{5} = 2.4,\qquad Z_{\max} = \frac{124}{5} = 24.8

Check: 1.6 + 2.4 = 4 and 3(1.6) + 8(2.4) = 4.8 + 19.2 = 24 — both tight.

6. f(x) = x^3 + x - 3, f'(x) = 3x^2 + 1. f(1) = -1 < 0, f(2) = 7 > 0, so the root lies in (1,2).

n	x_n	f(x_n)	f'(x_n)	x_{n+1}
0	1.000000	-1.000000	4.000000	1.250000
1	1.250000	0.203125	5.687500	1.214286
2	1.214286	0.004738	5.423469	1.213412
3	1.213412	0.000003	5.417107	1.213412

x = 1.2134 \quad (4 \text{ d.p.})

Group C — Long answer (8 marks each)

1. (a) Describe the Gauss elimination method, explaining the roles of the pivot, the multiplier and back substitution. [3](b) Hence solve \;2x + y + z = 10, \;3x + 2y + 3z = 18, \;x + 4y + 9z = 16. [5]
2. (a) Write down the Gauss–Seidel iteration formulas for a 3 \times 3 system and state the condition for convergence. Contrast a direct with an iterative method. [3](b) Solve, correct to three decimal places, \;10x + y - z = 11.19, \;x + 10y + z = 28.08, \;-x + y + 10z = 35.61. [5]
3. (a) Explain what slack variables are and why the optimum of an LPP must occur at a corner of the feasible region. [3](b) Maximise Z = 4x + 3y subject to 2x + y \le 10, x + 3y \le 15, x, y \ge 0 by the simplex method, and verify your answer by evaluating Z at every corner of the feasible region. [5]

Answers to Group C

1. (a) Write the system as the augmented matrix [A\,|\,b]. For each column k in turn, the diagonal entry a_{kk} is the pivot and its row is the pivot row; for every row i below it form the multiplier m_{ik} = a_{ik}/a_{kk} and replace R_i by R_i - m_{ik}R_k, which makes the entry in column k zero. If a pivot is zero, interchange with a lower row (partial pivoting). After all columns are treated the matrix is upper triangular, and back substitution solves the last equation for one unknown, substitutes it into the row above, and works upwards.

(b) Rows (2,1,1\,|\,10), (3,2,3\,|\,18), (1,4,9\,|\,16).

m_{21} = \frac32: R_2 - \frac32R_1 gives the row 0,\ \frac12,\ \frac32\ |\ 3. m_{31} = \frac12: R_3 - \frac12R_1 gives the row 0,\ \frac72,\ \frac{17}{2}\ |\ 11.

m_{32} = \frac72 \div \frac12 = 7: R_3 - 7R_2 = (0, 0, -2\,|\,-10).

Back substitution: -2z = -10 \Rightarrow z = 5; \frac12 y + \frac{15}{2} = 3 \Rightarrow \frac12 y = -\frac92 \Rightarrow y = -9; 2x - 9 + 5 = 10 \Rightarrow 2x = 14 \Rightarrow x = 7.

x = 7,\qquad y = -9,\qquad z = 5

Check in equation 3: 7 + 4(-9) + 9(5) = 7 - 36 + 45 = 16. Correct.

2. (a) With x^{(0)} = y^{(0)} = z^{(0)} = 0,

x^{(k+1)} = \frac{b_1 - a_{12}y^{(k)} - a_{13}z^{(k)}}{a_{11}},\quad y^{(k+1)} = \frac{b_2 - a_{21}x^{(k+1)} - a_{23}z^{(k)}}{a_{22}},\quad z^{(k+1)} = \frac{b_3 - a_{31}x^{(k+1)} - a_{32}y^{(k+1)}}{a_{33}}

Convergence is guaranteed when the coefficient matrix is diagonally dominant, |a_{ii}| > \sum_{j \ne i}|a_{ij}| for every i. A direct method such as Gauss elimination reaches the exact answer in a fixed, known number of arithmetic steps; an iterative method produces a sequence of approximations and is stopped when two successive ones agree to the wanted accuracy, so its cost depends on the accuracy demanded and on how fast it converges.

(b) All rows are dominant (10 > 2 in each). Using x = \frac{11.19 - y + z}{10}, y = \frac{28.08 - x - z}{10}, z = \frac{35.61 + x - y}{10}:

k	x^{(k)}	y^{(k)}	z^{(k)}
1	1.119000	2.696100	3.403290
2	1.189719	2.348699	3.445102
3	1.228640	2.340626	3.449801
4	1.229918	2.340028	3.449989
5	1.229996	2.340001	3.449999

x = 1.230,\qquad y = 2.340,\qquad z = 3.450

Check in equation 3: -1.23 + 2.34 + 10(3.45) = 35.61. Correct.

3. (a) A slack variable s_i \ge 0 is added to a \le constraint to turn it into an equation, a_ix + b_iy + s_i = d_i; its value is the amount of that resource left unused. Because every constraint is linear, the feasible region is the intersection of half-planes and is therefore a convex polygon. The objective Z = c_1x + c_2y is constant along a family of parallel lines; pushing such a line in the direction of increasing Z, the last point of the polygon it meets is either a vertex or (if the line is parallel to an edge) a whole edge whose endpoints are vertices. Either way a vertex attains the maximum, so only the corners need testing.

(b) Slacks: 2x + y + s_1 = 10, x + 3y + s_2 = 15, and Z - 4x - 3y = 0.

basic	x	y	s_1	s_2	RHS	ratio
s_1	\mathbf{2}	1	1	0	10	5 ←
s_2	1	3	0	1	15	15
Z	-4	-3	0	0	0

x enters, s_1 leaves, pivot 2. After R_1 \div 2, R_2 - R_1, Z + 4R_1:

basic	x	y	s_1	s_2	RHS	ratio
x	1	\frac12	\frac12	0	5	10
s_2	0	\mathbf{\frac52}	-\frac12	1	10	4 ←
Z	0	-1	2	0	20

y enters, s_2 leaves, pivot \frac52. After R_2 \times \frac25, R_1 - \frac12R_2, Z + R_2:

basic	x	y	s_1	s_2	RHS
x	1	0	\frac35	-\frac15	3
y	0	1	-\frac15	\frac25	4
Z	0	0	\frac95	\frac25	24

No negative entry remains, so x = 3, y = 4, Z_{\max} = 24.

Corner check. The lines 2x+y=10 and x+3y=15 meet where x = 3, y = 4. The corners are (0,0), (5,0), (3,4), (0,5):

corner	(0,0)	(5,0)	(3,4)	(0,5)
Z = 4x+3y	0	20	\mathbf{24}	15

The maximum is 24 at (3,4), agreeing with the simplex result.