#hw2.txt, Feb. 1 2014, Anthony Zaleski Help:=proc(): print('Help(),Hd1(f,d,x,x0),Hd(f,d,x,x0,N),B1(f,x,a,b),B(f,x,a,b,eps)'); #self-ref is fun! print('TestHouseholder(f,d,x,x0,eps,N),TestConvergence(L)'); end: ##########Problem 1: Garvan pp. 30-60########## #Note: the output isn't shown. #To see it, you must execute this. #Sequences x:=1,2,3; y:=4,5,6; x,y; seq(i^2,i=1..5); #Sets Mathematicians:={Wiles,Newton,Frenkel,Gauss}; Actors:={Cruise,Kidman,Bale,Frenkel}; MathActors:=Mathematicians intersect Actors; member(Newton,MathActors); #Lists x:='x': y:='y': L1:=[x,y,z,y]; L2:=[a,b,c]; L:=[op(L1),op(L2)]; M:=[seq(1/n^2,n=1..5)]; MM:=M[2..3]; #Tables MathRankings:=table([(princeton)=excellent,(brown)=good,(rutgers)=potatoes]); MathRankings[brown]; MathRankings[rutgers]; #Arrays F:=array(1..2,1..2,[[3,4],[5,6]]); #Data Types L:=[seq(sin(n*Pi/4),n=1..8)]; type(L,list); convert(L,set); whattype({John, Paul, George, Tuna}); #Functions f:=x->x^2; g:=y->exp(y); simplify((f@g)(w)); An:=n->1/n^2; evalf(sum(An(n),n=1..10000)-Pi^2/6); #Limits f:=x->sin(x)/x; limit(f(x),x=0); #Differentiation del2:=(f,x,y)->diff(f,x,x)+diff(f,y,y); f:=(x,y)->log(x^2+y^2); del2(f(x,y),x,y); simplify(%); #Extrema f:=(x,y,z)->z/(x^2+y^2); assume(z>0); minimize(f(x,y,z),x,y); #Integration f:=x->exp(-x^2); int(f(x),x=-infinity..infinity); int(x*f(x),x); #Changing Variables with(student); I1:=Int(x^7/sqrt(x^16-1),x); I2:=changevar(u=x^8,I1,u); value(%); subs(u=x^8,%); #Taylor Series f:=x->1/(1-x); taylor(f(x),x=0,5); ##########Problem 2: Householder's Method########## ##Hd1 performs 1 iteration. Note d must be >1. Hd1:=proc(f,d,x,x0) local D1,D2,RHS: option remember: D1:=diff(1/f,x$(d-1)): D2:=diff(1/f,x$d): RHS:=x+d*normal(D1/D2): subs(x=x0,RHS): end: ##Hd iterates Hd1 N times: Hd:=proc(f,d,x,x0,N) local xn,i: xn:=x0: for i from 1 to N do xn:=Hd1(f,d,x,xn): end: xn: end: ##########Problem 3: Testing Householder########## #First, we'll require the bisection code from class: #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: #Next, TestHouseholder inputs an expression #f in the variable x, 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: #Finally, TestConvergence(L) takes in the list from #TestHouseholder and outputs the list of ratios of logs #of successive errors. These ratios estimate the #order of convergence. TestConvergence:=proc(L) local N: N:=nops(L): [seq(log(L[i+1])/log(L[i]),i=1..N-1)]: end: ##########Convergence Results########## #Using TestHouseholder and TestConvergence, we find #that in the special cases below, convergence is approx. #order d+1 as expected: (* #d=2: TestHouseholder(x^2-2,2,x,1,10^(-300),5): evalf(TestConvergence(%)); [3.485317838, 3.140265809, 3.044666866, 3.014670527] #d=3: TestHouseholder(x^2-2,3,x,1,10^(-300),5): evalf(TestConvergence(%)); [4.519537718, 4.114826075, 4.027905449, 1.541074534] *)