% Harj2 Teht. 3 LV 1:n kuutiojuuret Newtonin menetelmällä
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab .m-tiedosto
% Voit kokeilla 
% use matlab
% matlab
% Siten vaikka vaikka vain cut/paste

p=[1 0 0 -1];   % Polynomi ilmaistaan kerroinvektorina , tässä z^3-1
%p=[1 0 0 0 -1]; % Yhtä hyvin z^4 -1 tms.
q=polyder(p);    % Derivaattapolynomin kertoimet
n=200;
x=linspace(-2,2,n); y=x;   % Tasavälinen jako n:ään osaan
[X,Y]=meshgrid(x,y);       % Pisteet (X(k),Y(j)) muod. xy-tason n x n - hilan
Z0=X+i*Y;Z=Z0;             % Z0 on kompleksihila ts. X(k)+i*Y(j)
                           % Homma tehtiin Matlabin tehokkailla matriisiope-
                           % raatioilla.
for k=1:10;                % Newton-iteraatio
   Z=Z-polyval(p,Z)./polyval(q,Z);  % polyval laskee p-kertoimien määräämän
end;                                % polynomin Z-pisteissä

%        Huomaa! X,Y,Z0,Z ovat kaikki n x n - matriiseja

juuret=roots(p);           % roots laskee polynomin (kerttiimet p) juuret.
            

j1=abs(juuret(1)-Z)<0.3;   % j1 on n x n-"bitti"matriisi, 1 siellä, missä ehto
j2=abs(juuret(2)-Z)<0.3;   %    on voimassa, samoin j2 ja j3
j3=abs(juuret(3)-Z)<0.3;


clf                      % Grafiikaruudun tyhjennys 
x1=j1.*X;y1=j1.*Y;       % Kerrotaan vastinalkioittain, x1 on n x n-matriisi,
                         % jossa X-arvo siellä, missä ehto voimassa, muuten 0
plot(x1(:),y1(:),'.b')   % Näin saadaan läpikäydyksi kaikki (x1(i),y1(j))-
                         % pisteet. Piirretään sinisellä pisteellä alkuperäisen
                         % hilan pisteet, jotka joutuivat 10 iteraatiolla 0.3:n
                         % päähän j1:stä

hold on                  % Pidetään edellinen grafiikka.

x2=j2.*X;y2=j2.*Y;       % Vastaavasti x2
plot(x2(:),y2(:),'.r')   % Piirretään nyt punaisella
x3=j3.*X;y3=j3.*Y;       % Vastaavasti x3
plot(x3(:),y3(:),'.g')   % Piirretään nyt vihreällä
                         % 
% Pisteet, joista lähdettäessä ei 10:llä iteraatiolla  päädytä 0.3:n päähän 
% mistään juuresta, jäävät valkoisiksi. (0.3 ja 10 ovat kokeiltavia arvoja.)

plot(real(juuret(1)),imag(juuret(1)),'*k')  % Piirretään vielä juuret 
                                            % (k, niinkuin black)
plot(real(juuret(2)),imag(juuret(2)),'*k')
plot(real(juuret(3)),imag(juuret(3)),'*k')
axis equal; axis square;                    % sama skaala
shg                                         % show graphics

%%%%%%%%% Tähän päättyy varsinainen perusajo. 
%         Lisää yksityiskohtia voi sitten hakea (vaikka lopun ikäänsä).


% Kuvaikkunasta voi valita suurennuslasin, hiiren vasemmalla zoomataan "sisään"
% ja oikealla "ulos" Mitä enemmän pikseleitä, sitä pitempi viive zoomaamisessa.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Kun on tarpeeksi zoomattu, alkavat pikselit loppua. Nyt voi halutulle
% alueelle lisätä pikseleitä. Jos ottaa samankokoisen hilan kuin äsken
% ja muuttaa vain rajoja, ei tule yllätyksiä laskenta-ajan suhteen.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=200;
x=linspace(-0.2,0.2,n); y=x;   % Tässä zoomataan keskelle kuvaa.
[X,Y]=meshgrid(x,y);
Z0=X+i*Y;Z=Z0;
for k=1:10;                    % Voi kannattaa nostaa hieman isommaksi.
   Z=Z-polyval(p,Z)./polyval(q,Z);
end;


j1=abs(juuret(1)-Z)<0.3;       % Voi kannattaa muutella näitä
j2=abs(juuret(2)-Z)<0.3;
j3=abs(juuret(3)-Z)<0.3;


hold on
x1=j1.*X;y1=j1.*Y;
plot(x1(:),y1(:),'.b')
hold on
x2=j2.*X;y2=j2.*Y;
plot(x2(:),y2(:),'.r')
x3=j3.*X;y3=j3.*Y;
plot(x3(:),y3(:),'.g') 
plot(real(juuret(1)),imag(juuret(1)),'*k')
plot(real(juuret(2)),imag(juuret(2)),'*k')
plot(real(juuret(3)),imag(juuret(3)),'*k')
shg

%%&
% Ja vielä lisää ...
%%%%%%%%%%%%%%%%%%%%%%
n=200;
x=linspace(-0.06,0.06,n); y=x;
[X,Y]=meshgrid(x,y);
Z0=X+i*Y;Z=Z0;
for k=1:10;
   Z=Z-polyval(p,Z)./polyval(q,Z);
end;


j1=abs(juuret(1)-Z)<0.3;
j2=abs(juuret(2)-Z)<0.3;
j3=abs(juuret(3)-Z)<0.3;


hold on
x1=j1.*X;y1=j1.*Y;
plot(x1(:),y1(:),'.b')
hold on
x2=j2.*X;y2=j2.*Y;
plot(x2(:),y2(:),'.r')
x3=j3.*X;y3=j3.*Y;
plot(x3(:),y3(:),'.g') 
plot(real(juuret(1)),imag(juuret(1)),'*k')
plot(real(juuret(2)),imag(juuret(2)),'*k')
plot(real(juuret(3)),imag(juuret(3)),'*k')
shg
