year = (1980:10:2010)'  
pop = [984.736; 1148.364; 1263.638; 1330.141];

t = year - 1980;
y = pop;

V = zeros(4,4);
for i = 1:4
    V(i,:) = [1 t(i) t(i)^2 t(i)^3];
end
V

c = V \ y

c = c(end:-1:1);       % reverse the ordering

polyval(c,2005-1980)   % apply the 1980 time shift

plot(1980+t,y,'o')

plot(1980+t,y,'o')
hold on
tt = linspace(0,30,300)';   % 300 times from 1980 to 2010
yy = polyval(c,tt);         % evaluate the cubic
plot(1980+tt,yy)
hold off

clf   % clear figure
plot(1980+t,y,'o')
hold on
tt = linspace(-10,50,300)';   
plot(1980+tt,polyval(c,tt))
title('Population of China')
xlabel('year'), ylabel('population (millions)')
legend('data','interpolant','location','northwest')

A = [ 1, 2, 3, 4, 5; 50 40 30 20 10
    pi, sqrt(2), exp(1), (1+sqrt(5))/2, log(3) ]

[m,n] = size(A)

x = [ 3; 3; 0; 1; 0 ]
size(x)

AA = [ A; A ]
B = [ zeros(3,2), ones(3,1) ]

A'

x'

row = 1:4              % start:stop
col = ( 0:3:12 )'      % start:step:stop

s = linspace(-1,1,5)'  % start,stop,number

a = A(2,end-1)

x(2)

A(1:2,end-2:end)    % first two rows, last three columns

A(:,1:2:end)        % all of the odd columns

B = diag( [-1 0 -5] )     % create a diagonal matrix

BA = B*A     % matrix product

try A*B, catch lasterr, end
disp('Error using  *')   % ignore this line
disp(lasterr.message)   % ignore this line

B^3    % same as B*B*B

C = -A;

elementwise = A.*C

xtotwo = x.^2

twotox = 2.^x

