%% Második gyakorlat
%% Normák, kondíciószám
%
% 1. feladat: Legyen x = (10^{200},1). Számítsuk ki x normáját az alábbi két módon:
% \Vert x \Vert = \sqrt{x_1^2 + x_2^2}
% \Vert x \Vert = c \cdot \sqrt{\left( \frac{x_1}{c}\right)^2 + \left( \frac{x_2}{c}\right)^2 }, ahol c= \mathrm{max} \left\{ |x_1|,|x_2| \right\}.

x=[10^200, 1];

x_norma1 = sqrt(x(1)^2 + x(2)^2)


c=max(abs(x(1)), abs(x(2)));
x_norma2 = c * sqrt((x(1)/c)^2 + (x(2)/c)^2)

% Látható, hogy csak az első esetben nem tudja kiszámolni a MATLAB az eredményt (mert 10^{400} már túl nagy érték).
%
% 2. feladat: A max és abs beépített függvények segítségével írjunk programot vektor_norma.m fájl néven, amely egy adott valós 𝑥 vektor felhasználó által kért 𝑝-normáját és maximumnormáját kiszámolja! Használjuk fel az előző feladat ötletét.

% függvénydef. a fájl alján
[norma_p, norma_max] = vektor_norma([1,-2,3], 1)

% Érdemes megjegyezni, hogy a MATLAB-ban létezik beépített norma számoló függvény is:

x=[-1 2 3]
norm(x,1)
norm(x,'Inf')

%
% 3. feladat: A switch-case szerkezet segítségével írjunk programot indukaltmatrixnorma.m fájl néven, mely adott négyzetes mátrixra és az 'egyes', 'kettes' vagy 'vegtelen' szavak input begépelése után kiszámolja az \Vert A  \Vert_1, \Vert A  \Vert_2 vagy \Vert A  \Vert_{\infty} normáját! Ha nem négyzetes mátrixot vagy a fenti sztringek valamilyék adjuk meg, akkor a programunk adjon vissza hibaüzenetet! Használjuk az előző feladatban említett 'norm' parancsot!

% függvénydef. a fájl alján
matrixnorma(eye(2), "egyes")

%
% 4. feladat: Tekintsük a félév végén tárgyalt, második (térbeli) derivált közelítésekor felbukkanó alábbi
% A_n = \left(\matrix{2 & -1 & 0 & \dots & 0\cr -1 & 2 & -1 & \dots & 0 \cr 0 & -1 & 2 & \dots & 0 \cr \vdots & \ddots & \ddots & \ddots& \vdots \cr 0 & \dots & 0 & -1 & 2}\right)
% mátrixot, valamint az alábbi, ún. Hilbert-mátrixot:
% H_n = \left(\matrix{1 & 1/2 & 1/3 & \dots & 1/n \cr 1/2 & 1/3 & 1/4 & \dots & 1/(n+1) \cr 1/3 & 1/4 & 1/5 & \dots & 1/(n+2) \cr \vdots & \vdots & \vdots & \vdots& \vdots \cr 1/n & 1/(n+1) & 1/(n+2) & \dots & 1/(2n-1)}\right)
% Írjunk programot, amely az alábbi két mátrix kettes kondíciószámát, valamint a condest értékét adja vissza adott n-re!
% Először tehát írjunk függvényt, ami adott n-re visszaadja az adott mátrixot,
% majd alkalmazzuk a megalkotott függvények kimenetére a cond, condest parancsokat.

% tfh. a dlap1d_matrix függvény szintén implementálva van (dlap1d_matrix.m)
% tfh. a hilbert_matrix függvény implementálva van (hilbert_matrix.m)

N=5;
% az első mátrix
A = dlap1d_matrix(N)

% a második mátrix
H = hilbert_matrix(N)

% Erre a mátrixra egyébként beépített parancs is van hilb néven (az input a dimenziója):
hilb(N)

% Használjuk a cond parancsot a kondíciószám kiszámítására!
% Alapesetben a cond(A) parancs is a 2-es normával számolt kondíciószámot
% adja vissza.
cond(A,2)
cond(H,2)

% Egy másik lehetőség, ha 1-es normában akarunk számolni, a condest
% parancs, ami becslést ad a kondíciószámra (ritka mátrixok esetén pontosabb).
cond(A,1)
cond(H,1)

condest(A)
condest(H)

