clear
reset(0)
echo on
%Solution to Problem 19
%
%The principle difficulty with finding the zeros
%of functions with repeated roots is that the
%higher derivatives vanish.  For the case of 
%double roots, this means that f'(x*) is zero.
%Recall the formula derived in the notes (I won't
%repeat that part of the derivation here):
%
% xs-x(i+1) = fpp(e)/2/fp(x(i))*(xs-x(i))^2
%
%where xs is the solution to the problem and e is
%in the interval [x(i),x(i+1)].  The difficulty is
%that fp(x(i)) vanishes as x(i) --> xs.  We can correct
%for this by expanding fp(x(i)) in a Taylor series
%about xs:
%
% fp(x(i)) = fp(xs) + (x(i)-xs)*fpp(eta)
%
%where eta is in the interval [x(i),xs].  Recognizing
%that fp(xs) = 0 for a double rooted function, and 
%substituting the second equation into the first, we get:
%
% xs-x(i+1) = -(fpp(e)/fpp(eta)/2)*(xs-x(i))
%
%           ~ -1/2 * (xs-x(i))
%
%provided the function does not have a triple root at xs.
%Thus Newton's method has a linear rate of convergence for
%double rooted functions with a constant C of 1/2.
%
%Now for the rest of the problem.  We use a modification
%of the example program used in class.  We are searching
%for the repeated roots of x*(1-x)^2 and x*(1-x)^3.  We
%use the function name fa.m for the first and fb.m for the
%second.  We also make use of the matlab command feval so
%that we can switch from one function to the other without
%having to rewrite things.
%
echo off
for flag=1:2

if flag == 1,
 disp('We first work with fa:')
 f='fa';
 fprime='faprime';
else
 disp('Now we work with fb:')
 f='fb';
 fprime='fbprime';
end

disp('Now for the bisection method:')
%We set the interval:
a=.5;
b=2;
y=[a:(b-a)/50:b];
figure(1)
plot(y,feval(f,y));
hold on
mid=(a+b)/2;
if feval(f,a)*feval(f,b)>0,
  disp('Bisection will not work with this interval')
else
  alpha(1)=a;
  beta(1)=b;
  mid(1)=(a+b)/2;
  plot(mid(1),feval(f,mid(1)),'o')
  for i=1:24,
    if feval(f,alpha)*feval(f,mid(i))<0,
      beta=mid(i);
    else,
      alpha=mid(i);
    end;
    mid(i+1)=(alpha+beta)/2;
    plot(mid(i),feval(f,mid(i)),'o')
	drawnow
  end
end
hold off

%Now let's see what the error looks like:
answer=1;
erbs=abs(mid-answer);
nbs=length(erbs);
if nbs>1,
coeff=polyfit(log(erbs(1:nbs-1)),log(erbs(2:nbs)),1);
figure(2)
plot(log(erbs(1:nbs-1)),log(erbs(2:nbs)),'o')
hold on
plot(log(erbs(1:nbs-1)),coeff(2)+coeff(1)*log(erbs(1:nbs-1)))
disp('The rate is')
rbs=coeff(1)
disp('The constant is')
cbs=exp(coeff(2))
hold off
else
rbs=NaN;
cbs=NaN;
end
pause

disp('Next we use Newton method.')
tol=1e-7;
xlast=a;
x(1)=(b+a)/2;
i=1;
figure(1)
plot(y,feval(f,y))
hold on
plot(x(i),feval(f,x(i)),'o')
pause(1)
while min(abs(feval(f,x(i))),abs(x(i)-xlast))>tol,
  x(i+1)=x(i)-feval(f,x(i))/feval(fprime,x(i));
  xlast=x(i);
  i=i+1;
  plot(x(i),feval(f,x(i)),'o')
  drawnow
end
hold off
%Let's look at the error as a function of iteration.
ernm=abs(x-answer);
nnm=length(ernm);
coeff=polyfit(log(ernm(4:nnm-1)),log(ernm(5:nnm)),1);
figure(2)
plot(log(ernm(1:nnm-1)),log(ernm(2:nnm)),'o')
hold on
plot(log(ernm(1:nnm-1)),coeff(2)+coeff(1)*log(ernm(1:nnm-1)))
disp('The rate is')
rnm=coeff(1)
disp('The constant is')
cnm=exp(coeff(2))
hold off
pause

disp('Now for the secant method:')
tol=1e-10;
xs(1)=(b+a)/2;
xs(2)=a+(b-a)/4;
i=2;
figure(1)
plot(y,feval(f,y))
hold on
plot(xs(1:2),feval(f,xs(1:2)),'o')
pause(1)
while min(abs(feval(f,xs(i))),abs(xs(i)-xs(i-1)))>tol,
  xs(i+1)=xs(i)-feval(f,xs(i))*...
          (xs(i)-xs(i-1))/(feval(f,xs(i))-feval(f,xs(i-1)));
  i=i+1;
  plot(xs(i),feval(f,xs(i)),'o')
  drawnow
end
hold off
%Let's look at the error as a function of iteration.
ersm=abs(xs-answer);
nsm=length(ersm);
coeff=polyfit(log(ersm(4:nsm-2)),log(ersm(5:nsm-1)),1);
figure(2)
plot(log(ersm(1:nsm-1)),log(ersm(2:nsm)),'o')
hold on
plot(log(ersm(1:nsm-1)),coeff(2)+coeff(1)*log(ersm(1:nsm-1)))
disp('The rate is')
rsm=coeff(1)
disp('The constant is')
csm=exp(coeff(2))
hold off
pause
disp('In summary, we have the following rates and constants');
disp('    bisection  newton  secant');
table=[rbs,rnm,rsm;cbs,cnm,csm];
disp(table);
end
echo off



___________________________________

cut out these programs and store them separately.

function y=fa(x)
y=x.*(x-1).^2;

function y=faprime(x)
y=(x-1).^2+2*x.*(x-1);

function y=fb(x)
y=x.*(x-1).^3;

function y=fbprime(x)
y=(x-1).^3+3*x.*(x-1).^2;