%%
%% Listing 7.5 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function [x, resid, info, output] = newton(nazwa_f, x0, atol, rtol, step, maxit)
%
%	[x, resid, info, output] = newton(nazwa_f, x0, atol, rtol, step, maxit)
%
% Rozwiązuje układ równań nieliniowych F(x) = 0, gdzie F jest zadana M-funkcją o 
% nazwie "nazwa_f".
% Parametry wejścia:
% 	nazwa_f - nazwa M-funkcji zadającej funkcję F(). 
% 		Ta funkcja *musi* przyjmować także argumenty macierzowe!
% 	x0 - przybliżenie początkowe rozwiązania
%	atol, rtol - tolerancja błędu residualnego
%	step - co ile iteracji uaktualniać macierz pochodnej; 
%		step = 1 daje standardową metodę Newtona, 
%		step > maxit daje metodę cięciw
%	maxit - maksymalna dopuszczalna liczba iteracji
% Wyjście:
%	x - końcowe przybliżenie rozwiązania
%	resid - wektor zawierający historię zbieżnosci, 
%		resid(i) - wartość normy residuum w i-tej iteracji
%	info - kod zakończenia: = 1 - sukces, = 0 - przekroczona maxit
%	output - struktura zawierająca statystyki

% kilka pomocniczych definicji
epsq = sqrt(eps);	% pierwiastek z podwojonej precyzji maszyny
unit = ones(size(x0));	% wektor złożony z samych jedynek

% początkowe rozwiązanie przybliżone
x = x0;
	
% wyznaczenie residuum początkowego
fx = feval(nazwa_f, x);
iter = 1;
resid(iter) = norm(fx);

while( (resid(iter) >= (atol + rtol*resid(1))) & (iter < maxit) )
	if (mod(iter-1, step) == 0)
		% aproksymacja pochodnej - wszystkie kolumny jednocześnie
		h = epsq*max(abs(x),unit);
		xh = x*unit' + diag(h);
		Df = (feval(nazwa_f, xh) - fx*unit')./(unit*h');
		% rozkład macierzy pochodnej do wykorzystania
		[L, U, p] = lu(Df, 'vector');
	end
	% poprawka przybliżonego rozwiązania, x = x - Df\fx;
	x = x - U\(L\fx(p));
	
	% wyznaczenie residuum
	fx = feval(nazwa_f, x);
	iter = iter + 1;
	resid(iter) = norm(fx);
end
info = ~(iter == maxit);
output.iterations = iter - 1;
end
