driver_hw.m

nonlinear eq routines in matlab

%
%  driver.m
% 
%  use bisect.m and newt1d.m to compute an approximation
%  of the root of  f(x) = cos(x)*cosh(x)+1
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  Jan 17, 2002
%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
disp(' Find depth at which water main will freeze')
disp(' ')

f     = 'freeze';
fplot(f,[0,2]);
xlabel(' depth  x  in degree meter')
ylabel(' temperature at depth x and 60 days ')

exportfig(gcf,'freeze.eps','FontMode','Fixed','FontSize',14,...
            'LineMode','Fixed','LineWidth',1.5);



a     = 0;
b     = 2;

[a, b, x, ithist, iflag] = bisect( f, a, b );

fprintf('\n  Bisection Method \n ')
fprintf(' iter         a              b              c           f(c) \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e  %13.6e  \n ', ithist(i,:))
end

disp(['The Bisection method returned with iflag = ',int2str(iflag)])
disp(' Depth estimated by the Bisection method:')
disp([' x= ', num2str(x),' meters'])
      

disp(' ')
x = 0.01;
[y,yp]=feval(f, x);
y1    = feval(f, x+sqrt(eps));
ypfd  = (y1-y)/sqrt(eps);
disp([' derivative at x = ',num2str(x),'                  : f''(x) = ', num2str(yp)])
disp([' finite difference approx. of derivative :          ', num2str(ypfd)])
disp([' error = ', num2str(abs(yp-ypfd)),' should be of order ', num2str(sqrt(eps))])


[x, ithist, iflag] = newt1d( f, x );

fprintf('\n  Newon''s Method \n ')
fprintf(' iter        x           f(x)       -f(x)/f''(x)  \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e   \n ', ithist(i,:))
end
disp(['Newton''s method returned with iflag = ',int2str(iflag)])
disp(' Depth estimated by Newton''s method:')
disp([' x= ', num2str(x),' meters'])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' Hit return to continue')
pause
disp(' ')
disp(' ')
disp(' Find steady state temperature' )
disp(' ')



f     = 'steadytemp';
fplot(f,[200,400]);
xlabel(' T  in  K')
ylabel(' f(T) ')
exportfig(gcf,'steadytemp.eps','FontMode','Fixed','FontSize',14,...
            'LineMode','Fixed','LineWidth',1.5);



a     = 200;
b     = 400;


[a, b, T, ithist, iflag] = bisect( f, a, b );

fprintf('\n  Bisection Method \n ')
fprintf(' iter         a              b              c           f(c) \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e  %13.6e  \n ', ithist(i,:))
end

disp(['The Bisection method returned with iflag = ',int2str(iflag)])
disp(' Steady state temperature estimated by the Bisection method:')
disp([' T= ', num2str(T),'K = ', num2str(T-273.15),'C'])
      


T = 200;
[T, ithist, iflag] = newt1d( f, T );

fprintf('\n  Newon''s Method \n ')
fprintf(' iter        T           f(T)       -f(T)/f''(T)  \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e   \n ', ithist(i,:))
end
disp(['Newton''s method returned with iflag = ',int2str(iflag)])
disp(' Steady state temperature estimated by Newton''s method:')
disp([' T= ', num2str(T),'K = ', num2str(T-273.15),'C'])



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' Hit return to continue')
pause
disp(' ')
disp(' ')
disp(' Number of years to pay off loan' )
disp(' ')



f     = 'payoff_years';
fplot(f,[0,20]);
xlabel(' n in years')
ylabel(' f(n) ')
exportfig(gcf,'payoff_years.eps','FontMode','Fixed','FontSize',14,...
            'LineMode','Fixed','LineWidth',1.5);


a     = 0;
b     = 20;

[a, b, n, ithist, iflag] = bisect( f, a, b );

fprintf('\n  Bisection Method \n ')
fprintf(' iter         a              b              c           f(c) \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e  %13.6e  \n ', ithist(i,:))
end

disp(['The Bisection method returned with iflag = ',int2str(iflag)])
disp(' Number of years estimated by the Bisection method:')
disp([' n= ', num2str(n),' years '])
      

n = 0;
[n, ithist, iflag] = newt1d( f, n );

fprintf('\n  Newon''s Method \n ')
fprintf(' iter        n           f(n)       -f(n)/f''(n)  \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e   \n ', ithist(i,:))
end
disp(['Newton''s method returned with iflag = ',int2str(iflag)])
disp(' Number of years estimated by Newton''s method:')
disp([' n= ', num2str(n),' years '])


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' Hit return to continue')
pause
disp(' ')
disp(' ')
disp(' Interest rate to pay off loan' )
disp(' ')



f     = 'payoff_rate';
fplot(f,[0.0001,0.1]);
xlabel(' r ')
ylabel(' f(r) ')
exportfig(gcf,'payoff_rate.eps','FontMode','Fixed','FontSize',14,...
            'LineMode','Fixed','LineWidth',1.5);


a     = 0.0001;
b     = 0.1;

[a, b, r, ithist, iflag] = bisect( f, a, b );

fprintf('\n  Bisection Method \n ')
fprintf(' iter         a              b              c           f(c) \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e  %13.6e  \n ', ithist(i,:))
end

disp(['The Bisection method returned with iflag = ',int2str(iflag)])
disp(' Interest rate estimated by the Bisection method:')
disp([' r= ', num2str(r*100),' percent '])
      

r = 1;
[r, ithist, iflag] = newt1d( f, r );

