Companion page: Chapter 1 — Linear system of equations — the notes that this page goes with.

A short list of commands for the material of Chapter 1 — linear systems, elimination, LU, pivoting and sparse matrices. Type them out and look at the output; that is the whole point. Everything below works in GNU Octave as well, so you do not need a licence to follow along. If you want Matlab itself, sign up for a MathWorks account with your smail ID and use Matlab Online.

Two habits to form from day one: end a line with a semicolon to suppress printing (essential once matrices get large), and use help <command> or doc <command> before searching online.

Making matrices and vectors

A = [2 1 1; 4 -6 0; -2 7 2];   % rows separated by ;
b = [5; -2; 9];                % a column vector
x = [5, -2, 9];                % a row vector
size(A)                        % dimensions
A(2,3)                         % the (2,3) entry
A(:,2)                         % all of column 2
A(3,:)                         % all of row 3
A(2:3, 1:2)                    % a submatrix block

Some standard matrices:

eye(4)          % identity
zeros(3,5)      % all zeros
ones(3)         % all ones, 3x3
rand(4)         % uniform random entries in [0,1]
diag([1 2 3])   % diagonal matrix from a vector
diag(A)         % the diagonal of A, as a vector
tril(A), triu(A)   % lower and upper triangular parts

Transpose — mind the dot

A'      % conjugate transpose (Hermitian), A^H
A.'     % plain transpose, A^T

For real matrices the two agree. For complex ones they do not, and using ' when you meant .' is a common source of silent errors in Matlab code. Try it:

z = [1+2i; 3-1i];
z'      % row vector with conjugated entries
z.'     % row vector, entries unchanged
z'*z    % a real number: the squared norm
z.'*z   % a complex number: not what you want

Solving a system

x = A\b        % solve Ax = b -- this is the right way
x = inv(A)*b   % works, but slower and less accurate. Avoid.

The backslash operator does not compute an inverse. It examines the matrix, picks a suitable factorization (LU with partial pivoting for a general square matrix), and does the triangular solves. Verify with the residual:

norm(A*x - b)      % should be at the level of round off, ~1e-16

Elimination and LU

[L,U,P] = lu(A)        % gives PA = LU, with partial pivoting
P*A - L*U              % should be zero to round off
[L2,U2] = lu(A)        % here L2 = P'*L is "psychologically lower triangular"
diag(U)                % the pivots that were actually used
det(P)*prod(diag(U))   % this equals det(A)

Two things to watch here. First, diag(U) reports the pivots that partial pivoting chose, and these need not match the ones you get eliminating by hand in the order the rows happen to be written. For the \(A\) above, by hand you get pivots \((2,-8,1)\) while Matlab reports \((4,4,1)\) — both are correct, they are just different factorizations of different row orderings. Second, the product of the pivots equals \(\det A\) only up to a sign, because \(\det(PA) = \det(P)\det(A)\) with \(\det P = \pm 1\). Try lu on a matrix needing a row exchange, say [0 1; 1 0], and look at P.

Related quantities:

det(A)      % determinant -- here just as a check on the pivots
rank(A)     % rank
inv(A)      % the inverse -- for looking at, not for solving with
cond(A)     % condition number

The cost of elimination

This is the \(O(n^3)\) claim from the notes, checked on your own machine.

for n = [500 1000 2000 4000]
    A = rand(n); b = rand(n,1);
    tic; x = A\b; t = toc;
    fprintf('n = %5d   time = %8.4f s\n', n, t);
end

Each doubling of \(n\) should multiply the time by roughly 8, since \(2^3 = 8\). You will see a smaller ratio because the library spreads the work over several cores; run maxNumCompThreads(1) first and the ratios climb to 5.7, 6.8, 7.2 and onwards towards 8.

Now compare a fresh solve against one that reuses factors you already have, for many right hand sides at once:

