13. Package Scipy and More Tools for Linear Algebra¶
This work is licensed under Creative Commons Attribution-ShareAlike 4.0 International
13.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
13.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.
13.1.2. Exercise D¶
Check this result by computing the residual vector \(r = A x - b\).
13.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|\).
13.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).
13.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.
13.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\).