// Exercice 1
//
A=[10   7    8    7
   7    5    6    5
   8    6   10    9
   7    5    9   10]
b=[32; 23; 33; 31]
x=A\b

[U, S, V]=svd(A)
cond(A)
db=1/5*V(:,4) 

xdx=A\(b+db)
dx=xdx-x
(norm(dx)/norm(x))/(norm(db)/norm(b))


// Exercice 2
//
W=triu(rand(50, 50));
A=W'*W;
R=chol(A);
R1=R+1e-6*triu(rand(A));
norm(W-R)/norm(W)
norm(W-R1)/norm(W)
norm(A-R'*R)/norm(A)
norm(A-R1'*R1)/norm(A)



// Exercice 4
//
m=100; n=15;
t=(0:m-1)'/(m-1); A=[]; 
for i=1:n,  A = [A t.^(i-1)]; end
b=exp(sin(4*t)); b=b/2006.7874531048527;

format(20)
x1=(A'*A)\(A'*b)
[Q,R]=qr(A); x2=R\(Q'*b)
[U,S,V]=svd(A); x3=V*(S\(U'*b))
x4=A\b
plot2d(t, [A*x1 A*x2 A*x3 A*x4])     

//Exercice 4
//
tmp = read('temp.dat', 145, 2); 
years=tmp(:,1)-1856; temp=tmp(:,2);
plot2d(years, temp, style=-1)
A2=[ones(temp)  years];                  
x2=A2\temp                               
plot2d(years, A2*x2, style=2, frameflag=0)
xset('window', 1)
plot2d(years, A2*x2-temp, style=2)        

A6=[];   
for i=1:6 ,  A6 = [A6 years.^(i-1)]; end    
x6=A6\temp
xset('window', 0)
plot2d(years, A6*x6, style=3, frameflag=0)
xset('window',1)
plot2d(years, A6*x6-temp, style=3, frameflag=0)