clear
t = (1955:5:2000)';
y = [ -0.0480; -0.0180; -0.0360; -0.0120; -0.0040;
    0.1180; 0.2100; 0.3320; 0.3340; 0.4560 ];
plot(t,y,'.')
xlabel('year'), ylabel('anomaly ({\circ}C)'), axis tight   % ignore this line

t = (t-1950)/10;  
V = t.^0;      % vector of ones
n = length(t);
for j = 1:n-1
    V(:,j+1) = t.*V(:,j);
end
c = V\y;
p = @(x) polyval(flipud(c),(x-1950)/10);
hold on, fplot(p,[1955 2000])

year = (1955:5:2000)';
y = [ -0.0480; -0.0180; -0.0360; -0.0120; -0.0040;
    0.1180; 0.2100; 0.3320; 0.3340; 0.4560 ];

t = year - 1955;             % better matrix conditioning later
V = [ t.^0 t ];    % Vandermonde-ish matrix
c = V\y;
f = @(x) polyval(c(end:-1:1),x-1955);
fplot(f,[1955 2000])
hold on, plot(year,y,'.')
xlabel('year'), ylabel('anomaly ({\circ}C)'), axis tight   % ignore this line

V = [ t.^0 t t.^2 t.^3];    % Vandermonde-ish matrix
c = V\y;   f = @(x) polyval(c(end:-1:1),x-1955);
fplot(f,[1955 2000])

k = (1:100)';
s = cumsum( 1./k.^2 );   % cumulative summation
p = sqrt(6*s);
plot(k,p,'.-')
xlabel('k'), ylabel('p_k')
title('Sequence converging to \pi')

ep = abs(pi-p);    % error sequence
loglog(k,ep,'.'), title('log-log convergence')
xlabel('k'), ylabel('error'),  axis tight

V = [ k.^0, log(k) ];    % fitting matrix
c = V \ log(ep)          % coefficients of linear fit

a = exp(c(1)),  b = c(2),

hold on, loglog(k,a*k.^b,'r'), title('power-law fit')
axis tight    % ignore this line

t = linspace(0,3,400)';
A = [ sin(t).^2, cos((1+1e-7)*t).^2, t.^0 ];
kappa = cond(A)

x = [1;2;1];
b = A*x;

x_BS = A\b;
observed_err = norm(x_BS-x)/norm(x)
max_err = kappa*eps

N = A'*A;
x_NE = N\(A'*b);
observed_err = norm(x_NE-x)/norm(x)
digits = -log10(observed_err)

A = magic(5);
A = A(:,1:4);
[m,n] = size(A)

[Q,R] = qr(A);
szQ = size(Q)
szR = size(R)

QTQ = Q'*Q
norm(QTQ - eye(m))

[Q,R] = qr(A,0);
szQ = size(Q)
szR = size(R)

Q'*Q - eye(4)

A = magic(6);
A = A(:,1:4);
[m,n] = size(A);

z = A(:,1);

v = z - norm(z)*eye(m,1);
P = eye(m) - 2/(v'*v)*(v*v');   % reflector

P*z

A = P*A

A(2:m,2:n)

z = A(2:m,2);
v = z - norm(z)*eye(m-1,1);
P = eye(m-1) - 2/(v'*v)*(v*v');   % reflector

A(2:m,2:n) = P*A(2:m,2:n)

for j = 3:n
    z = A(j:m,j);
    k = m-j+1;
    v = z - norm(z)*eye(k,1);
    P = eye(k) - 2/(v'*v)*(v*v');
    A(j:m,j:n) = P*A(j:m,j:n);
end

R = A
