echo on
%%%%% ch05_exa.m %%%%%
% Run script and study the effects of each line.

echo off
format compact;
format short;
disp(sprintf('**** Examples from Chapter 05 ****'));


disp(sprintf('**** Example 5.7 Interval Bisection *'));
echo on
%%%% Bisection method for solving nonlinear equation f(x)=x^2-4sin(x)=0
a=1; b=3;		% Initial bracket [a b] is [1 3]
tol = 10^10*eps;	% Tolerance

%%%% Implement the first bisection step.
fa=a^2-4*sin(a)	% f(a)
fb=b^2-4*sin(b)	% f(b)
m=a+(b-a)/2
fm=m^2-4*sin(m)	% f(m)
if (sign(fa)== sign(fm))
  a=m
else b=m
end
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('i	a		f(a)		b		f(b)'))
i=0;
while ((b-a)>tol & i<30)
   i=i+1;
   disp(sprintf('%d	%f	%f	%f	%f', i, a,a^2-4*sin(a),b,b^2-4*sin(b)));
	m=a+(b-a)/2;
	fa=a^2-4*sin(a);	% f(a)
	fm=m^2-4*sin(m);	% f(m)
	if (sign(fa)== sign(fm))
      a=m;
   else b=m;
   end
end
disp(sprintf('%d	%f	%f	%f	%f', i+1, a,a^2-4*sin(a),b,b^2-4*sin(b)));
echo on
fzero('x^2-4*sin(x)',2)	% The root found by MatLab.


echo off
disp(sprintf('**** Example 5.x Fixed Point Iteration *'));
echo on
%%%% One of the fixed point formulation for solving nonlinear equation 
		% f(x)=x^2-4sin(x)=0 is g(x)=sqrt(4*sin(x))
x0=3; 			% Initial approximation of a root.
tol = 10^10*eps;	% Tolerance

%%%% Implement the first iteration step.
g0=sqrt(4*sin(x0))	% Solution from the first iteration
x1=g0	
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('iter	xi		f(xi)'))
i=0;
while (abs(x1-x0)>tol & i<30)
   i=i+1;
	disp(sprintf('%d	%f	%f', i,x0,x1));
   x0=x1;
   x1=sqrt(4*sin(x0));
end
disp(sprintf('%d	%f	%f', i+1,x0,x1));
echo on


echo off
disp(sprintf('**** Example 5.10 Newton''s Method *'));
echo on
%%%% Needs function f(x)=x^2-4sin(x)  and 
		% derivative f'(x)=2*x-4*cos(x) evaluation. 
x0=3; 			% Initial approximation of a root.
tol = 10^10*eps;	% Tolerance

%%%% Implement the first iteration step.
fx0=x0^2-4*sin(x0)	% Function value in x0
dx0=2*x0-4*cos(x0)	% Derivative value in x0
h=fx0/dx0
x1=x0 - h 		% Implement iteration step.	
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('iter	xi		f(xi)		f''(xi)		h'))
i=0;
while (abs(x1-x0)>tol & i<30)
   i=i+1;
	disp(sprintf('%d	%f	%f	%f	%f', i,x0,fx0, dx0, h));
   x0=x1;
	fx0=x0^2-4*sin(x0);		% Function value in x0
	dx0=2*x0-4*cos(x0);		% Derivative value in x0
	h=fx0/dx0;
	x1=x0 - h;		% Implement iteration step.	
end
disp(sprintf('%d	%f	%f	%f	%f', i+1,x0,fx0, dx0, h));
echo on


echo off
disp(sprintf('**** Example 5.12 Secant Method *'));
echo on
%%%% Needs only one evaluation of function f(x)=x^2-4sin(x)=0   
x0=1; x1=3; 		% Two initial guesses near a root.
tol = 10^10*eps;	% Tolerance

%%%% Implement the first iteration step.
fx0=x0^2-4*sin(x0)	% Function value in x0
fx1=x1^2-4*sin(x1)	% Function value in x1
h=x1-x0
x2=x1 - fx1*h/(fx1-fx0)	% Implement iteration step.	
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('iter	xi		f(xi)		h'))
i=0;
disp(sprintf('%d	%f	%+f	', i,x0,fx0,h));
while (abs(x1-x0)>tol & i<30)
   i=i+1;
 	disp(sprintf('%d	%f	%+f	%+f', i,x1,fx1,h));
   x0=x1;
   x1=x2; 
	fx0=fx1;		% Already known function value
	fx1=x1^2-4*sin(x1);		% Function value in x1
	h=x1-x0;
	x2=x1 - fx1*h/(fx1-fx0); 		% Implement iteration step.	
