Introduction to scientific Computing A Matrix Vector Approach Using Matlab Written by Charles FVan Loan 陈文斌 Wbchen(fudan. edu. cn 复日大学
Introduction to Scientific Computing -- A Matrix Vector Approach Using Matlab Written by Charles F.Van Loan 陈 文 斌 Wbchen@fudan.edu.cn 复旦大学
Numerical Integration The Newton -Cotes rules Composite rules Adaptive Quadrature Special Topics Shared Memory Adaptive Quadrature
Numerical Integration • The Newton-Cotes Rules • Composite Rules • Adaptive Quadrature • Special Topics • Shared Memory Adaptive Quadrature
An m-point quadrature rule q for the definite integral b I=I f(x)dx is an approximation of the form Q=(b-a∑wf(x) k=1 weights abscissas Efficiency essentially depends upon the number offunction evaluations
= b a I f (x)dx An m-point quadrature rule Q for the definite integral is an approximation of the form = = − m k k k Q b a w f x 1 ( ) ( ) abscissas weights Efficiency essentially depends upon the number of function evaluations
The Newton-Cotes rules p(x)≈f(x) b p(x)x≈f( x)ax C M-point Newton-Cotes rule NC(m QPm(rdx where pm-l)interpolates f(x)at X.=a+ (b-a)
The Newton-Cotes Rules p(x) f (x) p x dx f x dx b a b a ( ) ( ) M-point Newton-Cotes rule Q p x dx b a NC m = m− ( ) ( ) 1 where pm-1 (x) interpolates f(x) at b a i m m i xi a ( ), 1: 1 1 − = − − = +
m=2 trapezoidal rule QC(2)=f(a)+ f(b)-∫(a) x-a b (b-af(a)+f(b) m=3 and c=(a+b/2 Simpsonrule f(b) -f(c f(c)-f(a) C x-a)+ 6-cc-a NC(3) x-ax-c c-a b-a a)+4/a+b f(b 2
m=2 trapezoidal rule = − + − − − = + ( ) 2 1 ( ) 2 1 ( ) ( ) ( ) ( ) ( ) (2) b a f a f b x a dx b a f b f a Q f a b a NC m=3 and c=(a+b)/2 Simpson rule + + + − = − − − − − − − − − + − − = + ( ) 2 ( ) 4 6 ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (3) f b a b f a f b a x a x c dx b a c a f c f a b c f b f c x a c a f c f a Q f a b a NC
m-1 ()=∑c∑ i=1 k Q Nc (m) pm(x)dk=∑∫I(x-x) Equal spacing, x=a+sh ONc(m)=h Pm-1(a+sh)ds=2ckh'Smk mk s-1+ C1=f1,C2=(f2-f)/h C3=(f3-2f2+f1)/(2h2) C4=(f4-3/3+2/2-f1)(3h)
( ) = − = − = − m k k i m k i p x c x x 1 1 1 1 ( ) = = − − = = − b a b a k i i m m QNC m pm x dx ck x x dx 1 1 1 ( ) 1 ( ) ( ) − = = − + = 1 0 1 ( ) 1 ( ) m m k m k k QNC m h pm a sh ds ck h S ( ) − − = = − + 1 0 1 1 1 m k i Smk s i ds Equal spacing, x=a+sh ( 3 2 )/(3! ) ( 2 )/(2 ) , ( )/ , 3 4 4 3 2 1 2 3 3 2 1 1 1 2 2 1 c f f f f h c f f f h c f c f f h = − + − = − + = = −
"1ds=(m-1) ml S=(m-1)2/2 0 3÷/m-1 S(s-1)d=(m-1)2(m-5/2)/3, 0 m4 S(s-1)(s-2)ds=(m-1)2(m-3)2/4 QC(4)=。(+32+3/+f4) =(b-a)(f1+3f2+3f3+f4)/8 Weight [133 1]8 for ONC4
− − − − = − − = − − = − = − − = = − = = − 1 0 2 2 4 1 0 2 3 1 0 1 0 2 1 1 ( 1)( 2) ( 1) ( 3) / 4 ( 1) ( 1) ( 5/ 2)/ 3, 1 ( 1), ( 1) / 2 m m m m m m m m S s s s ds m m S s s ds m m S ds m S sds m ( )( 3 3 )/8 ( 3 3 ) 8 3 1 2 3 4 (4) 1 2 3 4 b a f f f f f f f f h QNC = − + + + = + + + Weight [1 3 3 1]/8 for QNC(4)
function w= NCweights( m) ifm==2,w=[l1]/2; elseif n==3,w[141]/6; w(1:m)=w(m:-1:1) elseif n==4,w=[1331/8 elseif n==5,w=[73212327y/90 elseif n==6,w=[197550507519/288; elseif m--7,w=[41216272722721641/840; elseif m-8,w[751357713232989298913233577751/17280; elseif m=9,w=989588892810496-4540104969285888 989]/28350; elseif m==10,w=[28571574110801934457785778193441080 157412857]/89600; else w=[16067106300-48525272400-260550427368 2605502724004852510630016067/59872 end
function w = NCweights(m) if m==2, w=[1 1]'/2; elseif m==3, w=[1 4 1]'/6; elseif m==4, w=[1 3 3 1]'/8; elseif m==5, w=[7 32 12 32 7]'/90; elseif m==6, w=[19 75 50 50 75 19]'/288; elseif m==7, w=[41 216 27 272 27 216 41]'/840; elseif m==8, w=[751 3577 1323 2989 2989 1323 3577 751]'/17280; elseif m==9, w=[989 5888 -928 10496 -4540 10496 -928 5888 … 989]'/28350; elseif m==10, w=[2857 15741 1080 19344 5778 5778 19344 1080… 15741 2857]'/89600; else w=[16067 106300 -48525 272400 -260550 427368 … -260550 272400 -48525 106300 16067]'/598752; end w(1:m)=w(m:-1:1)
f(x1) Qw(m=(b-a∑w=(b-a)w,…,wm f(m) function numI=QNC(fname, a, b, m) w=NCweights(m) x= linspace(a, b,m) f= feval(fname, x) numl=(b-a) *(w f
= − = − = ( ) ( ) ( ) ( )[ , , ] 1 1 1 ( ) m m m i NC m i i f x f x Q b a w f b a w w function numI = QNC(fname,a,b,m) w = NCweights(m); x = linspace(a,b,m)'; f = feval(fname,x); numI = (b-a)*(w'*f);
m QNC('sin, 0, pi/2, m) 207853981633974483 31.0022798774922104 1.0010049233142790 50.9999915654729927 60.9999952613861668 7100000058372352 81.0000015829039 90.999999408976 100.9999999621676 1110000021
m QNC('sin', 0, pi/2,m ) 2 0.7853981633974483 3 1.0022798774922104 4 1.0010049233142790 5 0.9999915654729927 6 0.9999952613861668 7 1.0000000258372352 8 1.0000000158229039 9 0.9999999999408976 10 0.9999999999621676 11 1.0000000000001021