Root Finding (Encontrando as raízes de equações)

Prévia do material em texto

Name: Paulo Matos
MECH 2220 Lecture Homework 3
Lecture HW 03 Root Finding	
		Assigned: February 17, 2015
		Due: February 24, 2015 (11:00 a.m.)
�
Problem Statement #1:
Write a MATLAB script that will find the roots of a given equations using the bisection method. Format your output to look similar to the examples given. You should write your output to a file. Set the maximum number of iterations to 100 and the tolerance to 0.000005. 
Use a tolerance definition of 
 are the interval endpoints for the nth iteration
Use the script to solve the given equations. Use the interval endpoints as your initial a and b. 
Eq1	f(x) = ex - (x2 + 4) = 0 on [2,3] 
Eq2	f(x) = x - 2-x = 0 on [1/3,1] 
Eq3	f(x) = x3 - 1.8999 x2 + 1.5796 x - 2.1195 = 0 on [1,2]. 
Code for problem #1:
Function 1
function [out]=f1(x)
 out=(exp(x)-((x^2)+4));
end
Script 1
clc;clear;
%opening a file
file = fopen('hw3_eq1.txt','w');
%defining bounds
a(1)=2;
b(1)=3;
%checking bounds
if f1(a(1))*f1(b(1))>=0
 error('bad bounds');
end
%using the bissect interval method
n=1;
tolerance=10000;
while n<100 && tolerance>=0.000005
 if f1(a(n))*f1(b(n))<=0
 p(n)=(a(n)+b(n))/2;
 end
 if f1(a(n))*f1(p(n))<0
 a(n+1)=a(n);
 b(n+1)=p(n);
 elseif f1(p(n))*f1(b(n))<0
 a(n+1)=p(n);
 b(n+1)=b(n);
 elseif f1(p(n))==0
 fprintf('root is found');
 else
 error('I don''t understand');
 end
 tolerance=abs((b(n)-a(n))/2);
 n=n+1;
end
%printing
if length(p) >=1
fprintf(file,'The equation''s root for f(x)=(exp(x)-((x^2)+4))=0 on [%.2f,%.2f]\n',a(1),b(1));
fprintf(file,'is %.6f where f(%.6f)=%.6f\n',p(n-1),p(n-1),f1(p(n-1)));
fprintf(file,'There are %.0f iterations for the bisection, and they are:\n',(n-1));
fprintf(file,' n\t\t\t a\t\t\t b\t\t\t p\t\t\t f(p)\t\t\n');
z=1;
while z<=(n-1)
fprintf(file,'%3.0f \t\t %.6f\t %.6f\t %.6f\t %.6f\t\n', z, a(z), b(z),p(z), f1(p(z)));
z=z+1;
end
end
fclose(file);
Function 2
function [out] = f2(x)
out =(x-2^(-x));
end
Script 2
clc;clear;
%opening a file
file = fopen('hw3_eq2.txt','w');
%defining bounds
a(1)=(1/3);
b(1)=1;
%checking bounds
if f2(a(1))*f2(b(1))>=0
 error('bad bounds');
end
%using the bissect interval method
n=1;
tolerance=10000;
while n<100 && tolerance>=0.000005
 if f2(a(n))*f2(b(n))<=0
 p(n)=(a(n)+b(n))/2;
 end
 if f2(a(n))*f2(p(n))<0
 a(n+1)=a(n);
 b(n+1)=p(n);
 elseif f2(p(n))*f2(b(n))<0
 a(n+1)=p(n);
 b(n+1)=b(n);
 elseif f2(p(n))==0
 fprintf('root is found');
 else
 error('I don''t understand');
 end
 tolerance=abs((b(n)-a(n))/2);
 n=n+1;
end
%printing
if length(p) >=1
fprintf(file,'The equation''s root for f(x)=(x-2^(-x)) on [%.2f,%.2f]\n',a(1),b(1));
fprintf(file,'is %.6f where f(%.6f)=%.6f\n',p(n-1),p(n-1),f2(p(n-1)));
fprintf(file,'There are %.0f iterations for the bisection, and they are:\n',(n-1));
fprintf(file,' n\t\t\t a\t\t\t b\t\t\t p\t\t\t f(p)\t\t\n');
z=1;
while z<=(n-1)
fprintf(file,'%3.0f \t\t %.6f\t %.6f\t %.6f\t %.6f\t\n', z, a(z), b(z),p(z), f2(p(z)));
z=z+1;
end
end
fclose(file);
Function 3
function [out] = f3(x)
 out =(x^3 - 1.8999*x^2 + 1.5796*x - 2.1195);
end
Script 3
clc;clear;
%opening a file
file = fopen('hw3_eq3.txt','w');
%defining bounds
a(1)=1;
b(1)=2;
%checking bounds
if f3(a(1))*f3(b(1))>=0
 error('bad bounds');
