کسی بلده انتگرال sinx+x^2 از -1 تا 3 را به روش ذوزنقه ای در مطلب محاسبه کنه؟
اگر برنامه رو برای من ایمیل کنه ممنون میشم.shookhi10@yahoo.com
اگر برنامه رو برای من ایمیل کنه ممنون میشم.shookhi10@yahoo.com
function z = trapint(x, y) % % ------------------------------------------------------ % Reza Sepas Yar % http://www.avr.ir % ------------------------------------------------------ % % Trapezoidal integration of (x,y) data % UNEQUALLY SPACED DATA % % USAGE: z = trapint(x,y) % % input: x = vector of x values % y = vector of y values % % output: z = value of integral n = length(x); isum = 0; for i = 1:n-1 isum = isum + (y(i) + y(i+1))*(x(i+1) - x(i))/2; end z = isum;
>> x = [0:0.25:1]; >> y = 1./sqrt(1+x.^2); >> z = trapint(x, y) z = 0.8795
% % ------------------------------------------------------ % Reza Sepas Yar % http://www.avr.ir % ------------------------------------------------------ % % f is a function to integrating % a, the lower limit of integration % b, the upper limit of integration % n, the maximum number of Gauss points, 1-10 % Example1: % f= @(x) (x.*log(x))./(1+x); % a=0; % b=1; % n=2; % Example2: % f = @(x) 1./sqrt(2+(cos(x)).^2); % a=0; % b=1; % n=2; % Example3 % f = @(x) sin(x)./x; % a = 0 + 1e-10; % b = pi/2 + 1e-10; % n = 3; % Example3 % f = @(x) sin(x)./x; % a = 0 + 1e-10; % b = pi/2 + 1e-10; % n = 5; % Example4 % f = @(x) 1./sqrt(1+(sin(x)).^2); % a = 0; % b = 1; % n = 5; % Example4 % f = @(x) (x.*log(x))./(1+x); % a = 0; % b = 1; % n = 5; % Example4 % f = @(x) sin(x).^2; % a = 0; % b = pi/2; % n = 5; % Direct Way: quad(f, a, b) %********************************************************************** % Displays what inputs are being used disp(sprintf('\n\n********************************Input Data**********************************\n')) disp(sprintf(' f(x), integrand function')) disp(sprintf(' a = %g, lower limit of integration ',a)) disp(sprintf(' b = %g, upper limit of integration ',b)) disp(sprintf(' n = %g, number of Gauss Points, 1-10',n)) format long % Setup Gaussian Coefficients and Evaluation Points x_values = zeros(10,10) ; c_values = zeros(10,10) ; i=1 ; x_values(i,1) = 0.0 ; c_values(i,1) = 2.0 ; i=2 ; x_values(i,1) = -0.5773502691896260 ; x_values(i,2) = -x_values(i,1) ; c_values(i,1) = 1.0 ; c_values(i,2) = 1.0 ; i=3 ; x_values(i,1) = -0.7745966692414830 ; x_values(i,2) = 0.0 ; x_values(i,3) = -x_values(i,1) ; c_values(i,1) = 0.5555555555555560 ; c_values(i,2) = 0.8888888888888890 ; c_values(i,3) =c_values(i,1) ; i=4 ; x_values(i,1) = -0.8611363115940530 ; x_values(i,2) = -0.3399810435848560 ; x_values(i,3) = -x_values(i,2) ; x_values(i,4) = -x_values(i,1) ; c_values(i,1) = 0.3478548451374540 ; c_values(i,2) = 0.6521451548625460 ; c_values(i,3) =c_values(i,2) ; c_values(i,4) =c_values(i,1) ; i=5 ; x_values(i,1) = 0.9061798459386640 ; x_values(i,2) = 0.5384693101056830 ; x_values(i,3) = 0.0 ; x_values(i,4) = -x_values(i,2) ; x_values(i,5) = -x_values(i,1) ; c_values(i,1) = 0.2369268850561890 ; c_values(i,2) = 0.4786386704993660 ; c_values(i,3) = 0.5688888888890000 ; c_values(i,4) =c_values(i,2) ; c_values(i,5) =c_values(i,1) ; i=6 ; x_values(i,1) = -.9324695142032520 ; x_values(i,2) = -.6612093864662650 ; x_values(i,3) = -.2386191860831970 ; x_values(i,4) = -x_values(i,3) ; x_values(i,5) = -x_values(i,2) ; x_values(i,6) = -x_values(i,1) ; c_values(i,1) = 0.1713244923791700 ; c_values(i,2) = 0.3607615730481390 ; c_values(i,3) = 0.4679139345726910 ; c_values(i,4) =c_values(i,3) ; c_values(i,5) =c_values(i,2) ; c_values(i,6) =c_values(i,1) ; i=7 ; x_values(i,1) = -0.9491079123427590 ; x_values(i,2) = -0.7415311855993940 ; x_values(i,3) = -0.4058451513773970 ; x_values(i,4) = 0.0 ; x_values(i,5) = -x_values(i,3) ; x_values(i,6) = -x_values(i,2) ; x_values(i,7) = -x_values(i,1) ; c_values(i,1) = 0.1294849661688700 ; c_values(i,2) = 0.2797053914892770 ; c_values(i,3) = 0.3818300505051190 ; c_values(i,4) = 0.4179591836734690 ; c_values(i,5) =c_values(i,3) ; c_values(i,6) =c_values(i,2) ; c_values(i,7) =c_values(i,1) ; i=8 ; x_values(i,1) = -0.9602898564975360 ; x_values(i,2) = -0.7966664774136270 ; x_values(i,3) = -0.5255324099163290 ; x_values(i,4) = -0.1834346424956500 ; x_values(i,5) = -x_values(i,4) ; x_values(i,6) = -x_values(i,3) ; x_values(i,7) = -x_values(i,2) ; x_values(i,8) = -x_values(i,1) ; c_values(i,1) = 0.1012285362903760 ; c_values(i,2) = 0.2223810344533740 ; c_values(i,3) = 0.3137066458778870 ; c_values(i,4) = 0.3626837833783620 ; c_values(i,5) =c_values(i,4) ; c_values(i,6) =c_values(i,3) ; c_values(i,7) =c_values(i,2) ; c_values(i,8) =c_values(i,1) ; i=9 ; x_values(i,1) = -0.9681602395076260 ; x_values(i,2) = -0.8360311073266360 ; x_values(i,4) = -0.6133714327005900 ; x_values(i,4) = -0.3242534234038090 ; x_values(i,5) = 0.0 ; x_values(i,6) = -x_values(i,4) ; x_values(i,7) = -x_values(i,3) ; x_values(i,8) = -x_values(i,2) ; x_values(i,9) = -x_values(i,1) ; c_values(i,1) = 0.0812743883615740 ; c_values(i,2) = 0.1806481606948570 ; c_values(i,3) = 0.2606106964029350 ; c_values(i,4) = 0.3123470770400030 ; c_values(i,5) = 0.3302393550012600 ; c_values(i,6) =c_values(i,4) ; c_values(i,7) =c_values(i,3) ; c_values(i,8) =c_values(i,2) ; c_values(i,9) =c_values(i,1) ; i=10 ; x_values(i,1) = -0.9739065285171720 ; x_values(i,2) = -0.8650633666889850 ; x_values(i,3) = -0.6794095682990240 ; x_values(i,4) = -0.4333953941292470 ; x_values(i,5) = -0.1488743389816310 ; x_values(i,6) = -x_values(i,5) ; x_values(i,7) = -x_values(i,4) ; x_values(i,8) = -x_values(i,3) ; x_values(i,9) = -x_values(i,2) ; x_values(i,10) = -x_values(i,1) ; c_values(i,1) = 0.0666713443086880 ; c_values(i,2) = 0.1494513491505810 ; c_values(i,3) = 0.2190863625159820 ; c_values(i,4) = 0.2692667193099960 ; c_values(i,5) = 0.2955242247147530 ; c_values(i,6) = c_values(i,5) ; c_values(i,7) = c_values(i,4) ; c_values(i,8) = c_values(i,3) ; c_values(i,9) = c_values(i,2) ; c_values(i,10) = c_values(i,1) ; % simulation of the steps needed to perform the method begins here % disp(sprintf('\n*******************************Simulation*********************************\n')) % First do a lookup of the weights and evaluation locations disp('1) Lookup Gauss point constants and evaluation points from the table.') disp(' ') disp(sprintf(' Evaluation Locations')) % Used to display the gauss points and their weights for i=1:n disp(sprintf(' x%i=%1.15f',i,x_values(n,i))) end disp(' ') disp(sprintf(' Evaluation Constants')) for i=1:n disp(sprintf(' c%i=%1.15f',i,c_values(n,i))) end % Now transform the evaluation points to fit your interval % disp(' ') % disp('2) The evaluation points are definted for the interval [-1,1]. Transform them') % disp(' to fit [a,b]. This is done by the following formula') % disp(' ') % disp(' xnew = (b-a)/2 * x + (b+a)/2') % disp(' ') % disp(sprintf(' Transformed Evaluation Locations')) % Transform the locations for i=1:n xv(i)=(b-a)/2*x_values(n,i)+(b+a)/2 ; %disp(sprintf(' x%i=%1.15f',i,xv(i))); end % Apply the Gaussian procedure disp(' ') disp('2) The approximation is given as') disp(' ') disp(' c(i) * f(x(i))') approx = 0 ; disp(' ') % approximation via Gaussian integration for i=1:n-1 approx = approx + c_values(n,i)*f(xv(i)) ; disp(sprintf(' %1.15f * %1.15f = %1.15f',c_values(n,i),f(xv(i)), c_values(n,i)*f(xv(i) ))) end approx = approx + c_values(n,n)*f(xv(n)) ; disp(sprintf(' + %1.15f * %1.15f = %1.15f',c_values(n,n),f(xv(n)), c_values(n,n)*f(xv(n) ))) disp(sprintf(' ----------------------------')) disp(sprintf(' %1.15f',approx)) % the result must be scaled to fit the interval also disp(' ') disp('3) The approximation must also be transformed to fit the [a,b] interval.') disp(' ') disp(' Answer = (b-a)/2 * sum( c(i)* f(x(i)) ) ') disp(' ') disp(sprintf(' = %1.15f * %1.15f',(b-a)/2,approx)) % also transform the approximation approx = (b-a)/2 * approx ; disp(sprintf(' = %1.15f',approx)) % Display results section disp(sprintf('\n\n**************************Results****************************')) % The following finds what is called the 'Exact' solution exact = quad(f,a,b) ; disp(sprintf('\n Approximate = %1.15f',approx)) disp(sprintf(' Exact = %1.15f',exact))
% ------------------------------------------------------ % Reza Sepas Yar % http://www.avr.ir % ------------------------------------------------------ f= x.*log(x)./(1+x); % a, the lower limit of integration a=1e-10 ; % b, the upper limit of integration b=1 ; % n, the maximum number of Gauss points, 1-10 n=2 ; %********************************************************************** % Displays what inputs are being used disp(sprintf('\n\n********************************Input Data**********************************\n')) disp(sprintf(' f(x), integrand function')) disp(sprintf(' a = %g, lower limit of integration ',a)) disp(sprintf(' b = %g, upper limit of integration ',b)) disp(sprintf(' n = %g, number of Gauss Points, 1-10',n)) format short g % Setup Gaussian Coefficients and Evaluation Points x_values = zeros(10,10) ; c_values = zeros(10,10) ; i=1 ; x_values(i,1) = 0.0 ; c_values(i,1) = 2.0 ; i=2 ; x_values(i,1) = -0.5773502691896260 ; x_values(i,2) = -x_values(i,1) ; c_values(i,1) = 1.0 ; c_values(i,2) = 1.0 ; i=3 ; x_values(i,1) = -0.7745966692414830 ; x_values(i,2) = 0.0 ; x_values(i,3) = -x_values(i,1) ; c_values(i,1) = 0.5555555555555560 ; c_values(i,2) = 0.8888888888888890 ; c_values(i,3) =c_values(i,1) ; i=4 ; x_values(i,1) = -0.8611363115940530 ; x_values(i,2) = -0.3399810435848560 ; x_values(i,3) = -x_values(i,2) ; x_values(i,4) = -x_values(i,1) ; c_values(i,1) = 0.3478548451374540 ; c_values(i,2) = 0.6521451548625460 ; c_values(i,3) =c_values(i,2) ; c_values(i,4) =c_values(i,1) ; i=5 ; x_values(i,1) = -0.9061798459386640 ; x_values(i,2) = -0.5384693101056830 ; x_values(i,3) = 0.0 ; x_values(i,4) = -x_values(i,2) ; x_values(i,5) = -x_values(i,1) ; c_values(i,1) = 0.2369368850561890 ; c_values(i,2) = 0.4786386704993660 ; c_values(i,3) = 0.5688888888888890 ; c_values(i,4) =c_values(i,2) ; c_values(i,5) =c_values(i,1) ; i=6 ; x_values(i,1) = -.9324695142032520 ; x_values(i,2) = -.6612093864662650 ; x_values(i,3) = -.2386191860831970 ; x_values(i,4) = -x_values(i,3) ; x_values(i,5) = -x_values(i,2) ; x_values(i,6) = -x_values(i,1) ; c_values(i,1) = 0.1713244923791700 ; c_values(i,2) = 0.3607615730481390 ; c_values(i,3) = 0.4679139345726910 ; c_values(i,4) =c_values(i,3) ; c_values(i,5) =c_values(i,2) ; c_values(i,6) =c_values(i,1) ; i=7 ; x_values(i,1) = -0.9491079123427590 ; x_values(i,2) = -0.7415311855993940 ; x_values(i,3) = -0.4058451513773970 ; x_values(i,4) = 0.0 ; x_values(i,5) = -x_values(i,3) ; x_values(i,6) = -x_values(i,2) ; x_values(i,7) = -x_values(i,1) ; c_values(i,1) = 0.1294849661688700 ; c_values(i,2) = 0.2797053914892770 ; c_values(i,3) = 0.3818300505051190 ; c_values(i,4) = 0.4179591836734690 ; c_values(i,5) =c_values(i,3) ; c_values(i,6) =c_values(i,2) ; c_values(i,7) =c_values(i,1) ; i=8 ; x_values(i,1) = -0.9602898564975360 ; x_values(i,2) = -0.7966664774136270 ; x_values(i,3) = -0.5255324099163290 ; x_values(i,4) = -0.1834346424956500 ; x_values(i,5) = -x_values(i,4) ; x_values(i,6) = -x_values(i,3) ; x_values(i,7) = -x_values(i,2) ; x_values(i,8) = -x_values(i,1) ; c_values(i,1) = 0.1012285362903760 ; c_values(i,2) = 0.2223810344533740 ; c_values(i,3) = 0.3137066458778870 ; c_values(i,4) = 0.3626837833783620 ; c_values(i,5) =c_values(i,4) ; c_values(i,6) =c_values(i,3) ; c_values(i,7) =c_values(i,2) ; c_values(i,8) =c_values(i,1) ; i=9 ; x_values(i,1) = -0.9681602395076260 ; x_values(i,2) = -0.8360311073266360 ; x_values(i,4) = -0.6133714327005900 ; x_values(i,4) = -0.3242534234038090 ; x_values(i,5) = 0.0 ; x_values(i,6) = -x_values(i,4) ; x_values(i,7) = -x_values(i,3) ; x_values(i,8) = -x_values(i,2) ; x_values(i,9) = -x_values(i,1) ; c_values(i,1) = 0.0812743883615740 ; c_values(i,2) = 0.1806481606948570 ; c_values(i,3) = 0.2606106964029350 ; c_values(i,4) = 0.3123470770400030 ; c_values(i,5) = 0.3302393550012600 ; c_values(i,6) =c_values(i,4) ; c_values(i,7) =c_values(i,3) ; c_values(i,8) =c_values(i,2) ; c_values(i,9) =c_values(i,1) ; i=10 ; x_values(i,1) = -0.9739065285171720 ; x_values(i,2) = -0.8650633666889850 ; x_values(i,3) = -0.6794095682990240 ; x_values(i,4) = -0.4333953941292470 ; x_values(i,5) = -0.1488743389816310 ; x_values(i,6) = -x_values(i,5) ; x_values(i,7) = -x_values(i,4) ; x_values(i,8) = -x_values(i,3) ; x_values(i,9) = -x_values(i,2) ; x_values(i,10) = -x_values(i,1) ; c_values(i,1) = 0.0666713443086880 ; c_values(i,2) = 0.1494513491505810 ; c_values(i,3) = 0.2190863625159820 ; c_values(i,4) = 0.2692667193099960 ; c_values(i,5) = 0.2955242247147530 ; c_values(i,6) = c_values(i,5) ; c_values(i,7) = c_values(i,4) ; c_values(i,8) = c_values(i,3) ; c_values(i,9) = c_values(i,2) ; c_values(i,10) = c_values(i,1) ; % simulation of the steps needed to perform the method begins here disp(sprintf('\n*******************************Simulation*********************************\n')) % First do a lookup of the weights and evaluation locations disp('1) Lookup Gauss point constants and evaluation points from the table.') disp(' ') disp(sprintf(' Evaluation Locations')) % Used to display the gauss points and their weights for i=1:n disp(sprintf(' x%i=%1.15f',i,x_values(n,i))) end disp(' ') disp(sprintf(' Evaluation Constants')) for i=1:n disp(sprintf(' c%i=%1.15f',i,c_values(n,i))) end % Now transform the evaluation points to fit your interval disp(' ') disp('2) The evaluation points are definted for the interval [-1,1]. Transform them') disp(' to fit [a,b]. This is done by the following formula') disp(' ') disp(' xnew = (b-a)/2 * x + (b+a)/2') disp(' ') disp(sprintf(' Transformed Evaluation Locations')) % Transform the locations for i=1:n xv(i)=(b-a)/2*x_values(n,i)+(b+a)/2 ; disp(sprintf(' x%i=%1.15f',i,xv(i))) end % Apply the Gaussian procedure disp(' ') disp('3) The approximation is given as') disp(' ') disp(' approx = sum(i=1,n) c(i)*f(x(i))') approx = 0 ; disp(' ') % approximation via Gaussian integration for i=1:n-1 approx = approx + c_values(n,i)*f(xv(i)) ; disp(sprintf(' %1.15f * %g',c_values(n,i),f(xv(i)))) end approx = approx + c_values(n,n)*f(xv(n)) ; disp(sprintf(' + %1.15f * %g',c_values(n,n),f(xv(n)))) disp(sprintf(' ----------------------------')) disp(sprintf(' %g',approx)) % the result must be scaled to fit the interval also disp(' ') disp('4) The approximation must also be transformed to fit the [a,b] interval.') disp(' ') disp(' approx = (b-a)/2 * approx') disp(sprintf(' = %g * %g',(b-a)/2,approx)) % also transform the approximation approx = (b-a)/2 * approx ; disp(sprintf(' = %g',approx)) % Display results section disp(sprintf('\n\n**************************Results****************************')) % The following finds what is called the 'Exact' solution exact = quad(f,a,b) ; disp(sprintf('\n Approximate = %g',approx)) disp(sprintf(' Exact = %g',exact)) disp(sprintf('\n True Error = Exact - Approximate')) disp(sprintf(' = %g - %g',exact,approx)) disp(sprintf(' = %g',exact-approx)) disp(sprintf('\n Absolute Relative True Error Percentage')) disp(sprintf(' = | ( Exact - Approximate ) / Exact | * 100')) disp(sprintf(' = | %g / %g | * 100',exact-approx,exact)) disp(sprintf(' = %g\n\n',abs( (exact-approx)/exact )*100))
function area = simpson(h,y) % % ------------------------------------------------------ % Reza Sepas Yar % http://www.avr.ir % ------------------------------------------------------ % % area = simpson(h,y) % % Example: % x = [0:0.25:1] % y = 1/sqrt(1+x.^2) % simpson(0.25,y) % aiyi = % 1.000000000000000 3.880570000581328 1.788854381999832 3.200000000000000 0.707106781186547 % ans = % 0.881377596980642 % % Function Simpson computes the area under a function defined by % an odd number of ordinates in vector y spaced s units apart % by Simpson's First Rule. % n=length(y); aiyi(1) = y(1); aiyi(n) = y(n); for i=2:(n-1) if (mod(i,2) == 0) aiyi(i) = 4 * y(i); else aiyi(i) = 2 * y(i); end end aiyi if (mod(n,2) == 0) disp('Y must have odd element!'); area = 0; else area = y(1)+y(n); for i=2:(n-1) if (mod(i,2) == 0) area = area + 4 * y(i); else area = area + 2 *y(i); end end area = h*area/3; end end % Other Examples: % x = 0 + 1e-10:pi/16:pi/2+1e-10 % x = % Columns 1 through 6 % 0.000000000100000 0.196349540949362 0.392699081798724 0.589048622648086 0.785398163497448 0.981747704346810 % Columns 7 through 9 % 1.178097245196172 1.374446786045535 1.570796326894897 % y = sin(x)./x % y = % Columns 1 through 6 % 1.000000000000000 0.993586851137686 0.974495358391543 0.943165320725420 0.900316316132506 0.846927992473694 % Columns 7 through 9 % 0.784213303542454 0.713585487907166 0.636619772327053 % simpson(pi/16,y) % aiyi = % Columns 1 through 6 % 1.000000000000000 3.974347404550744 1.948990716783087 3.772661282901681 1.800632632265012 3.387711969894778 % Columns 7 through 9 % 1.568426607084908 2.854341951628665 0.636619772327053 % ans = % 1.370764076042494 %==================================================================== % x = [0+1e-10:1/8:1+1e-10] % x = % Columns 1 through 6 % 0.000000000100000 0.125000000100000 0.250000000100000 0.375000000100000 0.500000000100000 0.625000000100000 % Columns 7 through 9 % 0.750000000100000 0.875000000100000 1.000000000100000 % y = (x.*log(x))./(1+x) % y = % Columns 1 through 6 % -0.000000002302585 -0.231049060262061 -0.277258872232701 -0.267498887164168 -0.231049060150788 -0.180770626589236 % Columns 7 through 9 % -0.123292316717300 -0.062314649841909 0.000000000050000 % simpson(1/8,y) % aiyi = % Columns 1 through 6 % -0.000000002302585 -0.924196241048244 -0.554517744465402 -1.069995548656670 -0.462098120301577 -0.723082506356943 % Columns 7 through 9 % -0.246584633434600 -0.249258599367635 0.000000000050000 % ans = % -0.176238891495152 %==================================================================== % x = [0:pi/16:pi/2] % x = % Columns 1 through 6 % 0 0.196349540849362 0.392699081698724 0.589048622548086 0.785398163397448 0.981747704246810 % Columns 7 through 9 % 1.178097245096172 1.374446785945535 1.570796326794897 % y = sin(x.^2) % y = % Columns 1 through 6 % 0 0.038543592357938 0.153602060372138 0.340057724740335 0.578468789354558 0.821381311850125 % Columns 7 through 9 % 0.983323424736072 0.949766418160832 0.624265952639699 % simpson(pi/16,y) % aiyi = % Columns 1 through 6 % 0 0.154174369431752 0.307204120744276 1.360230898961340 1.156937578709116 3.285525247400499 % Columns 7 through 9 % 1.966646849472143 3.799065672643329 0.624265952639699 % ans = % 0.828205680955492 %==================================================================== % x = 0:1/8:1; % y = 1./sqrt(2+(cos(x)).^2); % simpson(1/8,y) % aiyi = % 0.57735 2.3154 1.1667 2.3628 1.2017 2.4536 1.2561 2.5762 0.66054 % ans = % 0.6071 %
دیدگاه