cos(pi*x')

A = magic(3)
b = [1;2;3];
x = A\b

residual = b - A*x

A = [0 1; 0 0];   % known to be singular
b = [1;2];
x = A\b

A = magic(5)
L = tril(A)

b = ones(5,1);
x = forwardsub(L,b)

b - L*x

alpha = 0.3;
beta = 2.2;
U = eye(5) + diag([-1 -1 -1 -1],1);
U(1,[4 5]) = [ alpha-beta, beta ];

x_exact = ones(5,1);
b = [alpha;0;0;0;1];

function x = backsub(U,b)
% BACKSUB   Solve an upper triangular linear system.
% Input:
%   U    upper triangular square matrix (n by n)
%   b    right-hand side vector (n by 1)   
% Output:
%   x    solution of Ux=b (n by 1 vector)

n = length(U);
x = zeros(n,1);
for i = n:-1:1
  x(i) = ( b(i) - U(i,i+1:n)*x(i+1:n) ) / U(i,i);
end

end

x = backsub(U,b);
err = x - x_exact

alpha = 0.3;
beta = 1e12;
U = eye(5) + diag([-1 -1 -1 -1],1);
U(1,[4 5]) = [ alpha-beta, beta ];
b = [alpha;0;0;0;1];

x = backsub(U,b);
err = x - x_exact

A = [
     2    0    4    3 
    -4    5   -7  -10 
     1   15    2   -4.5
    -2    0    2  -13
    ];

b = [ 4; 9; 29; 40 ];

S = [A, b]

mult21 = S(2,1)/S(1,1)
S(2,:) = S(2,:) - mult21*S(1,:)

mult31 = S(3,1)/S(1,1)
S(3,:) = S(3,:) - mult31*S(1,:);
mult41 = S(4,1)/S(1,1)
S(4,:) = S(4,:) - mult41*S(1,:);
S

for i = 3:4
    mult = S(i,2)/S(2,2);
    S(i,:) = S(i,:) - mult*S(2,:);
end
S

for i = 4
    mult = S(i,3)/S(3,3);
    S(i,:) = S(i,:) - mult*S(3,:);
end
S

U = S(:,1:4)
z = S(:,5)

x = backsub(U,z)

b - A*x

A = [2 0 4 3 ; -4 5 -7 -10 ; 1 15 2 -4.5 ; -2 0 2 -13];

I = eye(4);

mult21 = A(2,1)/A(1,1);
L21 = I - mult21*I(:,2)*I(:,1)';
A = L21*A

mult31 = A(3,1)/A(1,1);
L31 = I - mult31*I(:,3)*I(:,1)';
A = L31*A;

mult41 = A(4,1)/A(1,1);
L41 = I - mult41*I(:,4)*I(:,1)';
A = L41*A

function [L,U] = lufact(A)
% LUFACT   LU factorization (demo only--not stable!).
% Input:
%   A    square matrix
% Output:
%   L,U  unit lower triangular and upper triangular such that LU=A

n = length(A);
L = eye(n);   % ones on diagonal

% Gaussian elimination
for j = 1:n-1
  for i = j+1:n
    L(i,j) = A(i,j) / A(j,j);   % row multiplier
    A(i,j:n) = A(i,j:n) - L(i,j)*A(j,j:n);
  end
end

U = triu(A);
end

A = [2 0 4 3; -4 5 -7 -10; 1 15 2 -4.5; -2 0 2 -13];
[L,U] = lufact(A)

LtimesU = L*U

A - LtimesU

b = [4;9;29;40];
z = forwardsub(L,b);
x = backsub(U,z)

b - A*x

n = 6;
A = magic(n);
x = ones(n,1);
y = zeros(n,1);
for i = 1:n
    for j = 1:n
        y(i) = y(i) + A(i,j)*x(j);   % 2 flops
    end
end

n_ = (400:400:4000)';
t_ = 0*n_;
for i = 1:length(n_)
    n = n_(i);
    A = randn(n,n);  x = randn(n,1);
    tic  % start a timer
    for j = 1:10   % repeat ten times
        A*x;
    end
    t = toc;       % read the timer
    t_(i) = t/10;  % seconds per instance
end

% table(n_,t_,'variablenames',{'size','time'})  % warning: the 'table' function is not yet implemented in Octave
for i = 1:length(n_)
    display([num2str(n_(i)),'   ' ,num2str(t_(i))])
end

n_ = (400:200:6000)';
t_ = 0*n_;
for i = 1:length(n_)
    n = n_(i);
    A = randn(n,n);  x = randn(n,1);
    tic  % start a timer
    for j = 1:10  % repeat ten times
        A*x;
    end
    t = toc;  % read the timer
    t_(i) = t/10;  
end

loglog(n_,t_,'.-')
xlabel('size of matrix'), ylabel('time (sec)')
title('Timing of matrix-vector multiplications')

loglog(n_,t_,'.-')
xlabel('size of matrix'), ylabel('time (sec)')
title('Timing of matrix-vector multiplications')
hold on, loglog(n_,t_(1)*(n_/n_(1)).^2,'--')
axis tight
legend('data','O(n^2)','location','southeast')

n_ = (200:100:2400)';
t_ = 0*n_;
for i = 1:length(n_)
    n = n_(i);
    A = randn(n,n);  
    tic  % start a timer
    for j = 1:6,  [L,U] = lu(A);  end
    t = toc;  % read the timer
    t_(i) = t/6;  
end

loglog(n_,t_,'.-')
hold on, loglog(n_,t_(end)*(n_/n_(end)).^3,'--')
axis tight
xlabel('size of matrix'), ylabel('time (sec)')
title('Timing of LU factorization')
legend('lu','O(n^3)','location','southeast')

A = [2 0 4 3 ; -4 5 -7 -10 ; 1 15 2 -4.5 ; -2 0 2 -13]
b = [ 4; 9; 29; 40 ]

function x = forwardsub(L,b)
% FORWARDSUB   Solve a lower triangular linear system.
% Input:
%   L    lower triangular square matrix (n by n)
%   b    right-hand side vector (n by 1)   
% Output:
%   x    solution of Lx=b (n by 1 vector)

n = length(L);
x = zeros(n,1);
for i = 1:n
  x(i) = ( b(i) - L(i,1:i-1)*x(1:i-1) ) / L(i,i);
end

end

[L,U] = lufact(A);
x = backsub( U, forwardsub(L,b) )

A([2 4],:) = A([4 2],:);  
b([2 4]) = b([4 2]); 
x = A\b

[L,U] = lufact(A);
L

A = [ 2 0 4 3; -2 0 2 -13 ; 1 15 2 -4.5 ; -4 5 -7 -10 ];
b = [ 4; 40; 29; 9 ];

[L,U,P] = lu(A)

x = backsub( U, forwardsub(L,P*b) )

[PtL,U] = lu(A)

x = U \ (PtL\b)

x = A\b

x = [2;-3;1;-1];
twonorm = norm(x)         % or norm(x,2)
infnorm = norm(x,inf)
onenorm = norm(x,1)

A = [ 2  0; 1  -1 ]

twonorm = norm(A)

onenorm = norm(A,1)

max( sum(abs(A),1) )   % sum along rows (1st matrix dimension)

infnorm = norm(A,inf)
max( sum(abs(A),2) )   % sum along columns (2nd matrix dimension)

theta = linspace(0,2*pi,601);
x = [ cos(theta); sin(theta) ];  % 601 unit columns
subplot(1,2,1), plot(x(1,:),x(2,:)), axis equal
title('Unit circle in 2-norm')
xlabel('x_1'),  ylabel('x_2')

Ax = A*x;

subplot(1,2,2), plot(Ax(1,:),Ax(2,:)), axis equal
hold on, plot(twonorm*x(1,:),twonorm*x(2,:),'--')
title('Image of Ax, with ||A||')
xlabel('x_1'),  ylabel('x_2')

A = hilb(7);
kappa = cond(A)

x_exact = (1:7)';
b = A*x_exact;

randn('state',333);          % reproducible results 
dA = randn(size(A));  dA = 1e-12*(dA/norm(dA));
db = randn(size(b));  db = 1e-12*(db/norm(db));

x = (A+dA) \ (b+db); 
dx = x - x_exact;

rel_error = norm(dx) / norm(x_exact)

b_bound = kappa * 1e-12/norm(b)
A_bound = kappa * 1e-12/norm(A)

x = A\b;
rel_error = norm(x - x_exact) / norm(x_exact)
rounding_bound = kappa*eps

A = hilb(14);
kappa = cond(A)

rounding_bound = kappa*eps

x_exact = (1:14)';
b = A*x_exact;  x = A\b;

relative_error = norm(x_exact - x) / norm(x_exact)

A = [ 2 -1  0  0  0  0
      4  2 -1  0  0  0
      0  3  0 -1  0  0
      0  0  2  2 -1  0
      0  0  0  1  1 -1
      0  0  0  0  0  2 ]

diag_main = diag(A,0)'

diag_plusone = diag(A,1)'

diag_minusone = diag(A,-1)'

A = A + diag([5 8 6 7],2)

[L,U] = lufact(A);
subplot(1,2,1), spy(L), title('L')
subplot(1,2,2), spy(U), title('U')

[L,U,P] = lu(A);
subplot(1,2,1), spy(L), title('L')
subplot(1,2,2), spy(U), title('U')

n = 8000;
A = diag(1:n) + diag(n-1:-1:1,1) + diag(ones(n-1,1),-1);

tic, [L,U] = lu(A); toc

tic, [L,U] = lu(sparse(A)); toc

A = [  2     4     4     2
       4     5     8    -5
       4     8     6     2
       2    -5     2   -26 ];

L1 = eye(4); L1(2:4,1) = [-2;-2;-1];
A1 = L1*A

A2 = (L1*A1')'

A2 = A1*L1'

L2 = eye(4);  L2(3:4,2) = [0;-3]; 
A3 = L2*A2*L2'

L3 = eye(4); L3(4,3) = -1;
D = L3*A3*L3'

A = magic(4) + eye(4);
B = A+A'

R = chol(B) 

B = A'*A

R = chol(B)

norm( R'*R - B )

chol( triu(B) )