end
%using the bissect interval method
n=1;
tolerance=10000;
while n<100 && tolerance>=0.000005
 if f3(a(n))*f3(b(n))<=0
 p(n)=(a(n)+b(n))/2;
 end
 if f3(a(n))*f3(p(n))<0
 a(n+1)=a(n);
 b(n+1)=p(n);
 elseif f3(p(n))*f3(b(n))<0
 a(n+1)=p(n);
 b(n+1)=b(n);
 elseif f3(p(n))==0
 fprintf('root is found');
 else
 error('I don''t understand');
 end
 tolerance=abs((b(n)-a(n))/2);
 n=n+1;
end
%printing
if length(p) >=1
fprintf(file,'The equation''s root for f(x)=(x^3 - 1.8999*x^2 + 1.5796*x - 2.1195) on [%.2f,%.2f]\n',a(1),b(1));
fprintf(file,'is %.6f where f(%.6f)=%.6f\n',p(n-1),p(n-1),f3(p(n-1)));
fprintf(file,'There are %.0f iterations for the bisection, and they are:\n',(n-1));
fprintf(file,' n\t\t\t a\t\t\t b\t\t\t p\t\t\t f(p)\t\t\n');
z=1;
while z<=(n-1)
fprintf(file,'%3.0f \t\t %.6f\t %.6f\t %.6f\t %.6f\t\n', z, a(z), b(z),p(z), f3(p(z)));
z=z+1;
end
end
fclose(file);
Solution for problem #1:
Output 1
The equation's root for f(x)=(exp(x)-((x^2)+4))=0 on [2.00,3.00]
is 2.158726 where f(2.158726)=-0.000001
There are 18 iterations for the bisection, and they are:
 n a b p f(p) 
 1 2.000000 3.000000 2.500000 1.932494 
 2 2.000000 2.500000 2.250000 0.425236 
 3 2.000000 2.250000 2.125000 -0.142728 
 4 2.125000 2.250000 2.187500 0.127747 
 5 2.125000 2.187500 2.156250 -0.010732 
 6 2.156250 2.187500 2.171875 0.057680 
 7 2.156250 2.171875 2.164063 0.023269 
 8 2.156250 2.164063 2.160156 0.006218 
 9 2.156250 2.160156 2.158203 -0.002270 
 10 2.158203 2.160156 2.159180 0.001971 
 11 2.158203 2.159180 2.158691 -0.000151 
 12 2.158691 2.159180 2.158936 0.000910 
 13 2.158691 2.158936 2.158813 0.000380 
 14 2.158691 2.158813 2.158752 0.000115 
 15 2.158691 2.158752 2.158722 -0.000018 
 16 2.158722 2.158752 2.158737 0.000048 
 17 2.158722 2.158737 2.158730 0.000015 
 18 2.158722 2.158730 2.158726 -0.000001 
Output 2
The equation's root for f(x)=(x-2^(-x)) on [0.33,1.00]
is 0.641187 where f(0.641187)=0.000002
There are 18 iterations for the bisection, and they are:
 n a b p f(p) 
 1 0.333333 1.000000 0.666667 0.036706 
 2 0.333333 0.666667 0.500000 -0.207107 
 3 0.500000 0.666667 0.583333 -0.084087 
 4 0.583333 0.666667 0.625000 -0.023420 
 5 0.625000 0.666667 0.645833 0.006710 
 6 0.625000 0.645833 0.635417 -0.008338 
 7 0.635417 0.645833 0.640625 -0.000810 
 8 0.640625 0.645833 0.643229 0.002951 
 9 0.640625 0.643229 0.641927 0.001071 
 10 0.640625 0.641927 0.641276 0.000130 
 11 0.640625 0.641276 0.640951 -0.000340 
 12 0.640951 0.641276 0.641113 -0.000105 
 13 0.641113 0.641276 0.641195 0.000013 
 14 0.641113 0.641195 0.641154 -0.000046 
 15 0.641154 0.641195 0.641174 -0.000017 
 16 0.641174 0.641195 0.641184 -0.000002 
 17 0.641184 0.641195 0.641190 0.000006 
 18 0.641184 0.641190 0.641187 0.000002 