end
echo on

echo off
disp(sprintf('**** Example 5.12 Secant Method implemented by interpolation*'));
echo on
%%%% Interpolation function Fi(x)=ux+v=f(x) is a line, so
%fx0=u*x0+v and
%fx1=u*x1+v. This 2x2 system is solved for u and v and Fi(0) is calculated.
%This value becomes new x1 and old x1 becomes new x0. f(x1)=x1^2-4sin(x1)
%is calculated and new iteration step is started...
x0=1; x1=3; 		% Two initial guesses near a root.
tol = 10^10*eps;	% tolerance

%%%% Implement the first iteration step.
fx0=x0^2-4*sin(x0)	% Function value in x0
fx1=x1^2-4*sin(x1)	% Function value in x1
A=[x0 1; x1 1];
f=[fx0 fx1]';
coef=A\f;	%solve for coefficients
h=x1-x0;
x0=x1;
x1=-coef(2)/coef(1)
i=1;
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('iter	xi		f(xi)		h'))
disp(sprintf('%d	%f	%+f	', i,x1,fx1,h));
while (abs(x1-x0)>tol & i<30)
	fx0=fx1;		% Already known function value
	fx1=x1^2-4*sin(x1);		% Function value in x1
	A=[x0 1; x1 1];
	f=[fx0 fx1]';
	coef=A\f;	%solve for coefficients
	h=x1-x0;
	x0=x1;
	x1=-coef(2)/coef(1);
 	i=i+1;
 	disp(sprintf('%d	%f	%+f	%+f', i,x1,fx1,h));
end
echo on


echo off
disp(sprintf('**** Example 5.13 Inverse Quadratic Interpolation *'));
echo on
%%%% Needs one function evaluation/iteration for f(x)=x^2-4sin(x)=0, 
			% 3x3 system solution and some other arithmetic operations.
a=1; b=2; c=3; 		% Three initial guesses near a root.
tol = 10^10*eps;	% Tolerance

%%%% Implement the first iteration step.
fa=a^2-4*sin(a)		% Function value in a
fb=b^2-4*sin(b)		% Function value in b
fc=c^2-4*sin(c)		% Function value in c
u=fb/fc
v=fb/fa
w=fa/fc
p=v*(w*(u-w)*(c-b)-(1-u)*(b-a))
q=(w-1)*(u-1)*(v-1)
h=p/q
b1=b+p/q 		% Implement iteration step.
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('iter	xi		f(xi)		h'))
i=0;
disp(sprintf('%d	%f	%+f	', i,a,fa));
disp(sprintf('%d	%f	%+f	', i,b,fb));
disp(sprintf('%d	%f	%+f	', i,c,fc));
while (abs(b1-b)>tol & i<30)
   i=i+1;
   c=a;
	a=b;
	b=b1;
	fc=fa;		% Function value in c
	fa=fb;		% Function value in a
	fb=b^2-4*sin(b);		% Function value in b
	disp(sprintf('%d	%f	%+f	%+f', i,b,fb,h));
	u=fb/fc;
	v=fb/fa;
	w=fa/fc;
	p=v*(w*(u-w)*(c-b)-(1-u)*(b-a));
	q=(w-1)*(u-1)*(v-1);
	h=p/q;
	b1=b+p/q; 		% Implement iteration step.
end
echo on

echo off
disp(sprintf('**** Example 5.13 Inverse Quadratic Interpolation implemented by interpolation*'));
echo on
%%%% Use the inverse quadratic interpolation with the same method as
% in secant method for function f(x)=x^2-4sin(x)=0 
a=1; b=2; c=3; 		% Three initial guesses near a root.
tol = 10^10*eps;		% Tolerance

