1 / 11

Lecture 28: Comparison of different numerical integrators

Lecture 28: Comparison of different numerical integrators. Adaptive Simpson’s and Trapezoid Rules 2. Romberg Integration 3. Adaptive Gaussian Quadrature. adaptivesimpson.m. %test of adaptive Simpson's rule and comparison with composite Trapezoid and Composite Simpson's tol=10^(-10); a=1;

nusa
Download Presentation

Lecture 28: Comparison of different numerical integrators

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Lecture 28: Comparison of different numerical integrators • Adaptive Simpson’s and Trapezoid Rules • 2. Romberg Integration • 3. Adaptive Gaussian Quadrature

  2. adaptivesimpson.m %test of adaptive Simpson's rule and comparison with composite Trapezoid and Composite Simpson's tol=10^(-10); a=1; b=4; f1=inline('sin(x).^2') disp(['a=',num2str(a),' b=',num2str(b)]) %display limits of integration intexact=quad(f1,a,b,10^(-16)); [fadaptive funccount]=quad(f1,a,b,tol); disp(['Integration using adaptive simpson gives error=',num2str(fadaptive-intexact), ' in ',num2str(funccount), ' function evaluations']); disp(' '); disp(['Now compare with n step of composite simpson rule:']); i=1:10; nvec=2.^(i); ….

  3. >> adaptivesimpson f1 = Inline function: f1(x) = sin(x).^2 a=1 b=4 Integration using adaptive simpson gives error=-1.0503e-013 in 249 function evaluations Now compare with n step of composite simpson rule: n=2 Error in composite simpson rule =-0.12324 n=4 Error in composite simpson rule =-0.00075995 n=8 Error in composite simpson rule =-3.7687e-005 n=16 Error in composite simpson rule =-2.2363e-006 n=32 Error in composite simpson rule =-1.3801e-007 n=64 Error in composite simpson rule =-8.5986e-009 n=128 Error in composite simpson rule =-5.3699e-010 n=256 Error in composite simpson rule =-3.3555e-011 n=512 Error in composite simpson rule =-2.0945e-012 n=1024 Error in composite simpson rule =-1.299e-013 The same function with composite trapezoid rule: a=1 b=4 n=2 Error in composite trapezod rule =0.017886 n=4 Error in composite trapezod rule =0.0039016 n=8 Error in composite trapezod rule =0.00094713 n=16 Error in composite trapezod rule =0.0002351 n=32 Error in composite trapezod rule =5.8673e-005 n=64 Error in composite trapezod rule =1.4662e-005 n=128 Error in composite trapezod rule =3.665e-006 n=256 Error in composite trapezod rule =9.1623e-007 n=512 Error in composite trapezod rule =2.2906e-007 n=1024 Error in composite trapezod rule =5.7264e-008 n=2048 Error in composite trapezod rule =1.4316e-008 n=4096 Error in composite trapezod rule =3.579e-009

  4. romberg1.m % Romberg integration with number of function evaluations count funccount % Computes approximation to definite integral % Inputs: Matlab inline function specifying integrand f, % a,b integration interval, n=number of rows % Output: [r(n,n) funccount] function [r1 funccount] =romberg(f,a,b,n) h=(b-a)./(2.^(0:n-1)); r(1,1)=(b-a)*(f(a)+f(b))/2; funccount=2; for j=2:n subtotal = 0; for i=1:2^(j-2) subtotal = subtotal + f(a+(2*i-1)*h(j)); funccount= funccount+1; end r(j,1) = r(j-1,1)/2+h(j)*subtotal; for k=2:j r(j,k) = (4^(k-1)*r(j,k-1)-r(j-1,k-1))/(4^(k-1)-1); end end r1=r(n,n);

  5. rombergtest.m %test of romberg integration vs. adaptive Simpson's rule and comparison with composite Trapezoid and Composite Simpson's tol=10^(-10); %desired tolerance a=1; b=4; f1=inline('sin(x).^2') %function f(x) disp(['a=',num2str(a),' b=',num2str(b)]) %display limits of integration intexact=quadl(f1,a,b,10^(-16)); %calculation of "exact" value of integral [fadaptive funccount]=quad(f1,a,b,tol); disp(['Integration using adaptive simpson gives error=',num2str(fadaptive-intexact), ' in ',num2str(funccount), ' function evaluations']); disp(['Now compare with Romberg integration:']); i=1:8; nvec=1.*(i); for i=1:length(nvec) [intromber funccount]=romberg1(f1,a,b,nvec(i)); disp([' Error in Romberg integration =',num2str(intromber-intexact), ' in ',num2str(funccount), ' function evaluations']);% end …

  6. >> rombergtest f1 = Inline function: f1(x) = sin(x).^2 a=1 b=4 Integration using adaptive simpson gives error=-1.0525e-013 in 249 function evaluations Now compare with Romberg integration: Error in Romberg integration =0.44125 in 2 function evaluations Error in Romberg integration =-0.12324 in 3 function evaluations Error in Romberg integration =0.0074051 in 5 function evaluations Error in Romberg integration =-0.00010691 in 9 function evaluations Error in Romberg integration =3.8208e-007 in 17 function evaluations Error in Romberg integration =-3.4053e-010 in 33 function evaluations Error in Romberg integration =7.5717e-014 in 65 function evaluations Error in Romberg integration =2.2204e-016 in 129 function evaluations Now compare with n step of composite simpson rule: n=2 Error in composite simpson rule =-0.12324 n=4 Error in composite simpson rule =-0.00075995 n=8 Error in composite simpson rule =-3.7687e-005 n=16 Error in composite simpson rule =-2.2363e-006 n=32 Error in composite simpson rule =-1.3801e-007 n=64 Error in composite simpson rule =-8.5986e-009 n=128 Error in composite simpson rule =-5.3699e-010 n=256 Error in composite simpson rule =-3.3555e-011 n=512 Error in composite simpson rule =-2.0948e-012 n=1024 Error in composite simpson rule =-1.3012e-013 …

  7. rombergtest2.m %rombergtest2.m: test of romberg integration vs. adaptive Simpson's rule and comparison with composite Trapezoid and Composite Simpson's tol=10^(-10); %desired tolerance a=1; b=15; f1=inline('sin(x).^2') %function f(x) disp(['a=',num2str(a),' b=',num2str(b)]) %display limits of integration intexact=quadl(f1,a,b,10^(-16)); %calculation of "exact" value of integral [fadaptive funccount]=quad(f1,a,b,tol); disp(['Integration using adaptive simpson gives error=',num2str(fadaptive-intexact), ' in ',num2str(funccount), ' function evaluations']); disp(' '); disp(['Now compare with Romberg integration:']); disp(' '); …

  8. >> rombergtest2 f1 = Inline function: f1(x) = sin(x).^2 a=1 b=15 Integration using adaptive simpson gives error=-1.4833e-013 in 1233 function evaluations Now compare with Romberg integration: Error in Romberg integration =0.4423 in 2 function evaluations Error in Romberg integration =4.3003 in 3 function evaluations Error in Romberg integration =4.1559 in 5 function evaluations Error in Romberg integration =-2.6801 in 9 function evaluations Error in Romberg integration =0.23911 in 17 function evaluations Error in Romberg integration =-0.004795 in 33 function evaluations Error in Romberg integration =2.3442e-005 in 65 function evaluations Error in Romberg integration =-2.8473e-008 in 129 function evaluations Error in Romberg integration =8.6331e-012 in 257 function evaluations Error in Romberg integration =0 in 513 function evaluations Error in Romberg integration =-8.8818e-016 in 1025 function evaluations Error in Romberg integration =-7.9936e-015 in 2049 function evaluations Now compare with n step of composite simpson rule: n=2 Error in composite simpson rule =4.3003 n=4 Error in composite simpson rule =4.165 n=8 Error in composite simpson rule =-2.1522 n=16 Error in composite simpson rule =0.037939 n=32 Error in composite simpson rule =0.0016978 n=64 Error in composite simpson rule =9.8788e-005 n=128 Error in composite simpson rule =6.0685e-006 n=256 Error in composite simpson rule =3.7766e-007 n=512 Error in composite simpson rule =2.3579e-008 n=1024 Error in composite simpson rule =1.4733e-009

  9. gaussianadaptive.m %gaussianadaptive.m: test of gaussian adaptive integration vs. adaptive Simpson's rule and Romberg's integration tol=10^(-8); %desired tolerance a=1; b=4; f1=inline('sin(x).^2') %function f(x) disp(['a=',num2str(a),' b=',num2str(b)]) %display limits of integration intexact=quadl(f1,a,b,10^(-16)); %calculation of "exact" value of integral [fadaptive funccount]=quadl(f1,a,b,tol); disp(['Integration using adaptive Gaussian quadrature gives error=',num2str(fadaptive-intexact), ' in ',num2str(funccount), ' function evaluations']); disp(' '); [fadaptive funccount]=quad(f1,a,b,10^(-10)); disp(['Integration using adaptive simpson gives error=',num2str(fadaptive-intexact), ' in ',num2str(funccount), ' function evaluations']); disp(' '); disp(['Now compare with Romberg integration:']); …

  10. >> gaussianadaptive f1 = Inline function: f1(x) = sin(x).^2 a=1 b=4 Integration using adaptive Gaussian quadrature gives error=2.1849e-013 in 48 function evaluations Integration using adaptive simpson gives error=-1.0525e-013 in 249 function evaluations Now compare with Romberg integration: Error in Romberg integration =0.44125 in 2 function evaluations Error in Romberg integration =-0.12324 in 3 function evaluations Error in Romberg integration =0.0074051 in 5 function evaluations Error in Romberg integration =-0.00010691 in 9 function evaluations Error in Romberg integration =3.8208e-007 in 17 function evaluations Error in Romberg integration =-3.4053e-010 in 33 function evaluations Error in Romberg integration =7.5717e-014 in 65 function evaluations Error in Romberg integration =2.2204e-016 in 129 function evaluations Error in Romberg integration =-2.2204e-016 in 257 function evaluations Error in Romberg integration =-4.4409e-016 in 513 function evaluations Error in Romberg integration =4.4409e-016 in 1025 function evaluations Error in Romberg integration =0 in 2049 function evaluations

  11. >> adaptivesimpson f1 = Inline function: f1(x) = sin(x).^2 a=1 b=4 Integration using adaptive simpson gives error=-1.0503e-013 in 249 function evaluations Now compare with n step of composite simpson rule: n=2 Error in composite simpson rule =-0.12324 n=4 Error in composite simpson rule =-0.00075995 n=8 Error in composite simpson rule =-3.7687e-005 n=16 Error in composite simpson rule =-2.2363e-006 n=32 Error in composite simpson rule =-1.3801e-007 n=64 Error in composite simpson rule =-8.5986e-009 n=128 Error in composite simpson rule =-5.3699e-010 n=256 Error in composite simpson rule =-3.3555e-011 n=512 Error in composite simpson rule =-2.0945e-012 n=1024 Error in composite simpson rule =-1.299e-013 The same function with composite trapezoid rule: a=1 b=4 n=2 Error in composite trapezod rule =0.017886 n=4 Error in composite trapezod rule =0.0039016 n=8 Error in composite trapezod rule =0.00094713 n=16 Error in composite trapezod rule =0.0002351 n=32 Error in composite trapezod rule =5.8673e-005 n=64 Error in composite trapezod rule =1.4662e-005 n=128 Error in composite trapezod rule =3.665e-006 n=256 Error in composite trapezod rule =9.1623e-007 n=512 Error in composite trapezod rule =2.2906e-007 n=1024 Error in composite trapezod rule =5.7264e-008 n=2048 Error in composite trapezod rule =1.4316e-008 n=4096 Error in composite trapezod rule =3.579e-009

More Related