% % % % % % % % % % % % % % % % % % % % % % % %
% Loesung der Programmieraufgabe A11, WS06/07 %
% % % % % % % % % % % % % % % % % % % % % % % %

clear all
close all
% vorab)
A = [4 0 -2 -2; 0 1 3 0; -2 3 11 3; -2 0 3 14];
b = [-10; 11; 49; 63];

[L,x]=cholesky(A,b);

disp('Loesung Aufgabe 11 - Testmatrix aus Aufgabe 9:')
disp(' ')
disp('Die linke untere Dreiecksmatrix L ist L=');
disp(L);
disp('Loesungsvektor x=');
disp(x);
disp(' ')


% ===============================
% a)  f(t)=1, y_a=0, y_b=50, n=9:
% ===============================
a=0;                    % linke Intervallgrenze
b=10;                   % rechte Intervallgrenze
n=9;                    % Anzahl der Gitterpunkte im Intervall [a,b]
ya=0;                   % linke Randbedingung
yb=50;                  % rechte Randbedingung
h = (b-a)/(n+1);        % Schrittweite zwischen zwei Gitterpunkten
% Bilde Matrix A: 
% ---------------

% langsame Variante: 
A=zeros(n,n);
for j=1:n
    A(j,j)=2;           % Hauptdiagonale
    if j<n
        A(j,j+1)=-1;    % rechte Nebendiagonale
        A(j+1,j)=-1;    % linke Nebendiagonale
    end
end

% schnelle Variante: 
% A ist eine duennbesetzte Matrix. Matlab kann effizient mit duennbesetzten
% Matrizen rechnen, das spart (gerade bei großen Problemen) sehr viel
% Rechenzeit (ca. doppelt so schnell)
A=spalloc(n,n,3*n);     % A ist duennbesetzte Matrix, n Zeilen, n Spalten 
                        % maximal 3n Elemente, die von Null verschieden
                        % sind
Adiag=2*ones(n,1);      % Hauptdiagonale
Adiag1=-1*ones(n,1);    % Nebendiagonalen
A=spdiags([Adiag1 Adiag Adiag1],-1:1,n,n); % generiert Tridiagonalmatrix

% Bilde rechte Seite b_vec: 
% -------------------------
b_vec=-h^2*ones(n,1);   % innere Punkte von b_vec
b_vec(1)=b_vec(1)+ya;   % linke Randbedingung zufuegen 
b_vec(n)=b_vec(n)+yb;   % rechte Randbedingung zufuegen 

[L,y]=cholesky(A,b_vec);% Loese LGS

y_exakt=zeros(n,1);     % Initialisierung der exakten Loesung
for j=1:n
    y_exakt(j)=0.5*(a+h*j)^2;   % exakte Loesung, zum Vergleich  
end

y_plot=zeros(n+2,1);    % Randbedingungen sind fest und werden nicht mit-
                        % berechnet. Daher muessen sie als extra Punkte
                        % hinzugefuegt werden, baue einen neuen Vektor, in
                        % den ich die Randpunkte einfuegen kann.
y_plot(2:n+1)=y;        % innere Punkte der Loesung
y_plot(1)=ya;           % linke Randbedingung 
y_plot(n+2)=yb;         % rechte Randbedingung

y_exakt_plot=zeros(n+2,1);  % analog zu y_plot, hier fuer die exakte Lsg.
y_exakt_plot(2:n+1)=y_exakt;
y_exakt_plot(1)=ya;
y_exakt_plot(n+2)=yb;

err=y_plot-y_exakt_plot;    % absoluter Fehler der numerischen Loesung

% grafische Ausgabe: 
% ------------------

figure(1)      
plot(a:h:b,y_plot')
xlabel('t')
ylabel('y(t)')
title('Loesung von y^{(2)}(t)=1, n=9 Gitterpunkte')

figure(2)
plot(a:h:b,err)
title('Abs. Fehler e(t) der num. Lsg. von y^{(2)}(t)=1, n=9')
xlabel('t')
ylabel('e(t)')


% bessere grafische Ausgabe:
% (nur besser fuer Ausdruck, fuer Uebungsgruppen irrelevant...)
scrsz = get(0,'ScreenSize');
figure('Position',[200 600 900 400])
subplot(1,2,1)
plot(a:h:b,y_plot)
axis square
xlabel('t')
ylabel('y(t)')
title(['Loesung von y ^{(2)}(t)=1, n=',num2str(n),' Gitterpunkte'])
subplot(1,2,2)
plot(a:h:b,err,'-')
axis square
title('Absoluter Fehler e(t)')
xlabel('t')
ylabel('e(t)')
print -deps A12_a.eps %Export als EPS


% =============================================
% b)  f(t)=3t^2, y_a=0, y_b=2500, n=9,19,49,99:
% =============================================


