# 15. Package Scipy and More Tools for Linear Algebra#

This work is licensed under Creative Commons Attribution-ShareAlike 4.0 International

## 15.1. Introduction#

The package package *Scipy* provides a a great array of funtions for scinetific computing;
her we wil just exlore one part of it: some additional tools for linear algebra from module *linalg* within the package *Scipy*.
This provides tools for solving simultaneous linear equations, for variations on the LU factorization seen in a numerical methods course, and much more.

This module has the standard standard nickname `la`

, so import it using that:

```
import scipy.linalg as la
```

SciPy usually needs stuff from NumPy, so let’s import that also:

```
import numpy as np
```

### 15.1.1. Exercise C#

Use the Scipy function `solve`

to solve \(A x = b\) for \(x\),
where

and

as in Exercise A of section Numpy Array Operations and Linear Algebra.

### 15.1.2. Exercise D#

Check this result by computing the *residual vector* \(r = A x - b\).

#### 15.1.2.1. Exercise D-bonus#

Check further by computing the maximum error, or *maximum norm*, or *infinity norm* of this: \(\|r\|_\infty = \|A x - b\|_\infty\);

that is, the maximum of the absolute values of the elements of \(r\), \(\max_{i=1}^n |r_i|\).

### 15.1.3. Exercise E#

Next use the Scipy function `lu`

to compute the \(PA = LU\) factorization of \(A\).
Check you work by verifying that:

\(P\) is a permutation matrix (a one in each row and column; the rest all zeros)

\(PA\) is got by permutating the rows of \(A\)

\(L\) is (unit) lower triangular (all zeros above the main diagonal; ones on the main diagonal)

\(U\) is upper triangular (all zeros below the main diagonal)

The products \(PA\) and \(LU\) are the same (or very close; there might be rounding error).

The “residual matrix” \(R = PA - LU\) is zero (or very close; there might be rounding error).

### 15.1.4. Exercise F#

Use the above arrays \(b\), \(P\), \(L\) and \(U\) to solve \(AX = b\) via \(LUx = PAx = Pb\), as follows:

Solve \(Lc = Pb\) for \(c\).

Solve \(Ux = c\) for \(x\).

At each step verify by looking at the errors, or *residuals,*: \(Lc - Pb\), \(Ux - c\), \(Ax - b\).
Better yet, compute the maximum norm of each of these errors.

### 15.1.5. Optional Exercise G (on the matrix norm, seen in a numerical methods course)#

Try to work out how to compute the matrix norm \(\| A \|_\infty\).

Test your method on \(A\), \(L\) and \(U\), and compare \(\| A \|_\infty\), \(\| PA \|_\infty\), and \(\| L \|_\infty \| U \|_\infty\).