%Pascal/MATLAB-Kurs Oktober 2000, Gerhard Opfer
%Experiment mit g(x):=x^4+1. Berechung der komplexen Einzugsbereiche der vier
%verschiedenen Lösungen n1,n2,n3,n4 von g(x)=0 bei Anwendung des
%Newtonverfahrens. Die Lösungen sind 
%n1=0.5*sqrt(2)*( 1+i)= 0.7071 + 0.7071i, 
%n2=0.5*sqrt(2)*(-1+i)=-0.7071 + 0.7071i,
%n3=0.5*sqrt(2)*(-1-i)=-0.7071 - 0.7071i, 
%n4=0.5*sqrt(2)*( 1-i)= 0.7071 - 0.7071i.
n1=exp(i*pi/4);
n2=exp(3*i*pi/4);
n3=exp(5*i*pi/4);
n4=exp(7*i*pi/4);
erlaubt=100, %erlaubte Iterationszahl beim Newtonverfahren.
             %Kleine Zahlen (ca 10) führen dazu, daß weit von den Lösungen
             %entfernte Startpunkte innerhalb der erlaubten Iterationszahlen
             %nicht konvergieren. Bei größerer Zahl (ca 50) führen alle 
             %von den Achsen verschiedene Startpunkte zu einem Grenzwert,
             %allerdings ist nicht a-priore klar zu welchem.
klein=1e-4;  %Abbruchkriterium (könnte auch größer gewählt werden)
M=100;N=M;   %Diskretisierung
e=linspace(0,0,2*M+1); %gibt 0 0...0 (2*M+1)mal
id=[];
for m=1:2*M+1
  id=[id;e]; 
end; %for %Vorbesetzung einer ((2*M+1)x(2*M+1))-Matrix  mit der Zahl 0.  
size(id)
tic
for m=1:2*M+1
  for n=1:2*N+1
    zeler=0;
    x=(n-N-1-(m-M-1)*i)/M; %Das durch 2M+1 Punkte diskretisierte Quadrat
                           %[-1,1]x[-1,1] wird untersucht
    krit=10;
    while   (krit>klein)&(zeler<erlaubt)
      testvektor=[abs(x-n1),abs(x-n2),abs(x-n3),abs(x-n4)];
      krit=min(testvektor);
      zeler=zeler+1;
      if m~=M+1&n~=N+1 %x- oder y-Achse werden ausgelassen
        x=newton(x);else x=0;
      end; %if
    end; %while
    k=find(testvektor<=klein);
    if ~isempty(k)
      id(m,n)=k;
    end; %if
  end; %n
end; %m
toc %auf Ultra1 dauert 201 x 201 Matrix etwa 235 Sekunden bei klein = 1e-4
save newton201fein.erg id -ascii