%
% 5. feladat: Állítsuk elő azt az A \in \mathbb{R}^{100 \times 100} mátrixot, melnyek elemei az alábbi értékek:
% a_{ij}=\begin{cases}
% 1, \qquad \text{ha} \; i=j, \\
% -1, \qquad \text{ha} \; i<j, \\
% 0 \qquad \text{egyébként.}
% \end{cases}
% A lineáris egyenletrendszer jobboldali vektora b=(-98, -97, \dots, 0, 1)^{T}. Perturbáljuk a b
%  vektort úgy, hogy annak utolsó eleme 1 helyett 1.0001 legyen! Az előadáson tanultak fényében döntsük el, hogy a lineáris algebrai egyenletrendszer rosszul kondícionált-e, vagy sem, majd az alábbi formula használatával adjunk becslést a perturbáció következtében a megoldásban keletkezett hiba relatív mértékére.
% \mathrm{cond}_p(A) \frac{\Vert \delta b \Vert_p}{\Vert b \Vert_p} \geq \frac{\Vert \delta x\Vert_p}{\Vert x \Vert_p}
% Oldjuk is meg a feladatot a két különböző jobboldal esetén és nézzük meg, hogy mekkora a tényleges relatív hiba.

N=32; % Legyen először N kicsi, hogy lássunk is valamit

A = eye(N) + triu(-ones(N), 1) % Egy diagonális és egy felsőháromszög mátrix összege

b = (-N+2:1:1)';

delta_b=[zeros(N-1, 1); 
         0.0001]

size(A)
size(b)
size(delta_b)

[becsles_1, hiba_1] = becsles_es_hiba(A, b, delta_b, 1)
[becsles_2, hiba_2] = becsles_es_hiba(A, b, delta_b, 2)
[becsles_inf, hiba_inf]  = becsles_es_hiba(A, b, delta_b, 'inf')

%
% Szorgalmi feladat: Tekintsük az alábbi szabály szerint definiált mennyiséget, mely egy A mátrixhoz egy valós számot rendel:
% \Vert A \Vert_* = \sqrt[4]{\sum_{i=1}^n \sum_{j=1}^n a_{ij}^4 }
% (a) Bizonyítsuk be, hogy ez valóban norma!
% (b) Lehet ezt a normát vektornormából származtatni?
% (c) Írjunk programot, amely egy inputként kapott négyzetes A mátrix esetén kiszámítja a fenti mennyiséget úgy, hogy nem négyzetes mátrix esetén hibaüzenetet ad a program!
%
% A 2. feladat megoldása:

function [norma_p, norma_max] = vektor_norma(x,p)
    % x: a bemeneti vektor
    % p: a p-norma paramétere

    norma_max = max(abs(x));
    
    norma_p = norma_max * (sum(abs(x/norma_max).^p))^(1/p);
end

% A 3. feladat megoldása:

% Az indukaltmatrixnorma.m tartalma
function norma = matrixnorma(A,str)
    % A: bemeneti mátrix
    % str: a kért norma típusa

    % négyzetes a mátrix, és a str változó tipusa string
    %
    % utóbbi azért, mert MATLAB nem tud ebben a formában 
    % összehasonlitani string-et és double-t és
    % szokatlan hibaüzenetet ad ha valaki meghivja pl. 'str=2'-vel
    if size(A,1)==size(A,2) && isa(str, "string")
        switch str
            case 'egyes'
                norma=norm(A,1);
            case 'kettes'
                norma=norm(A,2);
            case 'vegtelen'
                norma=norm(A,'Inf');
            otherwise
                warning('Csak az "egyes", "kettes" és "vegtelen" szavak lehetnek a második input mezőben!')
        end
    
    else
        warning('Adj meg negyzetes matrixot!')
    
    end
end

% A 4. feladat függvényei:

function A = dlap1d_matrix(N)
    e = ones(N,1);
    A = spdiags([-e 2*e -e], [-1 0 1], N, N);
    A = full(A);
end

function H = hilbert_matrix(N)
    H = zeros(N);
    for i=1:N
        for j=1:N
            H(i,j) = 1/(i+j-1);
        end
    end
end

% Az 5. feladat függvénye(i):

function [becsles, hiba] = becsles_es_hiba(A, b, delta_b, p)    
    assert(any(b))
    becsles = cond(A,p)*norm(delta_b,p)/norm(b,p);
    hiba = norm(A\delta_b, p) / norm(A\b, p);
end