%% Solving Linear Systems
% This script is a tutorial for how to solve linear systems using both
% direct and iterative methods in Matlab.
%%
% A linear system of equations, of the form Ax=b, can be solved both
% directly and iteratively. Matlab has a number of different methods built into the
% function "\". To solve the system Ax=b you would call "A\b". This is by
% far the quickest way to directly solve a linear system since Matlab
% chooses the best method depending on the classification of the matrix.
N = 10;
A1=rand(N);
b1=rand(N,1);
tic
x=A1\b1;
toc
%%
% The other way of solving a linear system is by using iterative methods.
% This minimizes the residual or error to within a specified tolerance, solving the
% system.
% We can use the built in function "bicgstab" to use Biconjugate gradients
% stabilized method.
tic
bicgstab(A1,b1,1e-5,1000);
toc
%%
% Consider now if we change N to be much larger.
N = 1000;
A2=rand(N);
b2=rand(N,1);

tic
x=A2\b2;
toc

tic
bicgstab(A2,b2,1e-6,1000);
toc
%%
% As you can see, for large dense matrices, the iterative method does not
% perform anywhere close to the matrix left divide function (in fact, in
% many cases it does not even converge). However, let's now consider a
% sparse matrix.
xmin = 0;
xmax = 50;
N = 2^10;
x = linspace(xmin,xmax,N+1);
dx=x(2)-x(1); dx2=dx*dx;
e=ones(N+1,1);
Dxx = spdiags([e -2*e e], -1:1, N+1, N+1);
Dxx=(1/dx2)*Dxx;
b3=ones(size(Dxx,1),1);
%%
% This is the second order differentiation matrix from an earlier script.
% It is sparse with only entries along the diagonal. Now let's try again:
tic
x=Dxx\b3;
toc

tic
bicgstab(Dxx,b3,1e-6,1000);
toc
%%
% In general iterative methods work better on sparse matrices and ones
% with diagonal dominance. In specific cases and for particular problems
% iterative methods will be better, but in most cases you should use the
% built in "backslash" solve since it has been optimized over decades to
% use the best method for a given matrix.
%%
% If you want to see exactly what "\" is doing under the hood, you can
% change the verbosity so that it prints its choices:
spparms('spumoni',2)
A2\rand(size(A2,1),1);
Dxx\b3;
sparse(rand(N).*round(rand(N)-0.2))\rand(N,1);