a=0;                    % linke Intervallgrenze
b=10;                   % rechte Intervallgrenze
n_vec=[9;19;49;99];     % Anzahl der Gitterpunkte im Intervall [a,b]
ya=0;                   % linke Randbedingung
yb=2500;                % rechte Randbedingung


for k=1:4               % vier verschiedene n's
    n=n_vec(k);         % Wahl von n
    h = (b-a)/(n+1);    % Schrittweite zwischen zwei Gitterpunkten
    
    % schnelle Variante: 
    % A ist eine duennbesetzte Matrix. Matlab kann effizient mit duennbesetzten
    % Matrizen rechnen, das spart (gerade bei großen Problemen) sehr viel
    % Rechenzeit!
    A=spalloc(n,n,3*n);     % A ist duennbesetzte Matrix, n Zeilen, n Spalten 
                            % maximal 3n Elemente, die von Null verschieden
                            % sind
    Adiag=2*ones(n,1);      % Hauptdiagonale
    Adiag1=-1*ones(n,1);    % Nebendiagonalen
    A=spdiags([Adiag1 Adiag Adiag1],-1:1,n,n); % generiert Tridiagonalmatrix

    % Bilde rechte Seite b_vec: 
    % -------------------------
    b_vec=zeros(n,1);           % Initialisieren
    for j=1:n
        b_vec(j)=3*(a+j*h)^2;   % innere Punkte von b_vec
    end
    b_vec=-h^2*b_vec;
    b_vec(1)=b_vec(1)+ya;       % linke Randbedingung zufuegen 
    b_vec(n)=b_vec(n)+yb;       % rechte Randbedingung zufuegen 

    [L,y]=cholesky(A,b_vec);    % Loese LGS

    y_exakt=zeros(n,1);         % Initialisierung der exakten Loesung
    for j=1:n
        y_exakt(j)=0.25*(a+h*j)^4;   % exakte Loesung, zum Vergleich  
    end

    y_plot=zeros(n+2,1);    % Randbedingungen sind fest und werden nicht mit-
                            % berechnet. Daher muessen sie als extra Punkte
                            % hinzugefuegt werden, baue einen neuen Vektor, in
                            % den ich die Randpunkte einfuegen kann.
    y_plot(2:n+1)=y;        % innere Punkte der Loesung
    y_plot(1)=ya;           % linke Randbedingung 
    y_plot(n+2)=yb;         % rechte Randbedingung

    y_exakt_plot=zeros(n+2,1);  % analog zu y_plot, hier fuer die exakte Lsg.
    y_exakt_plot(2:n+1)=y_exakt;
    y_exakt_plot(1)=ya;
    y_exakt_plot(n+2)=yb;

    err=y_plot-y_exakt_plot;  % absoluter Fehler der numerischen Loesung
   
    % Speicher die Fehler und Schrittweiten fuer einen spaeteren Plot
    if k==1
        h1=h;
        err1=err;
    elseif k==2
        h2=h;
        err2=err;
    elseif k==3
        h3=h;
        err3=err;
    elseif k==4
        h4=h;
        err4=err;
    end

end

% grafische Ausgabe: 
% ------------------

figure(4)      
plot(a:h:b,y_plot)
xlabel('t')
ylabel('y(t)')
title(['Loesung von y ^{(2)}(t)=3t^2, n=,',num2str(n),' Gitterpunkte'])
figure(5)
plot(a:h1:b,err1,'-')
hold on
plot(a:h2:b,err2,'--')
plot(a:h3:b,err3,'-.')
plot(a:h4:b,err4,':')
hold off
title(['Absoluter Fehler e(t) der numerischen Loesung von  y ^{(2)}(t)=3t^2'])
legend('n=9','n=19','n=49','n=99',2)
xlabel('t')
ylabel('e(t)')


% bessere grafische Ausgabe:
% (nur besser fuer Ausdruck, fuer Uebungsgruppen irrelevant...)
figure('Position',[200 600 900 400])
subplot(1,2,1)
plot(a:h:b,y_plot)
axis square
xlabel('t')
ylabel('y(t)')
title(['Loesung von y ^{(2)}(t)=3t^2, n=',num2str(n),' Gitterpunkte'])
subplot(1,2,2)
plot(a:h1:b,err1,'-')
hold on
plot(a:h2:b,err2,'--')
plot(a:h3:b,err3,'-.')
plot(a:h4:b,err4,':')
hold off
axis square
title(['Absoluter Fehler e(t) der numerischen Loesung von y^{(2)}(t)=3t^2'])
legend('n=9','n=19','n=49','n=99',2)
xlabel('t')
ylabel('e(t)')
print -deps A12_b.eps %Export als EPS
