#Kafung Mok
#Homework 2
#Feb.2, 2014

#Note:HD1 uses the program NR1, and TestHouseholder uses B(which uses B1)
#thus, I've also included it in here. 

Help:=proc() print(`NR1(f,x,x0), HD1(f,d,x,x0), HD(f,d,x,x0,N), TestHouseholder(f,d,x,x0,eps,N)`):
print(`B1(f,x,a,b), B(f,x,a,b,eps)`):
end:

############
#Problem 1.#
############

#a:={1,2,3};
#b:={2,3,4};
#a union b; 

#f:=x->x^2 ;
#g:=x->x+1 ;
#(f@g)(x);

#Product(1-x^i, i=1..3);
#maximize(sin(x)+cos(x));

############
#Problem 2.#
############

#NR1(f,x,x0): inputs an expression f in x,
#a variable x, and an initial guess x0
#and peforms ONE iteration in the Newton-Raphson method.

NR1:=proc(f,x,x0)

#x0->x0-f(x0)/f'(x0))
x0- subs(x=x0,f)/subs(x=x0,diff(f,x)):

end:


#HD1(f,d,x,x0): inputs an expression f in x, a variable,
#a positive integer d, and an initial guess x0 and 
#performs ONE iteration using the Householder Method.
HD1:=proc(f,d,x,x0)

if d=1 then 
   NR1(f,x,x0):
else
   x0+normal((d*subs(x=x0,diff(1/f,x$(d-1))))/(subs(x=x0,diff(1/f,x$d)))):
fi:

end:


#HD(f,d,x,x0,N): performs N iterations using the Householder Method.
HD:=proc(f,d,x,x0,N) local x1,i:

x1:=x0:

for i from 1 to N do
 x1:=HD1(f,d,x,x1):
od:

x1:

end:


############
#Problem 3.# 
############

#B1(f,x,a,b): inputs an expression (like x^3-2)
#in a variable x, a variable x, and two (RATIONAL!) NUMBERS
#and outputs FAIL if subs(x=a,f) and subs(x=b,f) have
#the same sign a,a or b,b if of them is zero
#and either a,c or c,b (where c=(a+b)/2) is where
#f changes sign
B1:=proc(f,x,a,b) local c:

if subs(x=a,f)=0 then
  RETURN(a,a):
elif subs(x=b,f)=0 then
  RETURN(b,b):
elif sign(subs(x=a,f)*subs(x=b,f))=1 then
 RETURN(FAIL):

else

  c:=(a+b)/2:

    if sign(subs(x=a,f)*subs(x=c,f))=-1 then
       RETURN(a,c):
    else
       RETURN(c,b):
    fi:

fi:

end:


B:=proc(f,x,a,b,eps) local frank:

frank:=B1(f,x,a,b):

if frank=FAIL or frank[1]-frank[2]=0 then
 RETURN(frank):
fi:

while abs(frank[1]-frank[2])>eps do
 frank:=B1(f,x,frank):
od:

frank:

end:


#TestHouseholder(f,d,x,x0,eps,N): inputs an expression
#f in the variable x, a positive integer d, an intial guess x0,
#a small eps, and a pos. integer N, and outputs the successive errors
#abs(xn-alpha) for n from 1 to N where alpha is the
#appx. to the root of f(x)=0 found by the bisection method
TestHouseholder:=proc(f,d,x,x0,eps,N) local i,L,alpha,x1:

alpha:=B(f,x,x0-2,x0+2,eps):

if alpha=FAIL then
 RETURN(FAIL):
fi:

L:=[]:

x1:=x0:


for i from 1 to N do
 x1:=HD1(f,d,x,x1):
 L:=[ op(L), abs(x1-alpha[1])]:
od:
L:
end:

#Doing some examples in Maple shows that Householder does indeed converge
#at a rate of d+1.