## Simultaneous Linear Equations

In many problems of practical interest, we need to solve a set of linear equations simul-taneously for a given number of unknowns as in the Eq. 16 above. In this section, we discuss some numerical methods to solve systems of simultaneous linear equations (which are amenable to hand computation).

**Non-linear Fitting**

Consider the temperature dependence of the linear thermal expansion (α) of BN on tem-perature in the range 77 to 1205 K as shown in Table 1.

Table 1: The variation of linear thermal expansion with temperature in BN in the tem-perature regime (77 – 1205 K). The data is taken from [6].

Table 2: The variation of linear thermal expansion with temperature in BN in the tem-perature regime (77 – 1205 K).

Using Table. 2, one can see that the matrix reduces to

**Gaussian Elimination**

The method that we used above can also be made more algorithmic using the Gaussian elimination with partial pivoting. In this method, we make the augmented matrix which contains the matrix along with the column of the right hand side vector appended to it as shown below:

The first step is to reduce the augmented matrix to upper triangular form using the

following elimination operations:

• Interchange any rows of the augmented matrix; and

• Replace any row with a linear combination of itself and any other two;

Note that these operations leave the solution to the given set of equations unchanged. Let us consider the first row, for example. By dividing all the terms by 6, the 11 term of the matrix will become unity. Similarly, we can replace the second row using the linear combination (first row - second row divided by 3810); this leads to a zero at the 21 position and so on. This procedure can be continued till the 11, 22 and 33 become unity and the 21, 31, and 32 become zero. Then, by back substitution, we can obtain the unknowns a_{2} first, then a_{1} and finally, a_{0}.

Note that this procedure fails if any of the terms that we use for dividing the rows is zero. To avoid this, we use what is known as partial pivoting. In pivoting, in each column, we interchange rows in such a way that the largest number appears at the a_{ii}-th position. Even if the a_{ii}-th term is not zero, if pivoting is not done, the Gaussian elimination algorithm can becomes unstable – that is, the round-off errors can accumulate and lead to loss of accuracy. This is because small pivot values, when used in division, can lead to large multiplication factors which in turn produce loss of significance.

Considering the augmented matrix, thus, we first interchange the first row with the last one:

We divide the first row by 3416468; we replace the second row by (first row - second row/3810) and the third by (first row - third row/6). This leads to

Now, we exchange the second and third rows, and again divide the second row by 365.80 and replace the third row by (second row - third row/104.09). This leads to

Now, from the third row, a_{2} is directly obtained; from second row, using the a_{2}, we can obtain a_{1}; from the first row, using a_{2} and a_{1}, we obtain a_{0}. Thus, the fit is obtained to be

As compared to Eq. 23, we can see that the parameters obtained using Gaussian elimina- tion are slightly different.

**Iterative Improvement**

It is possible to improve the accuracy of the solution obtained by the Gaussian elimination method using the following iterative procedure:

• Let Ax = b be the equation for which iterative improvement is sought. Calculate

r = b − Aˆx

• Solve the equation Ae = r

• Replace ˆx by ˆx e

For example, for the Eq. 26, let us assume an approximate solution of ˆx = (−400, 0.1, 0.01).

We see that r is (−7504.5800,−9.4223, 1.3149).

Solving Ae = r gives e = (−11.071,−0.021223,−0.0070834).

This in turn leads to the improved solution, x = (−411.07, 0.078777, 0.0029166).

**Gauss-Seidel Algorithm**

There is an iterative procedure (known as Gauss-Seidel procedure) which can be used to solve the set of equations represented in Eq. 20. This procedure works better on a computer (and the algorithm looks slighlty different). However, since we will be using pen, paper and calculator to carry out the calculations, Gauss-Seidel can be carried out as follows.

Consider the Eq. 20; it can be re-written as

Now, we can start with a initial guess, say (0, 0, 0) for (a_{0}, a_{1}, a_{2}). This leads to a new set of values, namely (1299.7, 1.0899, 3.3562× 10−04). Using this as the new guess value, we get (416.45, 1.7390, 5.5709 × 10−04). We can continue in this manner till the guess value that we use and the solution that we obtain differ only by a small, predetermined error value; that is, they converge within some given error. In this case, after 482 iterations we obtain the solution, namely, (−411.39, 0.079944, 0.0029158) when the square of the distance between the old and new guess is less than 10^{−6}.

Note that the large number of iterations means that the methodology is better suited for computers than for hand computations. Further, the convergence of the solution is guaranteed only if the matrix is diagonally dominant or symmetric and positive definite.

**Design Matrix**

There is another matrix that can be generated to solve nonlinear least square fitting

problems. This matrix is called the design matrix [7]. The design matrix for the above

quadratic fit problem, for example, is as follows:

By solving the above equation for the unknowns a_{i}, one can again find the fitting param-eters.

The design matrix is not a square matrix. So, Gaussian elimination method can be used in which we make several of the rows identically zero and the rest as upper triangular – along with iterative refinement, if needed.

Design matrix is also useful in cases in which we need to fit to a specific functional form which is not a generic polynomial. For example, in the case of stress-strain curve, we want to fit the data to y = mx and not to y = mx c or, in the case of specific heat as a function of temperature at low temperatures, we want to fit the data to c_{p} = aT bT^{3}.

It is easier to set up the design matrix in such cases and solve them.

**Revisiting the Eigenvector Problem**

Let us consider the Eq. 16 for the eigenvalue 1:

**Further Reading**

Gerald and Wheatley [5] has a nice discussion on shifting with the power method which can help get the intermediate eigenvalues and the corresponding eigenvectors.