%% H2T18 Matriisin ratkaisujen virheherkkyys, häiriöalttius  

%Ylempi yhtälösysteemi muodossa Ax=b :

% kerroinmatriisi:
A1=[ 10 7 8 7;
    7 5 6 5; 
    7 6 10 9;
    7 5 9 10
    ];
% oikean puolen tekijät pystyvektorina:
b1=[32 23 33 31]';

%   Alempi vastaavasti muodossa A2x=b2 :
format compact
A2=[2 1 5 1;
    1 1 -3 -1;
    3 6 -2 1;
    2 2 2 -3
    ];
b2=[9 -5 8 3]';

%   a)
%   ratkaisut näille:
x1=A1\b1;
x2=A2\b2;
x1t=x1'    % Tilan säästämiseksi 
x2t=x2'    %  ... transponoidaan 
%{
x1t =
   -0.1111    2.8889    0.4444    1.3333
x2t =
   -1.5000    2.5208    1.7292    0.8333
%}
%%  b)

c1=[32.1 22.9 32.9 31.1]'
c2=[32.01 22.99 32.99 31.01]'

u1=A1\c1
u2=A1\c2

[x1 u1 u2]

%%
% pannaan merkille että ensimmäinen muutos vaikuttaa tulokseen melko
% paljon, etumerkitkin muuttuvat

d1=[9.1 -5.1 7.9 3.1]'
d2=[9.01 -5.01 7.99 3.01]'
v1=A2\d1
v2=A2\d2
%   huomataan, että matriisi A2 ei ole läheskään yhtä herkkä datan
%   muutoksille
[x2 v1 v2]
%%  c)
%   lisätään satunnaisuutta
Ar1 = A1+0.1*rand(4,4);
Ar2= A2+0.1*rand(4,4);
%   ja ratkaistaan alkuperäisillä b:illä
Ar1\b1
Ar2\b2

%%  d)  Suhteellisen virheen kasvu ja cond(A)
%{
Lasketaan "SVS"-epäyhtälön molemmat puolet matriiseilla A1 ja A2,
kun b=b1 ja muutettu b=c1, ts delb=c1-b1
%} 

k = cond(A1)
b=b1
delb=c1-b1
x=x1
delx=u1-x1

ratksuhtvirhe=norm(delx)/norm(x)
datansuhtvirhe=norm(delb)/norm(b)
%{
ratksuhtvirhe =
    0.5342
datansuhtvirhe =
    0.0033
%}
k=cond(A1)
ratksuhtvirhe/datansuhtvirhe
%{
k =
  455.4940
ans =
  160.3286
%}
%{
Voidaan sanoa, että suuruusluokka on oikea, pahin tapaus, jota häiriöluku
edustaa, on noin 3-kertainen todelliseen suht. virheeseen nähden ja molemmat
ovat "satoja".
%}
%%
%
k2 = cond(A2)
%{
k2 =
    9.2491
%}

%   Ensimmäinen matriisi on siis huomattavasti häiriöalttiimpi
%   Jo tarkastelemalla ensimmäisen matriisin kertoimia, voidaan huomata
%   näiden olevan lähellä lineaarista riippuvuutta.

% Tässä voidaan tehdä aivan vastaavat laskut ja todeta huomattavasti
% parempi käytös (ja epäilemättä k2:n luokkaa oleva virhesuhde).
% Jätetään yksityiskohdat oppilaalle itselleen.

% publish('H2T18R','pdf')