%% Numerik D-MAVT, RW (Nipp), SS06
%
% This is the M-file solution to the fifth exercise sheet.
% To run, save this file to your account and ensure that the "Current directory"
% in MATLAB is the same as that where the file is saved.
%
%
% Then to run the file type "s5" in the command window

function s5
%Aufgabe 3
disp(['Aufgabe 3 Newton Method']);
rtol=1e-10;
% Newton
x=[0.5;0.4]; i=0;
while (1)
  f=f3(x(1),x(2)); J=df3(x(1),x(2)); 
  xold=x;  x=x-J\f; 
  if (norm(x-xold)<rtol*norm(x)) break; end
  i=i+1;
end
disp(['Newton converged after ',int2str(i),' iterations:']);
disp(['  x=',num2str(x(1)),',y=',num2str(x(2)),...
      ' with ||f||_2=',num2str(norm(f))]);
% output:
% Aufgabe 3 Newton Method
% Newton converged after 4 iterations:
%   x=0.31076,y=0.2295 with ||f||_2=6.2287e-13

% Quasi Newton
disp(['Aufgabe 3 Quasi Newton Method'])
x=[0.5;0.4]; J0=df3(x(1),x(2)); [L,R]=lu(J0); i=0;
while (1)
  f=f3(x(1),x(2));
  xold=x; x=x-R\(L\f);
  if (norm(x-xold)<rtol*norm(x)) break; end
  i=i+1;
end
disp(['Quasi-Newton converged after ',int2str(i),' iterations:']);
disp(['  x=',num2str(x(1)),',y=',num2str(x(2)),...
      ' with ||f||_2=',num2str(norm(f))]);
% output:
% Aufgabe 3 Quasi Newton Method
% Quasi-Newton converged after 17 iterations:
%    x=0.31076,y=0.2295 with ||f||_2=6.2424e-11

%Aufgabe 4
disp(['Augabe 4 f(x) a0=0.1 b0=2'])
format compact; 
TOL0=1e-1; TOL=1e-10;

% *** bisection and secants for f(x)
a0=.1; b0=2; 
[a,b]=bisec(a0,b0,TOL0,'f');
x=secants(a,b,TOL,'f');
disp(['solution with bisec and secants: x = ',num2str(x,6),...
      ' (f(x) = ',num2str(feval('f',x),6),')']);
% output:
% Augabe 4 f(x) a0=0.1 b0=2
% bisec:   I = [0.1,2]
% bisec:   I = [1.05,2]
% bisec:   I = [1.525,2]
% bisec:   I = [1.525,1.7625]
% bisec:   I = [1.525,1.6437]
% bisec converged:   I = [1.525,1.5844]
% secants: x = 1.525
% secants: x = 1.5573
% secants: x = 1.5571
% secants: x = 1.5571
% secants converged: x = 1.5571
% solution with bisec and secants: x = 1.55715 (f(x) = -2.22045e-16)

% *** bisection and secants for g(x): different initial intervals
disp(['Augabe 4 g(x) a0=-12 b0=-1.7'])
a0=-12; b0=-1.7;
[a,b]=bisec(a0,b0,TOL0,'g');
x=secants(a,b,TOL,'g');
disp(['solution with bisec and secants: x = ',num2str(x,6),...
      ' (g(x) = ',num2str(feval('g',x),6),')']);
% output:
% Augabe 4 g(x) a0=-12 b0=-1.7
% bisec:   I = [-12,-1.7]
% bisec:   I = [-6.85,-1.7]
% bisec:   I = [-4.275,-1.7]
% bisec:   I = [-2.9875,-1.7]
% bisec:   I = [-2.3438,-1.7]
% bisec:   I = [-2.0219,-1.7]
% bisec:   I = [-2.0219,-1.8609]
% bisec converged:   I = [-2.0219,-1.9414]
% secants: x = -2.0219
% secants: x = -1.9982
% secants: x = -1.9999
% secants: x = -2
% secants: x = -2
% secants converged: x = -2
% solution with bisec and secants: x = -2 (g(x) = 0)

disp(['Augabe 4 g(x) a0=-1.5 b0=0'])
a0=-1.5; b0=0;
[a,b]=bisec(a0,b0,TOL0,'g');
x=secants(a,b,TOL,'g');
disp(['solution with bisec and secants: x = ',num2str(x,6),...
      ' (g(x) = ',num2str(feval('g',x),6),')']);
