%Weiterbildung WS 1998/1999, Thema cg-Verfahren, Programm G. Opfer
%Es wird das Modellproblem delta u=0 auf [0,1]^2 bei gegebenen Randwerten
%mit dem cg-Verfahren gelöst. Gleichzeitig wird damit eine quadratische Form
%minimiert. Die Lösung u und alle weiteren Größen wie Residuum, Richtung haben
%die Form einer Matrix. Das zu lösende lin. Glsys. 
%                         r:=bb-Au=0 
%hat hier die Besonderheit, daß r selbst ohne Kenntnis von A 
%einfach ausgerechnet werden kann. 

%Es werden die folgenden Unterprogramme verwendet:
%r=residuum(x); %r=b-Ax
%v=prodd(h,b);  %Das Matrixprodukt Ah wird ausgerechnet unter der Voraussetzung,
%               %daß tatsächlich nur das Residuum r=b-Ax berechnet werden kann.
%x=scalar(A,B); %Zwei gleich große Matrizen A, B werden als Vektoren aufgefaßt
                %und skalar miteinander multipliziert.
%Man vgl. auch die dortigen Kommentare.

%=============================================================================
%Dimensionierung des Problems:
n=10;hh=1/(n-1);
%Das Verfahren kommt auch bei h=1/500 (294001 Unbekannte) nach 811 Schritten
%sicher zum Ende und braucht auf einer Sun4 1708 Sekunden also ca 28 Minuten.
%Das entstehende Bild hat aber viel zu viel Information und kann kaum 
%irgendwo gedruckt werden. 
%Bei h=1/50 (2401 Unbekannte) werden 82 Schritte und etwa 3.4'' benötigt.
unbek=(n-2)^2;
%=============================================================================
%0. Schritt, Vorbesetzung mit den Randwerten cos(6*pi*x)*cos(4*pi*y).
%Das Quadrat wird mit kleinen hh x hh Quadraten diskretisiert.
u=zeros(n,n);
x=linspace(0,1,n);y=linspace(0,1,n);
u(n,:)=cos(6*pi*x);   %I   x-Achse unten
u(:,n)=cos(4*pi*y');  %II  y-Achse rechts
u(1,:)=cos(6*pi*x);   %III x-Achse oben
u(:,1)=cos(4*pi*y');  %IV  y-Achse links
disp(['Rechnung für Maschenweite 1/h = ',num2str(round(1/hh)),': ',num2str(unbek),...
' Unbekannte']);
%=============================================================================
%1. Schritt, Initialisierung, erste Richtung h
r=residuum(u);
bb=r;h=r;      tic %Stoppuhr läuft
for j=1:unbek
%=============================================================================
%2. Schritt, Abfrage ob schon angekommen, Berechnung von v=Ah und alpha
r2=scalar(r,r);
if r2<0.000000000001 %1e-12, starke Forderung!
  disp(' ');
  disp(['Lösung nach ',num2str(j-1),' Schritten gefunden. Residuum = ',...
         num2str(r2)]);
  toc %Stoppuhr stoppt
  [X,Y]=meshgrid(0:hh:1);
  surf(X,Y,u); %Bild von der Lösung wird gezeichnet
  return;
end; %if 
v=prodd(h,bb);
r2h=scalar(h,v);
alpha=r2/r2h;
%=============================================================================
%3. Schritt, neue Lösung
u=u+alpha*h;
%=============================================================================
%4. Schritt, neues Residuum
r=r-alpha*v;
%=============================================================================
%5. Schritt, beta
r2neu=scalar(r,r);
beta=r2neu/r2; 
%=============================================================================
%6. Schritt, neue Richtung, dann bei Schritt 2 fortsetzen
h=r+beta*h;
%=============================================================================
if rem(j,50)==0  %alle 50 Schritte ein Zwischenergebnis
  disp(['Schritt ',num2str(j),' ist zu Ende, Residuum = ',num2str(r2neu)]);
end; %if rem
end; %for j %zum Ende der Schleife kommt es normalerweise nie
return        %cg-Programm ist hier zu Ende







