%%%% Implement the first iteration step by interpolant Fi(x)=uf(a)^2+vf(a)+c=a.
fa=a^2-4*sin(a)		% Function value in a
fb=b^2-4*sin(b)		% Function value in b
fc=c^2-4*sin(c)		% Function value in c
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('iter	xi		f(xi)		h'))
i=0;
disp(sprintf('%d	%f	%+f	', i,a,fa));
disp(sprintf('%d	%f	%+f	', i,b,fb));
disp(sprintf('%d	%f	%+f	', i,c,fc));
i=i+1;
A=[fa^2 fa 1;fb^2 fb 1;fc^2 fc 1;];	% function values in matrix because of inverse interpolation
points=[a b c]';	% approximation points, 
coef=A\points;
c=b;
b=a;
a=coef(3);	% the value of Fi(0) is new approximation point on x
h=a-b;
disp(sprintf('%d	%f	%+f	%+f', i,a,fa,h));
while (abs(a-b)>tol & i<30)
	fc=fb;		% Function value in c
	fb=fa;		% Function value in b
	fa=a^2-4*sin(a);		% Function value in a
	A=[fa^2 fa 1;fb^2 fb 1;fc^2 fc 1;];	% function values in matrix because of inverse interpolation
	points=[a b c]';	% approximation points, 
	coef=A\points;
	c=b;
	b=a;
   a=coef(3);	% the value of Fi(0) is new approximation point a on x
   h=a-b;
	i=i+1;
	disp(sprintf('%d	%f	%+f	%+f', i,a,fa,h));
end
echo on


echo off
disp(sprintf('**** Example 5.14 Linear Fractional Interpolation *'));
echo on
%%%% Needs one function evaluation/iteration for f(x)=x^2-4sin(x)=0 and 
			% some other arithmetic operations.
a=1; b=2; c=3; 		% Three initial guesses near a root.
tol = 10^10*eps;	% Tolerance

%%%% Implement the first iteration step.
fa=a^2-4*sin(a)		% Function value in a
fb=b^2-4*sin(b)		% Function value in b
fc=c^2-4*sin(c)		% Function value in c
p=(a-c)*(b-c)*(fa-fb)*fc	% Numerator
q=(a-c)*(fc-fb)*fa - (b-c)*(fc-fa)*fb	% Denominator
h=p/q
c1=c+h 		% Implement iteration step.
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('iter	xi		f(xi)		h'))
i=0;
disp(sprintf('%d	%f	%+f	', i,a,fa));
disp(sprintf('%d	%f	%+f	', i+1,b,fb));
disp(sprintf('%d	%f	%+f	', i+2,c,fc));
i=2;
while (abs(c1-c)>tol & i<30)
   i=i+1;
   a=b;
	b=c;
	c=c1;
	fa=a^2-4*sin(a);		% Function value in a
	fb=b^2-4*sin(b);		% Function value in b
	fc=c^2-4*sin(c);		% Function value in c
	disp(sprintf('%d	%f	%+f	%+f', i,c,fc,h));
	p=(a-c)*(b-c)*(fa-fb)*fc;	% Numerator
	q=(a-c)*(fc-fb)*fa - (b-c)*(fc-fa)*fb;	% Denominator
	h=p/q;
	c1=c+h; 		% Implement iteration step.
end
echo on


echo off
disp(sprintf('**** Example 5.14 Linear Fractional Interpolation implemented by interpolation*'));
echo on
%%%% Needs one function evaluation/iteration for f(x)=x^2-4sin(x)=0  
% 3x3 system solution and some other arithmetic operations.
a=1; b=2; c=3; 		% Three initial guesses near a root.
tol = 10^10*eps;	% Tolerance

%%%% Implement the first iteration step by interpolant Fi(x)=(x-u)/(wx-v)=0,
%it follows that x=u and this is the next point c, old c is new b and old b
%is new a.
fa=a^2-4*sin(a)		% Function value in a
fb=b^2-4*sin(b)		% Function value in b
fc=c^2-4*sin(c)		% Function value in c
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('iter	xi		f(xi)		h'))
i=0;
disp(sprintf('%d	%f	%+f	', i,a,fa));
disp(sprintf('%d	%f	%+f	', i,b,fb));
disp(sprintf('%d	%f	%+f	', i,c,fc));
i=i+1;
A=[1 a*fa -fa;1 b*fb -fb;1 c*fc -fc;];	% function values in matrix because of inverse interpolation
points=[a b c]';	% approximation points, 
coef=A\points;
a=b;
b=c;
c=coef(1);	% the value of Fi(0) is new approximation point c on x
h=c-b;
disp(sprintf('%d	%f	%+f	%+f', i,c,fc,h));
while (abs(c-b)>tol & i<30)
	fa=fb;		% Function value in a
	fb=fc;		% Function value in b
	fc=c^2-4*sin(c);		% Function value in c
