Two Incomplete Matlab Codes for Numerical Integration
Related to a Computer Project
Note: You need to modify the codes to get them work.
% integral1.m: comparison of numerical integration rules a = 0; b = pi/2; M = 10; for n = 1:M h = (b-a)/n; rl = 0; rm = 0; for i = 1:n rl = rl + sin(a+(i-1)*h); % left endpont rectangle rule rm = rm + sin( ? ); % midpoint rectangle rule end rl = h*rl; erl(n) = log(abs(rl-1)); rm = h*rm; erm(n) = log(abs(rm-1)); t = sin(a) + sin(b); % trapezoid rule for i = 1:n-1 t = t + ?; end t = 0.5*h*t; et(n) = log(abs(t-1)); s = sin(a) + sin(b); % Simpson's rule for i = 1:n-1 s = s + 2*sin(a+i*h); end for i = 1:n s = s + ?*sin( ? ); end s = s*h/6; es(n) = log(abs(s-1)); g2 = 0; % Two-point Gaussian quadrature for i = 1:n g2 = g2 + sin(a+(i-0.5)*h-0.5*h/sqrt(3)) ... + sin(a+(i-0.5)*h+0.5*h/sqrt(3)); end g2 = 0.5*h*g2; eg2(n) = log(abs(g2-1)); disp(sprintf('%3.0f %12.8f %12.8f %12.8f %12.8f %12.8f',n,rl,rm,t,s,g2)) end x = []; for n = 1:M x(n) = log(n); end plot(x, erl, 'o', x, erm, '--', x, et, '-', x, es, '-.', x, eg2, ':') print error1.ps
% integral2.m : comparison of the trapezoid rule, the extrapolated % trapezoid rule, and the Simpson rule. a = 0; b = pi/2; M = 8; for n = 1:M+1 k = 2^(n-1); h = (b-a)/k; t = sin(a) + sin(b); % trapezoid rule for i = 1:k-1 t = t + 2*sin(a+i*h); end t = 0.5*h*t; tra(n) = t; s = sin(a) + sin(b); % Simpson's rule for i = 1:k-1 s = s + 2*sin(a+i*h); end for i = 1:k s = s + 4*sin(a+(i-0.5)*h); end s = s*h/6; sim(n) = s; end for n = 1:M % Richardson extrapolation etr(n) = ? ; end for n = 1:M k = 2^(n-1); disp(sprintf('%3.0f %20.16f %20.16f %20.16f',k,tra(n),etr(n),sim(n))) end disp(sprintf('%3.0f %20.16f %41.16f',256,tra(9),sim(9))) disp(sprintf('')) for n = 1:M k = 2^(n-1); eetr(n) = log(abs(etr(n)-1)); etra(n) = log(abs(tra(n)-1)); esim(n) = log(abs(sim(n)-1)); disp(sprintf('%3.0f %20.16f %20.16f %20.16f',k,etra(n),eetr(n),esim(n))) end x = []; for n = 1:M x(n) = log(2^(n-1)); end plot(x, etra, ':', x, eetr, '-', x, esim, 'o') print error2.ps