Author of Code: Charles Cadena
Overall, the format of your website is very easy to follow.
The size of the screenshots of the problem is too large compared to the rest of the text. I would suggest decreasing its size.
For the image under Part C of the Numerical Solution, is there something that was cut off?
In the text in the numerical solution part B, there is a typo of the word "chance."
I would recommend elaborating on how you choice the bounds for the secant method and why.
I would suggest relating the roots and solutions back to the problem and how it relates to the problem. For example, for a root you could say that it represents the tension of the cable caused by forces A, B, and C.
In the numerical solution part A, there is an unnecessary period present.
I would suggest writing where you got the sample calculation equation from in part B of the analytical solution. Without it, the reader may not understand its origin.
I would suggest expanding on the reasoning in part D. Although Brent's Method and the fzero function both found the exact same root, it would be helpful to know why.
I would recommend re-reading the error and analysis section for typos. It might be a good idea to copy and paste it into Word to catch the mistakes.
The code is fully functioning, so there isn't much that needs to be changed. Most of the feedback provided is nitpicking.
For the plots, I recommend slightly increasing the thickness of the lines and markers using the 'LineWidth' command.
I would also suggest increasing the size of the markers using the 'MarkerSize' syntax.
The true percent error plot shows the difference in error between the two methods for each iteration perfectly.
The format of the code is very spaced out with defined sections which allows the reader to stay focused easier.
The code architecture is consistent throughout with only a few arbitrary spaces.
The syntax all seems proper.
Lines 124-125, 196, 205-106: Comment how you determined the bounds.
Line 302: I would suggested being more descriptive with the axis title. For example, the label could be "y-component of tension at A."
Lines 17, 300: I would include a black, horizontal line representing the x-axis because just glancing at the plot, the x-axis and root is not obvious.
Line 16: Instead of writing out 'Color,' I recommend placing the 'r' before the '-'. This should change the line color to red without the need for the 'Color' command.
Lines 10, 198, 200: To stay consistent with the spacing style, the stray spaces could be removed.
In conclusion, the website needs a bit more work than the code. The code works perfectly and there are only a few minor aspects that should be changed. In the website, I would recommend elaborating more and relating the solutions back to the context of the problem. There are multiple typos throughout the text, so I suggest using a proof-reading strategy to fix them. Overall, great job!
%% Part A and B (Brent's Method)
% Brent's Method will be used to solve Parts A, B, and C of the problem.
% The Brent Method will then be compared to the secant method.
% The pseudocode on pg. 165 of "Numerical Methods" was used to code the
% Brent Method.
% Set Constants and Equations
eps=2.2204460*10^(-16); %error theshold
f=@(x) -12 - 21*x+x.^2-2.4*x.^3; %governing equation
x0=[-1 1]; % x bounds for initial guesses
z=fzero(f, x0); %analytical root
i=-1:0.01:1;
figure(1);
hold all;
plot(i, f(i),'-', 'Color', 'r');
plot(z, f(z), 'bo');
title('Analytical Solution', 'Fontsize', 18);
xlabel('x','Fontsize', 15);
ylabel('y','Fontsize', 15);
ylim([-1 1]);
grid on;
hold off;
%establish upper and lower bounds
a=-10;
b=10;
%evaluate func at a and b
fa=f(a);
fb=f(b);
c=a; %new point equal to a
fc=fa; %fc same value as fa
d=b-a;
e=d; %another point equal to e
brentmethoderror=zeros(1,100);
count=0;
while (1)
count = count+1;
if fb==0
break;
end
if (sign(fa) == sign(fb))
a=c;
fa=fc;
d=b-c;
e=d;
end
if(abs(fa)<abs(fb))% x values then switched to the new points made (a and b)
c=b;
b=a;
a=c;
fc=fb;
fb=fa;
fa=fc;
end
m=0.5*(a-b);
tolerance=2*(eps)*max(abs(b),1);
if (abs(m)<tolerance || fb==0)
break;
end
if (abs(e)>=tolerance && abs(fc)>abs(fb)) %open method used for for speed. if unacceptable answer is generated, secant method is chosen
s=fb/fc;
if (a==c)
p=2*m*s;
q=1-s;
else %if secant method generates unacceptable asnwer, inverse quadratic interpolation (IQI) is then used
q=fc/fa;
r=fb/fa;
p=s*(2*m*q*(q-r)-(b-c)*(r-1));
q=(q-1)*(r-1)*(s-1);
end
if(p>0)
q=-q;
else
p=-p;
end
if ((2*p<3*m*q-abs(tolerance*q)) && (p < abs(0.5*e*q)))
e=d;
d=p/q;
else
d=m;
e=m;
end
else %if IQI generates unacceptable answer, bisection method is used
d=m;
e=m;
end
c=b;
fc=fb;
if(abs(d)>tolerance)
b=b+d;
else
b=b-sign(b-a)*tolerance;
end
fb=f(b);
brentmethoderror(count)=abs((b-z)/z)*100;
end
disp("Part B and D");
disp("Brent's Method");
disp("Number of Iterations:");
disp(count);
disp("Root located at x= ");
disp(b);
disp("Fzero Command");
disp("Root located at x= ");
disp(z);
%% Part B continued (Secant Method)
% The Secant method was then used to compare to the Brent Method.
%establish upper and lower bounds
low=-1;
up=1;
% Set Constants and Equations
f=@(x) -12- 21*x+x.^2-2.4*x.^3; %governing equation
iter=200000000; %iteration limit
tol=0.1*10^(-12); %tolerance value taken from textbook
secantcount=0; %set counter to enter loop
secanterror = 100; %initialize error
SEvec=zeros(1, 100); %space for error vector
secantx=0;
while (1)
secantcount=secantcount+1;
if (f(low)==f(up))
break;
end
secantx=up-((f(up))/((f(up)-f(low))))*(up-low); %generate new approx.
secanterror=abs((secantx-z)./z)*100;
SEvec(secantcount)=secanterror;
if (secantcount > iter) %if max iterations reached, break
disp('Max number of iterations has been exceeded');
break
end
if (secanterror < tol) %if the error tolerance met, break
break;
else
low=up;
up=secantx;
end
end
disp("Part B");
disp("Secant Method");
disp("Number of Iterations");
disp(secantcount);
disp("Root located at x= ");
disp(secantx);
% display error from both methods
xaxis=1:1:7;
figure(3);
hold all;
plot(xaxis, brentmethoderror(1:1:7), '-o', 'Color', 'r');
plot(xaxis, SEvec(1:1:7), '-o', 'Color', 'b');
title("Error for Secant and Brent's Method",'Fontsize', 18);
xlabel('Number of Iterations','Fontsize', 15);
ylabel('True Error %','Fontsize', 15);
legend( "Brent's Method Error","Secant Method Error","Fontsize", 15);
grid on;
hold off;
%% Part C and D (Catenary Cable)
%Set Constant and Equation
w=10;
y0=5;
y=15;
x=50;
tolerance=0.001;
eps=2.2204460*10^(-16);
func= @(Ta) ((Ta./w).*cosh((w./Ta).*x)+y0-(Ta./w))-y; %governing equation
%Part D (Fzero Command Method)
x0=[1000 1500]; % set upper and lower bounds
fzd= fzero(func, x0); %analytical root
disp("Part C");
disp ('Fzero Command');
disp("Root located at x= ");
disp(fzd);
%Set upper and lower bounds from brent method
a=1500;
b=1000;
%evaluate func at a and b
fa=func(a);
fb=func(b);
c=a; %new point equal to a
fc=fa; %fc same value as fa
d=b-a;
e=d; %another point equal to e
CCcounter=0; %counter for Catenery Cable
CCerror=zeros(1, 100); %space for error vector
while (1)
CCcounter=CCcounter+1;
if fb==0
break;
end
if (sign(fa) == sign(fb))
a=c;
fa=fc;
d=b-c;
e=d;
end
if(abs(fa)<abs(fb))% x values then switched to the new points made (a and b)
c=b;
b=a;
a=c;
fc=fb;
fb=fa;
fa=fc;
end
m=0.5*(a-b);
tolerance=2*(eps)*max(abs(b),1);
if (abs(m)<tolerance || fb==0)
break;
end
if (abs(e)>=tolerance && abs(fc)>abs(fb)) %open method used for for speed. if unacceptable answer is generated, secant method is chosen
s=fb/fc;
if (a==c)
p=2*m*s;
q=1-s;
else %if secant method generates unacceptable asnwer, inverse quadratic interpolation (IQI) is then used
q=fc/fa;
r=fb/fa;
p=s*(2*m*q*(q-r)-(b-c)*(r-1));
q=(q-1)*(r-1)*(s-1);
end
if(p>0)
q=-q;
else
p=-p;
end
if ((2*p<3*m*q-abs(tolerance*q)) && (p < abs(0.5*e*q)))
e=d;
d=p/q;
else
d=m;
e=m;
end
else %if IQI generates unacceptable answer, bisection method is used
d=m;
e=m;
end
c=b;
fc=fb;
if(abs(d)>tolerance)
b=b+d;
else
b=b-sign(b-a)*tolerance;
end
fb=func(b);
CCerror(CCcounter)=(abs((fzd - b)/fzd)) * 100;
end
disp("Brent's Method Catenary Cable");
disp('Number of iterations');
disp(CCcounter);
disp("Root located at x= ");
disp(b);
figure(2); %tension function plot
hold all;
k=1000:5:1500;
plot (k, func(k),'Color', 'r');
plot(fzd, func(fzd), 'bo');
xlabel('TA [N]', 'Fontsize', 15);
ylabel('y(TA)', 'Fontsize', 15);
title('Analytical Solution (Cable)', 'Fontsize', 18);
grid on;
hold off;
figure(4); %error plot
plotx=1:1:7;
hold all;
plot(plotx, CCerror(plotx), '-o', 'Color', 'r');
xlabel('Iterations','Fontsize', 15);
ylabel('True Error %','Fontsize', 15);
title("Brent's Method Error for Catenary Cable",'Fontsize', 18);
grid on;
hold off;