% output:
% Augabe 4 g(x) a0=-1.5 b0=0
% bisec:   I = [-1.5,0]
% bisec:   I = [-1.5,-0.75]
% bisec:   I = [-1.125,-0.75]
% bisec:   I = [-1.125,-0.9375]
% bisec converged:   I = [-1.0312,-0.9375]
% secants: x = -1.0312
% secants: x = -1.0009
% secants: x = -0.99998
% secants: x = -1
% secants: x = -1
% secants converged: x = -1
% solution with bisec and secants: x = -1 (g(x) = -0) 

disp(['Augabe 4 g(x) a0=0 b0=5'])
a0=0; b0=5;
[a,b]=bisec(a0,b0,TOL0,'g');
x=secants(a,b,TOL,'g');
disp(['solution with bisec and secants: x = ',num2str(x,6),...
      ' (g(x) = ',num2str(feval('g',x),6),')']);
% output:
% Augabe 4 f(x) a0=0 b0=5
% bisec:   I = [0,5]
% bisec:   I = [0,2.5]
% bisec:   I = [0,1.25]
% bisec:   I = [0.625,1.25]
% bisec:   I = [0.9375,1.25]
% bisec:   I = [0.9375,1.0938]
% bisec converged:   I = [0.9375,1.0156]
% secants: x = 0.9375
% secants: x = 0.99916
% secants: x = 1
% secants: x = 1
% secants: x = 1
% secants converged: x = 1
% solution with bisec and secants: x = 1 (g(x) = 0)

% *** last part: special parameters for g(x)
disp(['Augabe 4 g(x) a0=-0.9 b0=1.1'])
a0=-.9; b0=1.1; TOL0=3; TOL=1e-10;
[a,b]=bisec(a0,b0,TOL0,'g');
x=secants(a,b,TOL,'g');
% output:
% Augabe 4 f(x) a0=-0.9 b0=1.1
% bisec converged:   I = [-0.9,1.1]
% secants: x = -0.9
% secants: x = -0.41395
% ??? Error using ==> secants
% x = -0.99191 out of interval




%-------------------------------------------------------------------

function f=f3(x,y)
f=[exp(x*y)+x^2+y-1.4;x^2+y^2+x-0.46];

% ------------------------------------------------------------------

function J=df3(x,y)
J=[ y*exp(x*y)+2*x  x*exp(x*y)+1;
    2*x+1           2*y];

% ------------------------------------------------------------------


function ff=f(x)
ff=x-2+log(x);

% ------------------------------------------------------------------

function gg=g(x)
gg=(x-1).*(x+1).*(x+2);

% ------------------------------------------------------------------

function [a,b] = bisec(a0,b0,tol,F)
% bisection to find zeros of F(x)
% in:   a0      initial interval is [a0,b0]
%       b0
%       tol     stop if |b-a| < tol
%       F       function to find zeros of (string)
% out:  a       final interval with |b-a| < tol and a<= x <=b
% auth: fab@math.ethz.ch
if (feval(F,a0)*feval(F,b0)>0) 
  error('function has same sign at a and b'); 
end;
a=a0; b=b0;
while (abs(b-a)>=tol)
  disp(['bisec:   I = [',num2str(a),',',num2str(b),']']);
  m = 0.5*(a+b);
  if (feval(F,a)*feval(F,m)<0) b=m;
  else a=m;
  end
end 
disp(['bisec converged:   I = [',num2str(a),',',num2str(b),']']);

% ------------------------------------------------------------------

function x = secants(x0,x1,rtol,F)
% method of secants to find zeros of F in interval [x0,x1]
% in:   x0,x1   starting values for iteration
%       rtol    relative tolerance: 
%               stop if abs(x1-x0)<=rtol*abs(x) 
%       F       function to find zeros of (string)
% out:  x       approximation for zero of F
% auth: fab@math.ethz.ch
a=x0; b=x1; x=x0;
while (abs(x-x1)>rtol*abs(x))
  disp(['secants: x = ',num2str(x)]);
  x0=x1; x1=x;
  x = (x0*feval(F,x1)-x1*feval(F,x0))/(feval(F,x1)-feval(F,x0));
  if (x>b)|(x<a) 
    error(['x = ',num2str(x),' out of interval']); 
  end;
end;
disp(['secants converged: x = ',num2str(x)]);













!DSPAM:4462ff6f195125260519334!
