%% Verschiedene Moeglichkeiten, die Differenzenmatrix der zweiten Ableitung zu erzeugen
% CWB SS 2014

%% Eindimensional: die (-1, 2, -1)-Tridiagonalmatrix
% Setze Matrixdimension

n = 8;
A = zeros(n);

%% Doppelschleife
%
for i=1:n
    for j=1:n
        if(i==j) 
            A(i,j)=2;
        elseif(abs(i-j)==1)
            A(i,j)=-1;
        else
            A(i,j)=0;
        end
    end
end

%% Einfache Schleife
%
A = zeros(n);
for i=1:n
    A(i,i) = 2; % Diagonalelemente
    if(i<n)
        A(i,i+1) = -1;
        A(i+1,i) = -1;
    end
end

%% Einfache Schleifen, ohne Verzweigungen
%

A = zeros(n);
A(1,1) = 2;
for i=2:n
    A(i,i) = 2;
    A(i-1,i) = -1;
    A(i,i-1) = -1;
end


%% Indexvektor
A = zeros(n);
A(1,1:2) = [2, -1];
for i=2:n-1
    A(i,i-1:i+1) = [-1, 2, -1];
end
A(n,n-1:n) = [-1, 2];

%% Sparse Matrix
%
% |Indizes und Werte der Hauptdiagonale
i = 1:n;
j = 1:n;
s = 2*ones(size(i));

% Indizes und Werte der Ost-Nachbarn kommen dazu
i = [i, 1:n-1] ; 
j = [j, 2:n];
s = [s, -ones(1,n-1)];

% Indizes und Werte der West-Nachbarn kommen dazu
i = [i, 2:n];  
j = [j, 1:n-1];
s = [s, -ones(1,n-1)];

A = sparse(i,j,s);

%% Befehl diag
e = ones(n,1);
A = diag(2*e) - diag(e(2:end),1)- diag(e(2:end),-1) ;

%% Befehl spdiags
e = ones(n,1);
A = spdiags([-e, 2*e, -e], -1:1, n, n);
 
%% Zweidimensional: die (-1, -1, 4, -1, -1)-Pentdiagonalmatrix
% Setze Gitter- und Matrixdimensionen
nx = 3;
ny = nx; 
n = nx*ny; % Anzahl der Unbekannten

e = ones(nx,1);
D = spdiags([-e, 4*e, -e], -1:1, nx, nx);
% Diagonale Block-Matrix

E = -speye(ny); % Nebendiagonalen-Bloecke

A = spalloc(n,n,5*n); % Platz zuweisen.
for i=1:ny
    i0 = (i-1)*nx + 1;   % Start- Index des Diagonal-blocks
    i1 =  i*nx       ; % End- Index des Diagonal-blocks 
    A(i0:i1,i0:i1) = D;
end

for i=1:ny-1
    i0 = (i-1)*nx + 1;   % Start- Index des Diagonal-blocks
    i1 =  i*nx       ; % End- Index des Diagonal-blocks 
    A(i0:i1,(i0+nx):(i1+nx)) = E;
end







%% Als Kronecker-Produkt
nx = 3;
ny = 4; 
n = nx*ny;

e = ones(nx,1);
D = spdiags([-e, 2*e, -e], -1:1, nx, nx);
E = speye(ny);

% kron(E,D) ersetzt Einsen in E durch D-Bloecke
% kron(D,E) ersetzt Eintraege in D durch d_ij mal Einheitsmatrizen: gibt
% Hauptdiagonalbeitrag 2 und Aussendiagonalen -1
A = kron(D,E) + kron(E,D);




