%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SEMESTRALNI PRACE MAT4
% Viktor Bubla
% FM, 2. rocnik, EIRS
% 27.4.2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Zadana matice:    dorr
% Zadana metoda:    Jacobi
% Zadana rovnice:   y' - x*(y^2 + 1) = 0
%                   y' -xy^2 -x = 0
% Poc. podminka:    y(0) = -1
% Zadany interval:  (0 ; 0.3)
% Hodnoty v bodech: 0.1, 0.2, 0.3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vymazani vsech promennych, zavreni oken
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vytvoreni matice
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vel = 20;                               % Velikost matice
A = gallery('dorr', vel, 0.1);          % Vygenerovani matice

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reseni Ax = b pomoci /
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x(1:vel,1) = 1;                         % (Jednotkovy) vektor x
b = A * x;                              % Vypocet prave strany
reseniA = A\b;                          % Primy vypocet reseni

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Jacobiova iteracni metoda
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = 2000;                               % Maximalni pocet iteraci

D = diag(diag(A));                      % Diagonala matice A
L = tril(-A,-1);                        % Dolni trojuhelnikova matice z A
U = triu(-A,1);                         % Horni trojuhelnikova matice z A

Tj = inv(D)*(L+U);                      % Jacobiho iteracni matice
cj = inv(D)*b;                          % Jacobiho konstantni vektor

tol = 1e-07;                            % Tolerance reseni
k = 0;                                  % Promenna cyklu
x = zeros(vel,1);                       % Pocatecni odhad

while k <= N
   k = k+1;
   x(:,k+1) = Tj*x(:,k) + cj;
   if norm(x(:,k+1)-x(:,k)) < tol
      disp('Vypocet byl uspesny!')
      str = sprintf('Podminka ||x^(k+1) - x^(k)|| < %d byla splnena po %d iteracich.', tol, k);
      disp(str);
      reseniB = x(:,k+1);
      disp('Reseni vypoctem x = A \ b');
      disp(reseniA);
      disp('Reseni iteracni metodou');
      disp(reseniB);
      break
   end
end

if norm(x(:,k+1) - x(:,k)) > tol || k > N
   str = sprintf('Byl dosazen maximalni pocet iteraci aniz by byla splnena podminka ||x^(k+1) - x^(k)|| < %f', tol);
   disp('Zkontrolujte sekvenci iteraci.')
   disp('V pripade, ze zjistite konvergenci, zvyste maximalni pocet iteraci.')
   disp('V pripade divergence je mozne, ze matice A neni diagonalne dominantni.')
   disp(x);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nalezeni vlastnich cisel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
A = gallery('dorr', 5, 0.1);            % Vygenerovani matice
vl_cisla = eigs(A);                     % Ziskani vlastnich cisel
% Poznamka: Je potreba misto eig pouzit funkci eigs, protoze matice je
% "sparse" (ridka).
disp('Vlastni cisla matice dorr:');
disp(vl_cisla);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nalezeni singularnich cisel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sing_cisla = svds(A);                   % Ziskani singularnich cisel
% Poznamka: Stejne jako u vlastnich cisel pouzijeme misto svd funkci svds.
disp('Singularni cisla matice dorr:');
disp(sing_cisla);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reseni pocatecni ulohy pro diferencialni rovnici
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

f = @(x,y) x*(y^2+1);                   % Definice diferencialni rovnice

x = 0 : 0.001 : 0.3;                    % Casove body reseni

[Xo45, Yo45] = ode45(f, x, -1);         % Reseni funkci ode45
[Xo113, Yo113] = ode113(f, x, -1);      % Reseni funkci ode113

plot(Xo45, Yo45, Xo113, Yo113);         % Srovnani vypoctenych prubehu

disp('Reseni v bode 0.1 obema funkcemi:');
disp(Yo45(101));
disp(Yo113(101));
disp('Reseni v bode 0.2 obema funkcemi:');
disp(Yo45(201));
disp(Yo113(201));
disp('Reseni v bode 0.3 obema funkcemi:');
disp(Yo45(301));
disp(Yo113(301));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vystup programu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PRVNI CAST

% Vypocet byl uspesny!
% Podminka ||x^(k+1) - x^(k)|| < 1.000000e-07 byla splnena po 1859 iteracich.
% Reseni vypoctem x = A \ b
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
% 
% Reseni iteracni metodou
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%     1.0000
%
% Vlastni cisla matice dorr:
%    15.4451
%    13.2177
%     8.4390
%     4.1823
%     0.7159
% 
% Singularni cisla matice dorr:
%    15.5867
%    13.3844
%     8.5347
%     4.1302
%     0.7014

% DRUHA CAST

% Reseni v bode 0.1 obema funkcemi:
%    -0.9900
% 
%    -0.9900
% 
% Reseni v bode 0.2 obema funkcemi:
%    -0.9608
% 
%    -0.9608
% 
% Reseni v bode 0.3 obema funkcemi:
%    -0.9138
% 
%    -0.9138
%