Output 3
The equation's root for f(x)=(x^3 - 1.8999*x^2 + 1.5796*x - 2.1195) on [1.00,2.00]
is 1.703129 where f(1.703129)=-0.000002
There are 18 iterations for the bisection, and they are:
 n a b p f(p) 
 1 1.000000 2.000000 1.500000 -0.649875 
 2 1.500000 2.000000 1.750000 0.185731 
 3 1.500000 1.750000 1.625000 -0.278558 
 4 1.625000 1.750000 1.687500 -0.058767 
 5 1.687500 1.750000 1.718750 0.060302 
 6 1.687500 1.718750 1.703125 -0.000016 
 7 1.703125 1.718750 1.710938 0.029946 
 8 1.703125 1.710938 1.707031 0.014916 
 9 1.703125 1.7070311.705078 0.007437 
 10 1.703125 1.705078 1.704102 0.003708 
 11 1.703125 1.704102 1.703613 0.001845 
 12 1.703125 1.703613 1.703369 0.000914 
 13 1.703125 1.703369 1.703247 0.000449 
 14 1.703125 1.703247 1.703186 0.000216 
 15 1.703125 1.703186 1.703156 0.000100 
 16 1.703125 1.703156 1.703140 0.000042 
 17 1.703125 1.703140 1.703133 0.000013 
 18 1.703125 1.703133 1.703129 -0.000002 
Problem Statement #2:
Do problem 5.9 (5.8 in the third edition) in the text
Code for problem #2:
Function
function [out] = eqc(TC,con)
TK=TC+273.15;
out=-log(con)-139.3441+(1.575701e5)/TK-(6.642308e7)/(TK^2)+(1.243800e10)/(TK^3)-(8.621949e11)/(TK^4);
end
Script
clc;
clear;
%First Part
%determine the number of iterations
A(1)=0;
B(1)=40;
%error
error1=abs(B(1)-A(1))/0.05;
iterations=ceil(log2(error1));
fprintf('The number of iterations for 0.05 error is %.0f\n',iterations)
 
 
%Second Part
d=[8,10,12]; %concentrations
f=1;
fmax=length(d);
while f<=fmax
 m=1;
 err=inf;
 if eqc(A(1),d(f))*eqc(B(1),d(f))>=0
 error('bounds are bad')
 end
 n=1;
 while n<=iterations+1 && err>=0.05
 if eqc(A(n),d(f))*eqc(B(n),d(f))<=0;
 P(n)=(A(n)+B(n))/2;
 end
 if eqc(A(n),d(f))*eqc(P(n),d(f))<0
 A(n+1)=A(n);B(n+1)=P(n);
 elseif eqc(B(n),d(f))*eqc(P(n),d(f))<0
 B(n+1)=B(n); A(n+1)=P(n);
 elseif eqc(P(n),d(f))==0
 fprintf('found root')
 else
 error('something is wrong')
 end
 err=abs(B(n)-A(n))/2;
 n=n+1;
 end
 if length(P)>=1
 fprintf('In the interval[%.1f,%.1f], the solution for the equation is \n',A(1),B(1));
 fprintf('with a concentration of %.2f mg/L\n',d(f));
 fprintf('is %.4f degrees Celsius\n\n', P(n-1));
 fprintf('There are %.0f iterations for the bisection, and they are:\n\n',(n-1));
 fprintf(' n\t A\t\t B\t\t P\t\t eqc(P)\t\t\n');
 g=1;
 while g<=(n-1)
 fprintf('%3.0f \t %.6f\t %.6f\t %.6f\t %.6f\t\n', g, A(g), B(g), P(g),eqc(P(g),d(f)));
 g=g+1;
 end
 end
 f=f+1;
 fprintf('\n');
end
 Solution for problem #2:
The number of iterations for 0.05 error is 10
In the interval[0.0,40.0], the solution for the equation is 
with a concentration of 8.00 mg/L
is 26.7578 degrees Celsius
There are 10 iterations for the bisection, and they are:
 n	 A		 B		 P		 eqc(P)		
 1 	 0.000000	 40.000000	 20.000000	 0.128010	
 2 	 20.000000	 40.000000	 30.000000	 -0.056720	
 3 	 20.000000	 30.000000	 25.000000	 0.032411	
 4 	 25.000000	 30.000000	 27.500000	 -0.012874	
 5 	 25.000000	 27.500000	 26.250000	 0.009578	
 6 	 26.250000	 27.500000	 26.875000	 -0.001694	
 7 	 26.250000	 26.875000	 26.562500	 0.003930	
 8 	 26.562500	 26.875000	 26.718750	 0.001115	
 9 	 26.718750	 26.875000	 26.796875	 -0.000290	
 10 	 26.718750	 26.796875	 26.757813	 0.000412	
In the interval[0.0,40.0], the solution for the equation is 
with a concentration of 10.00 mg/L
is 15.3516 degrees Celsius
There are 10 iterations for the bisection, and they are:
 n	 A		 B		 P		 eqc(P)		
 1 	 0.000000	 40.000000	 20.000000	 -0.095133	
 2 	 0.000000	 20.000000	 10.000000	 0.121160	
 3 	 10.000000	 20.000000	 15.000000	 0.008361	
 4 	 15.000000	 20.000000	 17.500000	 -0.044462	
 5 	 15.000000	 17.500000	 16.250000	 -0.018331	
 6 	 15.000000	 16.250000	 15.625000	 -0.005056	
 7 	 15.000000	 15.625000	 15.312500	 0.001634	
 8 	 15.312500	 15.625000	 15.468750	 -0.001716	
 9 	 15.312500	 15.468750	 15.390625	 -0.000042	
 10 	 15.312500	 15.390625	 15.351563	 0.000796	
