clear;
//// ***************************
//
//     Procedure SI ed SA
//
// ***************************
//
function x = SI(T,c)
//
// Funzione SI: per T matrice N x N triangolare superiore invertibile e c colonna
// di N componenti, calcola un'approssimazione della colonna soluzione del sistema
// di equazioni lineari
//
//       T x = c
//
// x: se T invertibile, la colonna x (di N numeri reali) approssima la soluzione
//    del sistema.
//
N = size(T,'r');
//
// *** Verifica se T invertibile. 
//
if find(diag(T)==0) ~= [] then
   error("SI: elementi nulli sulla diagonale.");
   else//
       // *** Calcolo soluzione
       //
       x = zeros(N,1);
       //
       x(N) = c(N) / T(N,N);
       for k = N-1:-1:1,
           x(k) = (c(k) - T(k,k+1:N)*x(k+1:N)) / T(k,k);
       end;
end;
endfunction
//
function x = SA(T,c)
//
// Funzione SA: per T matrice N x N triangolare inferiore invertibile e c colonna
// di N componenti, calcola un'approssimazione della colonna soluzione del sistema
// di equazioni lineari
//
//       T x = c
//
// x: se T invertibile, la colonna x (di N numeri reali) approssima la soluzione
//    del sistema.
//
N = size(T, 'r');
//
// *** Verifica se T invertibile. 
//
if find(diag(T)==0) ~= [] then
   error("SA: elementi nulli sulla diagonale.");
   else//
       // *** Calcolo soluzione
       //
       x = zeros(N,1);
       //
       x(1) = c(1) / T(1,1);
       for k = 2:N,
           x(k) = (c(k) - T(k,1:k-1)*x(1:k-1)) / T(k,k);
       end;
end
endfunction
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A = [2,3;3,5];
b = [5;6]; // [5;8]; [1;-1] 
b = b / norm(b,%inf);
//
[S,D,P] = lu(A);
c = SA(S,P*b);
x = SI(D,c);
//
Eps_b = 1d-1;
r_b = Eps_b * norm(b,%inf);
Eps_x = cond(A,%inf) * Eps_b;
r_x = Eps_x * norm(x,%inf);
//
tol = max(r_b,r_x) * 1.2;
xmin = min(b(1),x(1)) - tol;
ymin = min(b(2),x(2)) - tol;
xmax = max(b(1),x(1)) + tol;
ymax = max(b(2),x(2)) + tol;
R = [xmin,ymin,xmax,ymax];
//
N = 7;
db1 = r_b * linspace(-1,1,N)';
db2 = r_b * linspace(-1,1,N)';
//
clf();
subplot(121);
xtitle('b + db');
c = b;
r = r_b;
Ix = [c(1)-r,c(1)-r,c(1)+r,c(1)+r];
Iy = [c(2)-r,c(2)+r,c(2)+r,c(2)-r];
xpoly(Ix,Iy,'lines',1);
q=gce();
q.foreground=5;
subplot(122);
xtitle('x + dx');
c = x;
r = r_x;
Ix = [c(1)-r,c(1)-r,c(1)+r,c(1)+r];
Iy = [c(2)-r,c(2)+r,c(2)+r,c(2)-r];
xpoly(Ix,Iy,'lines',1);
q=gce();
q.foreground=5;
//
for k = 1:N,
    for j = 1:N,
        db = [db1(k);db2(j)];
        subplot(121)
        plot2d(b(1)+db(1),b(2)+db(2), style = -2,rect = R,frameflag = 3);
        c = SA(S,P*db);
        dx = SI(D,c);
        subplot(122)
        plot2d(x(1)+dx(1),x(2)+dx(2), style = -1,rect = R,frameflag = 3);
    end
end
//
subplot(121);
plot2d(b(1),b(2),style = -5,rect = R);
xgrid(2);
subplot(122);
plot2d(x(1),x(2),style = -5,rect = R);
xgrid(2);
