%% Solving non-singular linear systems %% Introduction % % In this lab, we will learn how to use Matlab to solve linear systems of % the form $$Ax = b$$, where the matrix % $A$ is an $N % \times N$ matrix, and $b$ is an % $N \times 1$ vector. The matrix % $A$ is assumed to be non-singular and so % the linear system has exactly one solution $x$. % clear all; format short g; %% Backslash operator : x = A\b % % Linear systems arise in all areas of engineering, science and % computational mathematics, and as a result, significant research has gone % into finding the most computationally efficient algorithms for solving % linear systems. In this tutorial, we will investigate some standard % solvers for square, non-singular systems. In a later tutorial, we will % investigate the solution to over-determined and under-determined systems. %

% The first linear system that most students encounter is the one % whose solution describes the point at the intersection of two lines. For % example, suppose we have the equations for two lines, given by % $$ % \begin{aligned} % 3x - y = 1 \\ % 2x + y = 7 % \end{aligned} % $$ % % Using matrix notation, we can write this as the $2 % \times 2$ linear system given by % $$ % \left[\begin{array}{rr} 3 & -1 \\ 2 & 1 \end{array}\right] % \left[\begin{array}{r} x \\ y \end{array}\right] = % \left[\begin{array}{r} 1 \\ 7 \end{array}\right]. % $$ % More generally, a linear system is written as % $$A{\bf x} = {\bf b}$$ % where $A$ is an $n % \times n$ matrix, and ${\bf x}$ and % ${\bf b}$ are both $n % \times 1$ vectors. Formally, the solution to a linear system can % be written as % $${\bf x} = A^{-1}{\bf b}$$ % In Matlab, we can solve this system using the "backslash" operator. % A = [3 -1; -2 -1] %% b = [1; -7]; %% x = A\b %% % % Checking this answer, we see that the residual error for this problem % is essentially zero % r = b - A*x %% The LU Decomposition % % Behind the scenes, what is the "backslash" operator doing? % If $A$ % is a square, non-singular matrix, the backslash operator uses Gaussian % Elimination, a procedure very similar to what you might do by hand to % solve a linear system. Basically, the computer eliminates % variables in order, until only a single equation for the last variable % remains. Then, working backwards, each variable is solved in order. This % procedure is carried out using the $LU$ decomposition. %

% Essentially, the backslash operator is doing the following behind the % scenes : % format short A = [1,2,3; 4,-1,5; 3,1,0]; b = [1;1;1]; %% [L,U,P] = lu(A); %% % % Here, $L$ is a lower triangular % matrix, $U$ is an upper triangular % matrix and $P$ is a permutation matrix. % These three matrices are related as % $$ % PA = LU % $$ %
% These matrices are % L %% U %% P*A %% L*U %% % % We can then use these to solve for $x$ % as follows. % format long y = L\(P*b); % "Forward substitution" - P*A = P*b x = U\y % "Backward substitution" %% % % We can compare this solution to the solution using the backslash % operator. % A\b %% % % The solutions are exactly the same. % %% Using linsolve and mldivide % % Most Matlab operators have equivalent functions that can be used when % using an operator ("backslash" in this case) is inconvenient. For % example, you cannot pass a "backslash" operator as an argument to a % function. In this case, the linsolve or the % mldivide functions may be more % convenient ways to obtain essentially the same solution as that obtained % using backslash. % A = [1,2,3; 4,-1,5; 3,1,0]; b = [1;1;1]; x = linsolve(A,b) % Return the same solution as backslash %% % % A second function, mldivide (for "matrix left divide") % is exactly the functional equivalent of backslash % A\b - mldivide(A,b) %% Solving sparse linear systems % % If the matrix $A$ is stored in sparse % matrix format, then the backslash operator will take advantage of the % sparse format and solve the system much more quickly that the standard % full matrix solve. To see this, will solve the trivial linear system % $Ix = b$ where $I$ is the indentity matrix. % I = eye(5000); b = ones(5000,1); %% % % We will time both solves using tic and toc and see which solve is faster. % % Time the solve using full matrix. tic x = I\b; toc %% % Create a sparse matrix that stores only the non-zero entries. I = sparse(eye(5000)); tic x = I\b; toc %% % % The sparse format is orders of magnitude times faster than the full format. This is % due to the fact that the sparse format is able to avoid most % multiplication by zeros. %

% An easy way to create banded sparse matrices is to use the spdiags command. % %% Why don't we just use inv? % % It is very tempting to use the inv function to % solve a linear system. After all, when we write down the solution to a % linear system using mathematical notation, we write $${\bf x} = A^{-1}{\bf b}$$ So, it seems % natural to use the inv function. % A = [1,2,3; 4,-1,5; 3,1,0]; b = [1;1;1]; x = inv(A)*b %% % % For this simple problem, the solution is identical to what we got % above using backslash. But there are examples of matrices for which % the solution using inv will be much worse than that using % the backslash operator. For this reason, the inv should be avoided. %