In the interval[0.0,40.0], the solution for the equation is 
with a concentration of 12.00 mg/L
is 7.4609 degrees Celsius
There are 10 iterations for the bisection, and they are:
 n	 A		 B		 P		 eqc(P)		
 1 	 0.000000	 40.000000	 20.000000	 -0.277455	
 2 	 0.000000	 20.000000	 10.000000	 -0.061161	
 3 	 0.000000	 10.000000	 5.000000	 0.062280	
 4 	 5.000000	 10.000000	 7.500000	 -0.000849	
 5 	 5.000000	 7.500000	 6.250000	 0.030354	
 6 	 6.250000	 7.500000	 6.875000	 0.014663	
 7 	 6.875000	 7.500000	 7.187500	 0.006885	
 8 	 7.187500	 7.500000	 7.343750	 0.003012	
 9 	 7.343750	 7.500000	 7.421875	 0.001080	
 10 	 7.421875	 7.500000	 7.460938	 0.000115	
 
Problem Statement #3:
Do problem 5.13 (same number for second and third edition) in the text.
Code for problem #3:
clear;clc;
E=50000;I=30000;w=2.5;L=600;
%derivative of the function
f=@(x) ((w)/(120*E*I*L))*(-(5*(x^4))+(2*(L^2)*(3*(x^2)))-((L^4)));
A=0;B=600;
tol=.00005;iter=100;
A(1)=A;B(1)=B;
if f(A(1))*f(B(1))>0
 error('Bad Bounds')
elseif abs(f(A))<10^(-14);
 r=A;
elseif abs(f(A))<10^(-14)
 r=B;
else
 P(1)=(A(1)+B(1))/2;
 n=1;
 err=(B(n)-A(n))/2;
 while n<=iter && abs(err)>tol
 if f(A(n))*f(P(n))<0
 A(n+1)=A(n);
 B(n+1)=P(n);
 P(n+1)=(A(n+1)+B(n+1))/2;
 elseif f(B(n))*f(P(n))<0
 A(n+1)=P(n);
 B(n+1)=B(n);
 P(n+1)=(A(n+1)+B(n+1))/2;
 elseif f(P(n))==0
 r=P(n);break
 else
 error('Something is wrong')
 end
 n=n+1;
 err=(B(n)-A(n))/2;
 end
 r=P(n);
end
x=r;
maximum_deflection=((w)/(120*E*I*L))*(-(x^5)+(2*(L^2)*(x^3))-((L^4)*x));
fprintf('The maximum deflection is %.4f cm and the it occurs at %.4f cm', maximum_deflection,r);
Solution for problem #3:
The maximum deflection is -0.5152 cm and the it occurs at 268.3282 cm
�
Sample Output for Equation 3:
Joe Ragan 
Bisection Routine
The root to the equation f= x^3-1.8999*x^2+1.5796*x-2.1195 =0 on [1,2] was found to be 1.703133 where f(1.703133)=0.000013.
The 17 iterations for the bisection routine were recorded as: 
n 	 a 		 b 		 p 		 f(p) 
1 	 1.000000 	 2.000000 	 1.500000 	 -0.649875 
2 	 1.500000 	 2.000000 	 1.750000 	 0.185731 
3 	 1.500000 	 1.750000 	 1.625000 	 -0.278558 
4 	 1.625000 	 1.750000 	 1.687500 	 -0.058767 
5 	 1.687500 	 1.750000 	 1.718750 	 0.060302 
6 	 1.687500 	 1.718750 	 1.703125 	 -0.000016 
7 	 1.703125 	 1.718750 	 1.710938 	 0.029946 
8 	 1.703125 	 1.710938 	 1.707031 	 0.014916 
9 	 1.703125 	 1.707031 	 1.705078 	 0.007437 
10 	 1.703125 	 1.705078 	 1.704102 	 0.003708 
11 	 1.703125 	 1.704102 	 1.703613 	 0.001845 
12 	 1.703125 	 1.703613 	 1.703369 	 0.000914 
13 	 1.703125 	 1.703369 	 1.703247 	 0.000449 
14 	 1.703125 	 1.703247 	 1.703186 	 0.000216 
15 	 1.703125 	 1.703186 	 1.703156 	 0.000100 
16 	 1.703125 	 1.703156 	 1.703140 	 0.000042 
17 	 1.703125 	 1.703140 	 1.703133 	 0.000013 
 
_1127540000.unknown