


1==1;



function [x,it,yf,er]=newton(f,df,x0,tol=1e-6,maxit=100,r=0)
#implementacja m Newtona do testow (zakladamy ze znamy rozwiazanie r tzn f(r)=0)
#Input
#f,df - wskazniki do funkcji obliczajacych wartosc funckji (f) i jeje pochodnej (df)
#x0 - przyblizenie poczatkowe
#r - znane rozwiazanie tzn f(r)=0
#tol - tolerancka - metoda zatrzymuje sie jak |f(x)|<=tol 
#maxit - max ilosc iteracji
#
#output
#x przyblizenie zera f
#it - ilosc iteracji
#yf - wartosc f(x)
#er - wektor z bledami er(k)=|x_k-r| jesli r nie podane to wektor z iteracjami (r=0) 

 er=zeros(maxit,1); 
 er(1)=x0;
 x=x0;
 it=0;
 if(nargin>5)
  e=x-r;
 endif
 yf=f(x);

 while((it<maxit)&&(abs(yf)>tol))
  x=x-yf/df(x);
  it++;
  yf=f(x);
  er(it+1)=x;
 endwhile
 er=er(1:it);
 if(nargin>5)
  er=er-r;
  erp=[1;er(1:it-1)];
  err=[0;er(2:it)];
  printf("nr iteracji,   blad,    blad/bl_popr, blad/(bl_popr)^2, blad/(bl_popr)^3\n");
  disp([(1:it)', er, err./erp, err./(erp.*erp), err./(erp.*erp.*erp)]);
  printf("W przypadku zbieznosci kwad - 3cia kolumna zbiega do stalej>0\n");
 endif
endfunction

f=@(x) x*x -2;
df=@(x) 2*x;
[y,it,yf]=newton(f,df,2)
1/(2*y) #do tego powinna zbiec 2ga kolumna

[y,it,yf]=newton(f,df,2 ,1e-6,100,r=sqrt(2))


#kolejne przyklady
disp("\ncos(x)=0\n");
[y,it,yf]=newton(@cos,@(x) -sin(x),2.3 ,1e-6,100,r=pi/2)
disp("(x^3-...=0 -zero 3-krotne rozw x=1");
[y,it,yf]=newton(@(x) x*x*x-3*x*x+3*x-1,@(x) 3*x*x-6*x+3, 1.8 ,1e-6,100,r=1.0)
#zbieznosc liniowa z wykladnikiem 2/3 tzn 3 kolumna powinna zbiec do 2/3
disp("\nx^8=256=0\n");
[y,it,yf]=newton(@(x) x^8-256,@(x) 8*x^7, 2.5 ,1e-6,100,r=2.0)
disp("\nexp(x)=2\n");
[y,it,yf]=newton(@(x) exp(x)-2.0,@exp, 1 ,1e-8,100,r=log(2))

#m Herona dla mantysy
m=1+rand; #m in [1,2)
[y,it,yf]=newton(@(x) x*x-m,@(x) 2*x,2 ,1e-6,100,sqrt(m))

function [x,it,yf]=newtonap(f,x0,del=1e-8,tol=1e-6,maxit=100,r=0)
#implementacja m Newtona do testow (zakladamy ze znamy rozwiazanie r tzn f(r)=0)
#pochodna przyblizamy ilorazem roznicowym z parametrem del tzn 
#df(x) \approx (f(x(1+del))-f(x))/del
#Input
#f - wskaznik do funkcji obliczajacej wartosc funckji (f) 
#x0 - przyblizenie poczatkowe
#r - znane rzowiazanie tzn f(r)=0
#
#tol - tolerancka - metoda zatrzymuje sie jak |f(x)|<=tol 
#maxit - max ilosc iteracji
#przyklad wywolania
#
#f=@(x) x*x -2;df=@(x) 2*x;[y,it,yf]=newtonap(f,2)
#1/(2*y) #do tego powinna zbiec 2ga kolumna
#f=@(x) x*x -2;[y,it,yf]=newtonap(f,2,1e-8,1e-6,100,r=sqrt(2))


x=x0;it=0;
if(nargin>5)
 e=x-r;
endif
yf=f(x);

while((it<maxit)&&(abs(yf)>tol))
Delt=x*del;
iDf=Delt/(f(x+Delt)-f(x)); #odwrotnosc il roznicowego przyblizajacego pochodna
x=x-yf*iDf;
it++;
yf=f(x);
if(nargin>5)
 ep=e;e=x-r;
 [it e e/ep e/(ep*ep) e/(ep*ep*ep)]
endif

endwhile

endfunction

[y,it,yf]=newtonap(@(x) x*x-2,x0=2,del=1e-8,tol=1e-6,maxit=100,r=sqrt(2))

#zmniejszamy parametr przy ilorazie przyblizajacym pochodna
[y,it,yf]=newtonap(@(x) x*x-2,x0=2,del=1e-10,tol=1e-6,maxit=100,r=sqrt(2))
#ciut lepiej - prosze porownac z wynikiem dokladnej met Newtona

[y,it,yf]=newtonap(@(x) x*x-2,x0=2,del=1e-12,tol=1e-6,maxit=100,r=sqrt(2))
#nie lepiej
[y,it,yf]=newtonap(@(x) x*x-2,x0=2,del=1e-14,tol=1e-6,maxit=100,r=sqrt(2))
#gorzej..

#zwiekszmy
[y,it,yf]=newtonap(@(x) x*x-2,x0=2,del=1e-4,tol=1e-6,maxit=100,r=sqrt(2))
#gorzej
[y,it,yf]=newtonap(@(x) x*x-2,x0=2,del=1e-2,tol=1e-6,maxit=100,r=sqrt(2))


#metoda siecznych


function [x,it,yf]=sieczne(f,x0,x1=x0-1e-8*f(x0)/(f(x0+1e-8)-f(x0)),del=1e-8,tol=1e-6,maxit=100,r=0)
#implementacja m siecznych do testow (zakladamy ze znamy rozwiazanie r tzn f(r)=0)
#pochodna przyblizamy ilorazem roznicowym  na 2 ostatnic hiteracjach  
#df(x) \approx (f(x(1+del))-f(x))/del
#Input
#f - wskaznik do funkcji obliczajacej wartosc funckji (f) 
#x0 - przyblizenie poczatkowe
#r - znane rzowiazanie tzn f(r)=0
#
#tol - tolerancka - metoda zatrzymuje sie jak |f(x)|<=tol 
#maxit - max ilosc iteracji
#przyklad wywolania
#
#f=@(x) x*x -2;df=@(x) 2*x;[y,it,yf]=sieczne(f,2)
#1/(2*y) #do tego powinna zbiec 2ga kolumna
#f=@(x) x*x -2;[y,it,yf]=sieczne(f,2,2.00001,1e-8,1e-6,100,r=sqrt(2))


x=x1;xp=x0;
it=0;
if(nargin>5)
 ep=xp-r;e=x-r;
endif
yfp=f(xp);
yf=f(x);
while((it<maxit)&&(abs(yf)>tol))
iDf=(x-xp)/(yf-yfp);
xp=x;x=x-yf*iDf;
it++;
yfp=yf;yf=f(x);
if(nargin>5)
 epp=ep;ep=e;e=x-r;
 [it e e/ep e/(ep*epp) e/(ep*ep) e/(ep*ep*ep)]
endif

endwhile

endfunction


[x,it,yf]=sieczne(@(x) x*x-2,2,x1=x0-1e-8*f(x0)/(f(x0+1e-8)-f(x0)),del=1e-8,tol=1e-6,maxit=100,r=sqrt(2))


#porownanie
xx0=[2,4,10,0.1,001];
M=length(xx0);
f=@(x) x*x-2;df=@(x) 2*x;
noit=zeros(M,4);
for k=1:M,
 x0=xx0(k); 
 [xs,its]=sieczne(f,x0);
 [xn,itn]=newton(f,df,x0);
 [xan,itan]=newtonap(f,x0);
 noit(k,:)=[x0,itn,itan,its];
endfor
[" x0" " newton" " newton z ap poch" " sieczne"]
noit