fprintf('\n  Newon''s Method \n ')
fprintf(' iter        r           f(r)       -f(r)/f''(r)  \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e   \n ', ithist(i,:))
end
disp(['Newton''s method returned with iflag = ',int2str(iflag)])
disp(' Interest rate estimated by Newton''s method:')
disp([' r= ', num2str(r*100),' percent '])


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' Hit return to continue')
pause
disp(' ')
disp(' ')
disp(' Payment amount to pay off loan' )
disp(' ')



f     = 'payoff_amount';
fplot(f,[0,100000]);
xlabel(' p ')
ylabel(' f(p) ')
exportfig(gcf,'payoff_amount.eps','FontMode','Fixed','FontSize',14,...
            'LineMode','Fixed','LineWidth',1.5);


a     = 0;
b     = 100000;

[a, b, p, ithist, iflag] = bisect( f, a, b );

fprintf('\n  Bisection Method \n ')
fprintf(' iter         a              b              c           f(c) \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e  %13.6e  \n ', ithist(i,:))
end

disp(['The Bisection method returned with iflag = ',int2str(iflag)])
disp(' Payment amount estimated by the Bisection method:')
disp([' p= $', num2str(p)])
      

p = 0;
[p, ithist, iflag] = newt1d( f, p );

fprintf('\n  Newon''s Method \n ')
fprintf(' iter        p           f(p)       -f(p)/f''(p)  \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e   \n ', ithist(i,:))
end
disp(['Newton''s method returned with iflag = ',int2str(iflag)])
disp(' Payment amount estimated by Newton''s method:')
disp([' p= $', num2str(p)])




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' Hit return to continue')
pause
disp(' ')
disp(' ')
disp(' Solve Kepler''s equation' )
disp(' ')



f     = 'kepler';

M     = 1;
e     = 0.5;

it    = 0;
itmax = 100;
E     = 0;
Eold  = 1;

fprintf('\n  Fixed point iteration \n ')
fprintf(' iter         E     E - (M + e*sin(E))   \n ')
while (abs(E-Eold) >= 1.e-7 &  it < itmax )
  Eold = E;
  E    = M + e*sin(Eold);
  fprintf('%4d  %13.6e  %13.6e  \n ', it, Eold, E-Eold)
  it   = it+1;
end


disp(' Eccentric anomaly estimated by the fixed point iteration:')
disp([' E= ', num2str(E)])
      

E = 0;
[E, ithist, iflag] = newt1d( f, E );

fprintf('\n  Newon''s Method \n ')
fprintf(' iter        p           f(p)       -f(p)/f''(p)  \n ')
for i = 1:size(ithist,1)
   fprintf('%4d  %13.6e  %13.6e  %13.6e   \n ', ithist(i,:))
end
disp(['Newton''s method returned with iflag = ',int2str(iflag)])
disp(' Eccentric anomaly estimated by Newton''s method:')
disp([' E = ', num2str(E)])



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' Hit return to continue')
pause
disp(' ')
disp(' ')
disp(' Solve a problem with roots of higher multiplicity' )
disp(' ')


global p
f     = 'x2lnx';


p     = 6;
maxit = 200;
tolf  = 1.e-10;
tolx  = 1.e-10;
x0     = 0.8;


disp([' Solve  (x^2-1)^p ln(x) = 0 with p = ',int2str(p)] )

m = 1;
[x, ithist, iflag] = mnewt1d( f, x0, m, tolf, tolx, maxit );

fprintf('\n  Newon''s Method \n ')
fprintf([' iter        x           f(x)       -f(x)/f''(x)        m',...
         '      |x_{k+1}-1|/|x_k-1|\n '])
for i = 1:size(ithist,1)-1
   fprintf('%4d  %13.6e  %13.6e  %13.6e  %13.6e  %13.6e \n ', ...
           ithist(i,:), abs(ithist(i+1,2)-1)/abs(ithist(i,2)-1) )
end
i = size(ithist,1);
fprintf('%4d  %13.6e  %13.6e   \n ', ithist(i,1:3) )

disp(['Newton''s method returned with iflag = ',int2str(iflag)])
disp('Solution estimated by modified Newton method:')
disp([' x = ', num2str(x)])



m = p+1;
[x, ithist, iflag] = mnewt1d( f, x0, m, tolf, tolx, maxit );

fprintf('\n  modified Newon''s Method \n ')
fprintf([' iter        x           f(x)       -f(x)/f''(x)        m',...
         '      |x_{k+1}-1|/|x_k-1|   \n '])
for i = 1:size(ithist,1)-1
   fprintf('%4d  %13.6e  %13.6e  %13.6e  %13.6e  %13.6e \n ', ...
           ithist(i,:), abs(ithist(i+1,2)-1)/abs(ithist(i,2)-1) )
end
i = size(ithist,1);
fprintf('%4d  %13.6e  %13.6e   \n ', ithist(i,1:3) )

disp(['Modified Newton method returned with iflag = ',int2str(iflag)])
disp('Solution estimated by modified Newton method:')
disp([' x = ', num2str(x)])



m = 0;
[x, ithist, iflag] = mnewt1d( f, x0, m, tolf, tolx, maxit );

fprintf('\n  adaptive modified Newon''s Method \n ')
fprintf([' iter        x           f(x)       -f(x)/f''(x)        m',...
         '      |x_{k+1}-1|/|x_k-1|   \n '])
for i = 1:size(ithist,1)-1
   fprintf('%4d  %13.6e  %13.6e  %13.6e  %13.6e  %13.6e \n ', ...
           ithist(i,:), abs(ithist(i+1,2)-1)/abs(ithist(i,2)-1) )
end
i = size(ithist,1);
fprintf('%4d  %13.6e  %13.6e   \n ', ithist(i,1:3) )

disp(['Adaptive modified Newton method returned with iflag = ',int2str(iflag)])
disp('Solution estimated by modified Newton method:')
disp([' x = ', num2str(x)])