1 / 12

Lecture 26: Numerical Integration

Lecture 26: Numerical Integration. Newton-Cotes Formulas. Trapezoid rule. Simpson's rule. Simpson's 3/8 rule. Boole’s rule. trapezoid1.m. function f1=trapezoid1(func1,a,b,n) %Integration function func1 %using composite trapezoid rule at n points between a and b

kalli
Download Presentation

Lecture 26: Numerical Integration

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 26: Numerical Integration Newton-Cotes Formulas Trapezoid rule Simpson's rule Simpson's 3/8 rule Boole’s rule

  2. trapezoid1.m function f1=trapezoid1(func1,a,b,n) %Integration function func1 %using composite trapezoid rule at n points between a and b h=(b-a)/n; f1=func1(a)/2; for i=1:n-1 f1=f1+func1(a+i*h); end f1=f1+func1(b)/2; f1=h*f1; trapezoid1test.m

  3. %examples of use of function f1=trapezoid1(func1,a,b,n) %Integration function func1 %using composite trapezoid rule at n points between a and b a=1; b=5; f1=inline('exp(-x)') disp(['a=',num2str(a),' b=',num2str(b)]) %display limits of integration intexact=-exp(-b)+exp(-a); i=1:10; nvec=2.^(i); for i=1:length(nvec) f1int=trapezoid1(f1,a,b,nvec(i)); disp(['n=',num2str(nvec(i)),' Error in composite trapezod rule =',num2str(f1int-intexact)]); %display value of function f(x) end

  4. >> trapezoid1test f1 = Inline function: f1(x) = exp(-x) a=1 b=5 n=2 Error in composite trapezod rule =0.11305 n=4 Error in composite trapezod rule =0.029605 n=8 Error in composite trapezod rule =0.0074926 n=16 Error in composite trapezod rule =0.001879 n=32 Error in composite trapezod rule =0.00047011 n=64 Error in composite trapezod rule =0.00011755 n=128 Error in composite trapezod rule =2.9389e-005 n=256 Error in composite trapezod rule =7.3474e-006 n=512 Error in composite trapezod rule =1.8369e-006 n=1024 Error in composite trapezod rule =4.5922e-007

  5. a=0; b=2*pi; f1=inline('exp(sin(x)^2)') disp(['a=',num2str(a),' b=',num2str(b)]) %display limits of integration intexact=trapezoid1(f1,a,b,200); i=1:5; nvec=2.^(i); for i=1:length(nvec) f1int=trapezoid1(f1,a,b,nvec(i)); disp(['n=',num2str(nvec(i)),' Error in composite trapezod rule =',num2str(f1int-intexact)]); %display value of function f(x) end

  6. f1 = Inline function: f1(x) = exp(sin(x)^2) a=0 b=6.2832 n=2 Error in composite trapezod rule =-4.7337 n=4 Error in composite trapezod rule =0.66447 n=8 Error in composite trapezod rule =0.0034145 n=16 Error in composite trapezod rule =7.8954e-009 n=32 Error in composite trapezod rule =-3.5527e-015 f1 = Inline function: f1(x) = exp(sin(x)^2) a=0 b=0.7854 n=2 Error in composite trapezod rule =0.021343 n=4 Error in composite trapezod rule =0.005305 n=8 Error in composite trapezod rule =0.0013228 n=16 Error in composite trapezod rule =0.00032898 n=32 Error in composite trapezod rule =8.0649e-005

  7. a=0; b=pi/4; f1=inline('exp(sin(x)^2)') disp(['a=',num2str(a),' b=',num2str(b)]) %display limits of integration intexact=trapezoid1(f1,a,b,200); i=1:5; nvec=2.^(i); for i=1:length(nvec) f1int=trapezoid1(f1,a,b,nvec(i)); disp(['n=',num2str(nvec(i)),' Error in composite trapezod rule =',num2str(f1int-intexact)]); %display value of function f(x) end

  8. simpson1.m function f1=simpson1(func1,a,b,n) %Integration function func1 %using composite simpson rule at n points between a and b h=(b-a)/n; f1=func1(a); for i=1:(n/2-1) f1=f1+4*func1(a+(2*i-1)*h)+2*func1(a+(2*i)*h); end f1=f1+4*func1(a+(n-1)*h)+func1(b); f1=(h/3)*f1; simpson1test.m

  9. %examples of use of function f1=simpson1(func1,a,b,n) %Integration function func1 %using composite Simpson's rule at n points between a and b; n must be even a=1; b=5; f1=inline('exp(-x)') disp(['a=',num2str(a),' b=',num2str(b)]) %display limits of integration intexact=-exp(-b)+exp(-a); i=1:10; nvec=2.^(i); for i=1:length(nvec) f1int=simpson1(f1,a,b,nvec(i)); disp(['n=',num2str(nvec(i)),' Error in composite trapezod rule =',num2str(f1int-intexact)]); %display value of function f(x) end …

  10. >> simpson1test f1(x) = exp(-x) a=1 b=5 n=2 Error in composite trapezod rule =0.021369 n=4 Error in composite trapezod rule =0.0017902 n=8 Error in composite trapezod rule =0.00012176 n=16 Error in composite trapezod rule =7.7793e-006 n=32 Error in composite trapezod rule =4.8892e-007 n=64 Error in composite trapezod rule =3.06e-008 n=128 Error in composite trapezod rule =1.9132e-009 n=256 Error in composite trapezod rule =1.1958e-010 n=512 Error in composite trapezod rule =7.4751e-012 n=1024 Error in composite trapezod rule =4.6674e-013 The same function with composite trapezoid rule: a=1 b=5 n=2 Error in composite trapezod rule =0.11305 n=4 Error in composite trapezod rule =0.029605 n=8 Error in composite trapezod rule =0.0074926 n=16 Error in composite trapezod rule =0.001879 n=32 Error in composite trapezod rule =0.00047011 n=64 Error in composite trapezod rule =0.00011755 n=128 Error in composite trapezod rule =2.9389e-005 n=256 Error in composite trapezod rule =7.3474e-006 n=512 Error in composite trapezod rule =1.8369e-006 n=1024 Error in composite trapezod rule =4.5922e-007

  11. Exercise Compute integral of function e^(-x^2) between a=0 and b=10 using composite trapezoid and composite Simpson’s rules for n=2,4,8,…, 1024. Compare with exact answer = sqrt(pi)/2 Use simpson1test.m as example.

  12. >> inclass26 f1 = Inline function: f1(x) = exp(-x^2) a=0 b=10 n=2 Error in composite trapezod rule =0.78044 n=4 Error in composite trapezod rule =-0.046459 n=8 Error in composite trapezod rule =-0.1186 n=16 Error in composite trapezod rule =-0.0010671 n=32 Error in composite trapezod rule =-6.2875e-012 n=64 Error in composite trapezod rule =0 n=128 Error in composite trapezod rule =2.2204e-016 n=256 Error in composite trapezod rule =5.5511e-016 n=512 Error in composite trapezod rule =0 n=1024 Error in composite trapezod rule =-1.3323e-015 The same function with composite trapezoid rule: a=0 b=10 n=2 Error in composite trapezod rule =1.6138 n=4 Error in composite trapezod rule =0.3686 n=8 Error in composite trapezod rule =0.0032014 n=16 Error in composite trapezod rule =1.8863e-011 n=32 Error in composite trapezod rule =0 n=64 Error in composite trapezod rule =-1.1102e-016 n=128 Error in composite trapezod rule =-1.1102e-016 n=256 Error in composite trapezod rule =-4.4409e-016 n=512 Error in composite trapezod rule =-6.6613e-016 n=1024 Error in composite trapezod rule =-4.4409e-016

More Related