####################################################################
##                          ONA 2001                              ##
##                                                                ##
##### Oto przyk³adowe rozwi±zanie zadania w³asnego Ay=lambda*y #####
##                                                                ##
##                       P Krzyzanowski                           ##
####################################################################

#
# Jest to skrypt Octave, aby wszystkie funkcje by³y w jednym miejscu;
# w praktyce, wszystkie poni¿ej zdefiniowane FUNKCJE umieszczaliby¶my
# w osobnych M-plikach.
#

#
# Definiujemy macierz A i obliczamy jej wymiar
#

A = diag([1,2,3,4,5]);
n = size(A,1);

#
# NOWO¦Æ: deklarujemy A i jej wymiar n jako zmienne GLOBALNE!
# Dziêki temu nie trzeba bêdzie podawaæ A jako parametru,
# co upro¶ci nam u¿ycie funkcji newton() bez modyfikacji!
# 

global A;
global n;

#
# Definiujemy funkcje my_f(x) oraz my_df(x).
# Zwykle robimy to w dwóch zewnêtrznych plikach, my_f.m i my_df.m,
# tym razem dla prostoty, zapiszemy je wprost w skrypcie;
# Zwróæmy uwagê na sposób dostêpu do zmiennych globalnych A i n!
#

function f=my_f(x)
global A;
global n;
y = x(1:n); lambda = x(n+1);
f = [ A*y - lambda*y ; 0.5*y'*y - 1.0 ];
endfunction

function df=my_df(x)
global A;
global n;
y = x(1:n); lambda = x(n+1);
df = [ A - lambda*eye(size(A)), -y;
       y'               ,  0 ];
endfunction


#
# Na wszelki wypadek, podajê tak¿e przyk³adow± implementacjê 
# funkcji newton() zmodyfikowanej tak, by dzia³a³a dla problemów 
# wielowymiarowych
# Proszê zwróciæ uwagê, ¿e ta funkcja nie ma ¿adnych innych parametrów,
# poza niezbêdnymi - w szczególno¶ci, nie odwo³uje siê wcale do macierzy A
# - z macierzy A korzystaj± jedynie funkcje definiuj±ce problem!
#

function [x, iter] = newton( f_name, df_name, x0, maxit, eps, flag)
x = x0;
iter = 1;
fx = feval(f_name,x); 
dfx = feval(df_name,x);
resid(iter) = norm(fx);
printf("[%d] f(x) = %e\n",iter-1,resid(iter));
while ((resid(iter) > eps) & (iter < maxit))
  x = x - dfx\fx;
  fx = feval(f_name,x);
  dfx = feval(df_name,x);
  iter = iter + 1;
  resid(iter) = norm(fx);
  printf("[%d] f(x) = %e\n",iter-1,resid(iter));
endwhile
if (flag > 0) 
  semilogy([0:iter-1],resid(1:iter),'-or;||f(x)||;');
endif
endfunction

# definiujemy wektor startowy x0 = [y, lambda]',
# gdzie      y = [1 1 0 0 0]
#       lambda = 10
#

x0 = [1 1 0 0 0 10]';

#
# rozwi±zujemy przy u¿yciu funkcji newton(), 
#

[x, iter] = newton( 'my_f', 'my_df', x0, 100, 1e-9, 0);

#
# oto moje wyniki:
#
# [0] f(x) = 1.204159e+01
# [1] f(x) = 7.046685e-01
# [2] f(x) = 3.898887e+01
# [3] f(x) = 9.641488e+00
# [4] f(x) = 2.279859e+00
# [5] f(x) = 4.297124e-01
# [6] f(x) = 3.690392e-02
# [7] f(x) = 4.031843e-04
# [8] f(x) = 5.175317e-08
# [9] f(x) = 9.272592e-16