n = 1000; A = rand(n); B = rand(n,50);
tic; X1 = A\B;                    t1 = toc;   % Matlab factors once
[L,U,P] = lu(A);
tic; X2 = U\(L\(P*B));            t2 = toc;   % reuse the factors
fprintf('%8.4f  %8.4f\n', t1, t2);

Round off and pivoting

The example from the notes, with the toy three-digit machine replaced by ordinary double precision. Since Matlab keeps about sixteen significant digits rather than three, the pivot has to be pushed down from \(10^{-4}\) to \(10^{-20}\) before the same failure appears. Carry out the elimination yourself, without exchanging rows, and watch it fail:

A = [1e-20 1; 1 1];  b = [1; 2];
m    = A(2,1)/A(1,1);            % the multiplier, 1e20
row2 = A(2,:) - m*A(1,:);        % 1 - 1e20 rounds to -1e20: the 1 is lost here
b2   = b(2)   - m*b(1);
v    = b2/row2(2)
u    = (b(1) - A(1,2)*v)/A(1,1)  % now divide by the tiny pivot
x    = A\b                       % backslash pivots, and gets it right

The hand elimination returns \((0,1)\) while the true answer is \((1,1)\) to every digit that double precision can hold. The damage is done on the line marked above: adding \(1\) to \(-10^{20}\) changes nothing at all in floating point, so everything the second equation had to say about \(u\) is thrown away before it is ever used. Exchange the two rows first and the multiplier becomes \(10^{-20}\), which is harmless.

Ill conditioning is a different problem, and no amount of pivoting cures it. Try the Hilbert matrix, which is the standard example:

H = hilb(12);  cond(H)
x_true = ones(12,1);
b = H*x_true;
norm(H\b - x_true)             % large, even though the algorithm is fine

The condition number comes out around \(10^{16}\) and the error in the recovered solution is a few parts in ten, even though we handed it the exact right hand side. Matlab will also warn that the matrix is "close to singular"; that warning is not a mistake on your part, it is the point of the example. We will come back to this once we have matrix norms.

Sparse and banded matrices

Building the tridiagonal matrix from the Poisson problem:

n = 8;
e = ones(n,1);
A = spdiags([-e 2*e -e], [-1 0 1], n, n);   % sparse tridiagonal
full(A)                                      % look at it as a dense matrix
spy(A)                                       % plot the sparsity pattern
nnz(A)                                       % number of nonzeros

Now confirm the two claims made in the notes — that the inverse of a sparse matrix is dense, and that the factors of a banded matrix stay banded:

spy(inv(full(A)))          % completely dense
[L,U] = lu(full(A));
subplot(1,2,1); spy(L);
subplot(1,2,2); spy(U);    % both still banded

And the pivots and determinant, against the formulae \(d_k = (k+1)/k\) and \(\det A = n+1\):

[L,U] = lu(full(A));
diag(U)'
det(full(A))

Finally, see how much sparse storage buys you. The dense solve should be dramatically slower:

n = 5000; e = ones(n,1);
As = spdiags([-e 2*e -e], [-1 0 1], n, n);
Ad = full(As);
b  = ones(n,1);
tic; xs = As\b; ts = toc;
tic; xd = Ad\b; td = toc;
fprintf('sparse %8.4f s   dense %8.4f s\n', ts, td);

Solving the boundary value problem

Putting the chapter together: solve \(V'' = -\rho/\epsilon_0\) on \([0,1]\) with \(V(0)=V(1)=0\) and a point charge in the middle.

n = 99; h = 1/(n+1); x = (h:h:1-h)';
e = ones(n,1);
A = spdiags([-e 2*e -e], [-1 0 1], n, n);
rho = zeros(n,1); rho(50) = 1;      % a charge at the centre
V = A\(h^2*rho);
plot(x, V, 'LineWidth', 1.5); xlabel('x'); ylabel('V(x)'); grid on

You should get a triangular potential profile with a kink at the charge, which is the Green’s function of this problem. Move the charge, add more charges, and see that the answer is the sum of the individual answers — linearity, in a picture.

Home