% Here is an example of the famous Hilbert matrix, % introduced by the German mathematician David Hilbert in 1894. This matrix % is an example of particular ill-conditioned matrix, meaning that % it is very hard to solve a linear systems accurately that involves this matrix. %

% In the following example, we compute the solution to a Hilbert system % using both backslash and the inverse operator, and compare the results to % that obtained using the exact solution. % format short A = hilb(10) b = ones(10,1); %% % % We can compare the exact solution (which we can obtain using the % analytically known exact inverse)to the solutions obtained using % backslash and inv. To compare the % differences between solutions, we use the norm % function norm(e,Inf), defined as % $$ % {\tt norm(e,Inf)} = \parallel \! u \! \parallel_\infty = \max_{i} |e_i| % $$ % x_true = invhilb(10)*b; error_backslash = norm(A\b-x_true,Inf)/norm(x_true,Inf); error_inv = norm(inv(A)*b - invhilb(10)*b,Inf)/norm(x_true,Inf); %% % % Here are the results. % fprintf('Error using backslash : %12.4e\n',error_backslash); fprintf('Error using inv : %12.4e\n',error_inv); %% % % The solution obtained from the backslash operator is several order of % magnitude better than the the solution obtained using the inverse % operator. %

% Using x = A^(-1)*b results in % the same error as inv operator, and should % also be avoided. % %% Lab exercises % %
% Solve for the unknowns in the following problems. %
    %
  1. $$\begin{eqnarray} % x - y & = & 5 \\ % 3x + 5y & = & 1 % \end{eqnarray} % $$
  2. %
  3. $$\begin{eqnarray} % x - y + 3z & = & 5 \\ % 3x + 5y - z & = & 1 \\ % 5y - 3z & = & -22 \\ % \end{eqnarray} % $$
  4. %
  5. % $$\begin{eqnarray} % a & = & 4(b - c) \\ % \frac{b-a}{2c} & = & 1 \\ % a + b & = & c-2 % \end{eqnarray} % $$
  6. %
  7. Solve the linear system % $$\begin{eqnarray} % ax_1^2 + b x_1 + c & = & y_1 \\ % ax_2^2 + b x_2 + c & = & y_2 \\ % ax_3^2 + b x_3 + c & = & y_3 % \end{eqnarray} % $$ % for unknown parameters % $a$, % $b$ and % $c$, % where points % $(x_1,y_1)$, % $(x_2,y_2)$, % $(x_3,y_3)$ are known and are given by % $(1,3)$, % $(2,5)$, % $(3,-1)$, respectively. %
  8. %
  9. Solve the linear system $A{\bf x} % = {\bf b}$, where $A$ is a $5 \times 5$ matrix with ones on the sub- and % super-diagonal, and -2 on the diagonal, i.e. % $$ % A = \left[\begin{array}{rrrrr} % -2 & 1 & 0 & 0 & 0 \\ % 1 & -2 & 1 & 0 & 0 \\ % 0 & 1 & -2 & 1 & 0 \\ % 0 & 0 & 1 & -2 & 1 \\ % 0 & 0 & 0 & 1 & -2 \\ % \end{array}\right]. % $$ % The right hand side vector ${\bf b}$ is a % $5 \times 1$ column vector of ones. % Create the system using Matlab functions (1) diag and (2) spdiag.
  10. %
  11. % Often, linear systems arise in the context of solving % recurrence relations. For example, the famous % Fibonacci sequence, named after the 13th century Italian mathematician % Fibonacci, is generated from the recurrence relation % $$ % x_{j+1} = x_j + x_{j-1}, \qquad j = 1,2,3,4,... % $$ % where $x_0 = 0$ and $x_1 = 1$. The first few numbers % in the sequence are $[0,1,1,2,3,5,8,13,21,34,...]$. Construct a linear % system that relates values $x_j$, % $j = 1,2,3,4,5...$ to the known % starting values $x_0$ and % $x_1$. Solve % the linear system to generate the first 10 Fibonacci numbers. % Try this using both diag and % spdiag. %

    % Hint: Write out the equations % for enough values of $j$ so that you % can see a pattern in the linear system. Then construct the system.
  12. %
  13. Suppose you only knew $x_0$ % and $x_{11}$. Could you still % generate the Fibonacci sequence? Modify the system you found above to % solve for the first 10 values of the sequence, using only the known % values $x_0 = 0$ and % $x_{11} = 89$.
  14. %
%
% %% % %

Compare your answers with the solutions.

%