A=[1 a*fa -fa;1 b*fb -fb;1 c*fc -fc;];	% function values in matrix because of inverse interpolation
points=[a b c]';	% approximation points, 
coef=A\points;
a=b;
b=c;
c=coef(1);	% the value of Fi(0) is new approximation point c on x
h=c-b;
	i=i+1;
	disp(sprintf('%d	%f	%+f	%+f', i,c,fc,h));
end
echo on



echo off
disp(sprintf('**** Example 5.15 Newton''s Method for n-dimensions *'));
echo on
%%%% Needs one evaluation of n functions F(x)=[x_1+2x_2-2 x_1^2+4x_2^2-4]'
			% and evaluation of Jacobian matrix J(x)=[1 2; 2x_1 8x_2]'
			% and solution of a system; all in each iteration.
tol = [10^10*eps 10^10*eps]';	% Tolerance
x0=[1 2]';  	% Initial guess       

%%%% Implement the first iteration step.
Jx0=[1 2; 2*x0(1) 8*x0(2)]	% Jacobian matrix for x0(df_i/dx_i)
Fx0=[x0(1)+2*x0(2)-2 x0(1)^2+4*x0(2)^2-4]'	% function value for x0 F(x)
s=Jx0\-Fx0		% Solve the system
x1=x0 + s 		% Update solution	
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('iter	x0		x1		f(x1)		f(x2)'))
i=0;
while (abs(x1-x0)>tol & i<30)
   i=i+1;
	disp(sprintf('%d	%f	%f	%f	%f', i,x0(1),x0(2), Fx0(1), Fx0(2)));
   x0=x1;
	Jx0=[1 2; 2*x0(1) 8*x0(2)];	% Jacobian matrix for x0(df_i/dx_i)
	Fx0=[x0(1)+2*x0(2)-2 x0(1)^2+4*x0(2)^2-4]';	% function value for x0 F(x)
	s=Jx0\-Fx0;		% Solve the system
	x1=x0 + s; 		% Update solution	
end
disp(sprintf('%d	%f	%f	%f	%f', i+1,x0(1),x0(2), Fx0(1), Fx0(2)));
echo on


echo off
disp(sprintf('**** Example 5.16 Broyden''s Method for n-dimensions *'));
echo on
%%%% The secant updating method, needs one evaluation of n functions F(x)=[x_1+2x_2-2 x_1^2+4x_2^2-4]'
				% and updating of approximate Jacobian J(x)=[1 2; 2x_1 8x_2]'
tol = [10^10*eps 10^10*eps]';	% Tolerance
x0=[1 2]';  	% Initial guess       

%%%% Implement the first iteration step.
Bx0=[1 2; 2*x0(1) 8*x0(2)]	% initial Jacobian approximation 
Fx0=[x0(1)+2*x0(2)-2 x0(1)^2+4*x0(2)^2-4]'	% function value for x0 F(x)
s=Bx0\-Fx0		% Solve the system
x1=x0 + s 		% Update solution	
echo off		% Implement iterations to the tolerance tol.
disp(sprintf('iter	x0		x1		f(x1)		f(x2)'))
i=0;
while (abs(x1-x0)>tol & i<30)
   i=i+1;
	disp(sprintf('%d	%f	%f	%f	%f', i,x0(1),x0(2), Fx0(1), Fx0(2)));
   Fx1=[x1(1)+2*x1(2)-2 x1(1)^2+4*x1(2)^2-4]';	% function value for x1 F(x)
	y0=Fx1-Fx0; 	
	Bx0=Bx0+(y0-Bx0*s)*s'/(s'*s);	% Update Jacobian
	s=Bx0\-Fx1;		% Solve the system
   x0=x1;
   Fx0=Fx1;
	x1=x0 + s; 		% Update solution	
end
disp(sprintf('%d	%f	%f	%f	%f', i+1,x0(1),x0(2), Fx0(1), Fx0(2)));
echo on