# Run the command "Execute Worksheet" in the Edit menu to see output.restart;Digits:=16; # Forces Maple to compute exactly 16 digits in every "evalf"interface(displayprecision=16); # Forces Maple to display all 16 of those digits2.1, #5: Use the Bisection method to find a solution accurate to within 10^(-5) for the following problems,# This program uses the bisection method on the given function f, starting from given endpoints
# left and right, and a given tolerance Tol. It computes the number of iterations required.
bisection:=proc(f,left,right,Tol)
local N,Count,RootFound,a,b,m,x:
N:=floor(ln((right-left)/Tol)/ln(2)):
Count:=1:
RootFound:=false:
a:=evalf(left):
b:=evalf(right):
print("Starting interval", [a,b]);
if(sign(f(a))<>sign(f(b))) then
while ( (RootFound=false) and (Count<=N) ) do
m := evalf(a+(b-a)/2):
if (evalf(f(m))=0) then
RootFound:=true:
x:=m:
else
if sign(evalf(f(m)))=sign(evalf(f(a))) then
a:=m:
else
b:=m:
end if:
x:=evalf(a+(b-a)/2):
end if:
Count:=Count+1:
print(Count, "Interval", [a,b], "Root", x);
end do:
else
print(`Function doesn't change sign.`):
end if:
end:prob2_1_5a:=x->x-2^(-x);prob2_1_5b:=x->exp(x)-x^2+3*x-2;prob2_1_5c:=x->2*x*cos(2*x)-(x+1)^2;prob2_1_5d:=x->x*cos(x)-2*x^2+3*x-1;bisection(prob2_1_5a,0,1,0.00001);bisection(prob2_1_5b,0,1,0.00001);bisection(prob2_1_5c,-3,-2,0.00001);bisection(prob2_1_5c,-1,0,0.00001);bisection(prob2_1_5d,0.2,0.3,0.00001);bisection(prob2_1_5d,1.2,1.3,0.00001);2.2, #3: The following four methods are proposed to compute 21^(1/3). Rank them in order, based on their apparent speed of convergence, assuming p0=1.# This program just iterates f, from a given initial condition, a total of Maxiter times.
iterator:=proc(f,initial,Maxiter)
local i,x:
x:=evalf(initial);
for i from 1 to Maxiter do
x:=evalf(f(x)):
print(x):
end do:
end:prob2_2_3a:=x->(20*x+21/x^2)/(21);prob2_2_3b:=x->x-(x^3-21)/(3*x^2);prob2_2_3c:=x->x-(x^4-21*x)/(x^2-21);prob2_2_3d:=x->(21/x)^(1/2);iterator(prob2_2_3a,1,10);iterator(prob2_2_3b,1,10);iterator(prob2_2_3c,1,4);iterator(prob2_2_3d,1,10);LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Ynevalf(21^(1/3)); # the true answerBy far the best method is from the function in part (b), which happens to be Newton's method. After 8 steps, it has already nailed the answer to 16 digits.The second best method is part (d), which gets the answer to three digits after 10 steps. The third best is part (a), which after 10 steps is only correct to one digit.The worst method is part (c), which immediately lands on the wrong fixed point and never leaves.2.2 #7. Use Theorem 2.2 to show that g(x) = Pi + 0.5*sin(x/2) has a unique fixed point on [0,2*Pi]. Use fixed-point iteration to find an approximation to the fixed point that is accurate to within 10^(-2). Use Corollary 2.4 to estimate the number of iterations required to achieve 10^(-2) accuracy, and compare this theoretical estimate to the number actually needed.g:=x->Pi+1/2*sin(x/2);maximize(D(g)(x),x=0..2*Pi);So the maximum of the derivative of g on the given interval is k = 1/4. Then Theorem 2.3 tells us the iterates of g converge to the unique fixed point.# Run the iterator program for 10 steps, with initial condition the midpoint
iterator(g,Pi,10);# Find the exact solution using prepackaged Maple program:
fsolve(g(x)=x,x=0..2*Pi);The iterator is correct to within 0.01 after two steps.Now Corollary 2.4 says the error is at most k^n/(1-k)*abs(p1-p0).fsolve(0.01=(1/4)^n/(1-1/4)*abs(Pi-g(Pi)),n);This tells us that we should require at least three steps, using the theoretical estimate. Of course, the theoretical estimate is always the worst possible upper bound, and we often get a little luckier than we'd expect.2.3 #5: Use Newton's method to find solutions accurate to within 10^(-4) for the following problems.newtonizer:=proc(f,initial,Tol,Maxiter)
local Count,AbsError,x,newx:
Count:=1:
x:=evalf(initial):
AbsError:=Tol+1:
print(Count,x):
while ( (AbsError>Tol) and (Count<=Maxiter) and (evalf(D(f)(x))<>0) ) do
newx:=evalf(x-f(x)/D(f)(x)):
Count:=Count+1:
AbsError:=evalf(abs(newx-x)):
x:=newx:
print(Count,x):
end do:
if (AbsError>Tol) then
print(`Couldn't reach tolerance.`):
end if:
if (evalf(D(f)(x))=0) then
print(`Derivative was zero.`):
end if:
end:prob2_3_5a:=x->x^3-2*x^2-5;prob2_3_5b:=x->x^3+3*x^2-1;prob2_3_5c:=x->x-cos(x);prob2_3_5d:=x->x-0.8-0.2*sin(x);Now to actually use Newton's method, we start with the midpoint of the interval and assume 10 steps are sufficient.newtonizer(prob2_3_5a,2.5,0.0001,10);newtonizer(prob2_3_5b,-2.5,0.0001,10);newtonizer(prob2_3_5c,Pi/4,0.0001,10);newtonizer(prob2_3_5d,Pi/4,0.0001,10);2.3 #7. Repeat Exercise 5 using the Secant Method.secantizer:=proc(f,left,right,Tol,Maxiter)
local Count,AbsError,firstx,secondx,newx,firsty,secondy:
Count:=1:
firstx:=evalf(left):
secondx:=evalf(right):
AbsError:=Tol+1:
print(0,firstx):
print(1,secondx):
firsty:=evalf(f(firstx)):
secondy:=evalf(f(secondx)):
while ( (AbsError>Tol) and (Count<=Maxiter) and (firsty<>secondy) ) do
newx:=evalf(secondx-secondy*(secondx-firstx)/(secondy-firsty)):
Count:=Count+1:
print(Count,newx):
AbsError:=evalf(abs(newx-secondx)):
firstx:=secondx:
firsty:=secondy:
secondx:=newx:
secondy:=evalf(f(secondx)):
end do:
if (AbsError>Tol) then
print(`Couldn't reach tolerance.`):
end if:
if (firsty=secondy) then
print(`Secant slope was zero.`):
end if:
end:To run the Secant Method, we start with the left and right endpoints of the interval and suppose it will be done in 20 iterations.secantizer(prob2_3_5a,1,4,0.0001,20);secantizer(prob2_3_5b,-3,-2,0.0001,20);secantizer(prob2_3_5c,0,Pi/2,0.0001,20);secantizer(prob2_3_5d,0,Pi/2,0.0001,20);2.3 #9. Repeat Exercise 5 using the method of False Position.# This program uses the bisection method on the given function f, starting from given endpoints
# left and right, and a given tolerance Tol. It computes the number of iterations required.
falsepositionizer:=proc(f,left,right,Tol,Maxiter)
local N,Count,AbsError,a,b,m,x,lefty,righty,midy:
N:=floor(ln((right-left)/Tol)/ln(2)):
Count:=1:
a:=evalf(left):
b:=evalf(right):
lefty:=evalf(f(a)):
righty:=evalf(f(b)):
if (sign(lefty)<>sign(righty)) then
print("Starting interval", [a,b]):
AbsError:=Tol+1:
while ( (AbsError>Tol) and (Count<=N) and (sign(righty)<>sign(lefty))) do
x := evalf(b-righty*(b-a)/(righty-lefty)):
midy:=evalf(f(x)):
if (midy=0) then
AbsError:=0:
elif sign(midy)=sign(lefty) then
AbsError:=abs(x-a):
a:=x:
lefty:=midy:
else
AbsError:=abs(x-b):
b:=x:
righty:=midy:
end if:
Count:=Count+1:
print(Count, "Interval", [a,b], "Root", x);
end do:
else
print(`Function doesn't change sign.`):
end if:
end:falsepositionizer(prob2_3_5a,1,4,0.0001,20);falsepositionizer(prob2_3_5b,-3,-2,0.0001,20);falsepositionizer(prob2_3_5c,0,Pi/2,0.0001,20);falsepositionizer(prob2_3_5d,0,Pi/2,0.0001,20);