Package Scipy and More Tools for Linear Algebra

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

\[\begin{split} A = \left[ \begin{array}{ccc} 4. & 2. & 1. \\ 9. & 3. & 1. \\ 25. & 5. & 1. \end{array} \right] \end{split}\]

and

\[b = [0.693147, 1.098612, 1.386294]\]

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:

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

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

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

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

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

  6. 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:

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

  2. 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\).

Contents