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
?? ba dxxfI )(
An m-point quadrature rule Q for the definite integral
is an approximation of the form
?
?
??
m
k
kk xfwabQ
1
)()(
abscissasweights
Efficiency essentially depends upon the
number of function evaluations.
The Newton-Cotes Rules
)()( xfxp ?
dxxfdxxp b
a
b
a ??
? )()(
M-point Newton-Cotes rule
dxxpQ b
a mmNC ? ?
? )(1)(
where pm-1(x) interpolates f(x) at
miab
m
iax
i,1),(1
1 ??
?
???
m=2 trapezoidal rule
?
?
?
?
?
?
???
?
?
?
?
?
?
?
?
?
?? ?
)(
2
1
)(
2
1
)(
)(
)()(
)(
)2(
bfafab
dxax
ab
afbf
afQ
b
a
NC
m=3 and c=(a+b)/2 Simpson rule
?
?
?
?
?
?
?
?
??
?
?
?
?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
??
?
?
?? ?
)(
2
4)(
6
))((
)()()()(
)(
)()(
)(
)3(
bf
ba
faf
ab
dxcxax
ab
ac
afcf
cb
cfbf
ax
ac
afcf
afQ
b
a
NC
? ?? ?
?
?
?
? ?
?
??
?
? ?? m
k
k
i
ikm xxcxp
1
1
1
1 )(
? ? ?? ???
?
???
? ??? ?
??
?
b
a
b
a
k
i
i
m
m
kmmNC dxxxcdxxpQ
1
11
1)( )()(
? ?
?
?
? ???
1
0
1
1)( )(
m m
k
mk
k
kmmNC ShcdsshaphQ
? ?? ??
?
?
???
?
???
? ??? 1
0
1
1
1
m k
i
mk dsisS
Equal spacing,x=a+sh
)!3/()23(
)2/()2(
,/)(,
3
12344
2
1233
12211
hffffc
hfffc
hffcfc
????
???
???
?
?
? ?
?
?
? ?
??????
?????
??????
1
0
22
4
1
0
2
3
1
0
1
0
2
11
4/)3()1()2)(1(
,3/)2/5()1()1(
2/)1(),1(1
m
m
m
m
m m
mm
mmdssssS
mmdsssS
ms d sSmdsS
8/)33)((
)33(
8
3
4321
4321)4(
ffffab
ffff
h
Q NC
?????
????
Weight [1 3 3 1]/8 for QNC(4)
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)
?
?
?
?
?
?
?
?
?
?
???? ?
?
)(
)(
],,)[()(
1
1
1
)(
m
m
m
i
iimNC
xf
xf
wwabfwabQ ??
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 )
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
If f(x) and its first four derivatives are
continuours on [a,b],then
4
5
)3( 2880
)()( MabQdxxfb
a NC
????
Where M4 is an upper bound on |f(4)(x)| on [a,b]
Note that if f(x) is a cubic polynomial,then f(4)(x)=0 and so
Simposn's rule is exact,This is somewhat surprising because
The rule is based on the integration of a quadratic interpolant.
2
)1(
)( 1)()(
?
? ?
?
?
?
?
?
?
???
?
d
d
mmNC
b
a m
abfcQdxxf ?
where cm is a small constant,and
?
?
? ?
?
e v e n is if
e v e n is if1
mm
mm
d
2
1)( 1)(
?
? ?
?
?
?
?
?
?
???
?
d
dmmNC
b
a m
abMcQdxxf
1
)1( )( If
?
? ?
d
d Mxf
function error = NCerror(a,b,m,M)
if m==2,d=1; c = -1/12;
elseif m==3,d=3; c = -1/90;
elseif m==4,d=3; c = -3/80;
elseif m==5,d=5; c = -8/945;
elseif m==6,d=5; c = -275/12096;
elseif m==7,d=7; c = -9/1400;
elseif m==8,d=7; c = -8183/518400;
elseif m==9,d=9; c = -2368/467775;
elseif m==10,d=9; c = -173/14620;
else d=11; c = -1346350/326918592;
end;
h = (b-a)/(m-1);
error = abs(c*M*h^(d+2));
m QNC(m) Error Error Bound
2 0.7853981633974483 2.146e-001 3.230e-001
3 1.0022798774922104 2.280e-003 3.321e-003
4 1.0010049233142790 1.005e-003 1.476e-003
5 0.9999915654729927 8.435e-006 1.219e-005
6 0.9999952613861668 4.739e-006 6.867e-006
7 1.0000000258372352 2.584e-008 3.714e-008
8 1.0000000158229039 1.582e-008 2.277e-008
9 0.9999999999408976 5.910e-011 8.466e-011
10 0.9999999999621676 3.783e-011 5.417e-011
11 1.0000000000001021 1.021e-013 1.460e-013
? 2/0 )s in (? dxx
Open Newton-Cotes Rules
Close Newton-Cotes Rules
For m=3,5,6,7,…,the open rules involve negative weights,
a feature that can unermine the numerical stability of the
rule,(Numerical cancellation can occur with positive
integrands,and this can jeopardize the relative error.)
Composite Rules
? ??
?
??
n
i
z
z
b
a
i
i
dxxfdxxf
1
1 )()(
Composite Simpson rule
?
?
?? ??
?? n
i
iii
i zfzfzfQ
1
12/1 ))()(4)((6
2/)(,12/11 ??? ????? iiiiii zzzzz
numI=0
For I=1:length(z)-1
numI=numI+QNC(fname,z(I),z(I+1),m)
end
n
abiaz
i
??????? )1(
numI=0; Delta=(b-a)/n
For I=1:n
numI=numI+QNC('f',a+(I-1)*Delta,a+I*Delta,m)
end
))17:13()13:9()9:5()5:1(()4( )5( fwfwfwfwQ TTTTNC ?????
f(1:n)
)( )(n mNCQ
? ? ? ?
function numI = CompQNC(fname,a,b,m,n)
Delta = (b-a)/n; h = Delta/(m-1);
x = a+h*(0:(n*(m-1)))'; w = NCWeights(m);
x = linspace(a,b,n*(m-1)+1)'; f = feval(fname,x);
numI = 0; first = 1;
last = m;
for i=1:n
numI = numI + w'*f(first:last);
first = last;
last = last+m-1;
end
numI = Delta*numI;
? ??? ??
?
?
??
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
???? ? n
i
d
ii
i
d
mi
n
i
z
z
b
a m
zzfcQdxxfdxxf i
i 1
2
1)1(
1 1
)()( 1 ?
?
?
?
n
i
i
n
mNC QQ
1
)(
)( n
abzz
ii
?????
? 1
?
?
?? ?
n
i
d
i
d ff
n 1
)1()1( )()(1 ??
)(
)1(
)( )1(
2
)(
)( ?
?
?
???
?
???
?
?
???? dd
m
n
mNC
b
a
nf
mn
abcQdxxf
1
2
1
)(
)(
1
)1(
)( ?
?
?
?
?
?
?
?
?
?
?
???
?
???
?
?
???
? d
d
dm
n
mNC
b
a nm
abMcQdxxf
2
1)( 1)(
?
? ?
?
?
?
?
?
?
???
?
d
dmmNC
b
a m
abMcQdxxf
1
2
1
)(
)(
1
)1(
)( ?
?
?
?
?
?
?
?
?
?
?
???
?
???
?
?
???
? d
d
dm
n
mNC
b
a nm
abMcQdxxf
t o l
n
abM ??
?
??
?
? ?
4
5
4
1
290
1
4 4
2 8 8 0
)()(
to l
abMabn
?
???
0 5 10 15 20 25 30 35 40
10-15
10-10
10-5
100
m = 3
m = 5
m = 7
Example 1,QNC(m,n) error for integral of sin from 0 to pi/2
n = number of subintervals.
Er
ro
r in
Q
NC
(m
,n)
Adaptive Quadrature
Adaptive quadrature methods,to discover where the
integrand is ill behaved and shortening the subintervals
accordingly.
If the magnitude of the derivative varies widely across
the interval of integration,then the error control process
may result in an unnecessary number of function
evaluations.
An Adaptive Newton-Cotes Procedure
function numI=AdaptQNC(fname,a,b,…)
Compute the integral from a to b in two ways,Call the values A1 and
A2 and assume that A2 is better.
Estimate the error in A2 base on |A1-A2|.
If the error is sufficiently small,then
numI=A2;
Else
mid=(a+b)/2;
numI=AdaptQNC(fname,a,mid,…)+AdaptQNC(fname,mid,b,….);
end
Details:
Set ）（）（, 2
)(21 )(1 mNCmNC QAQA ??
CAI ?? 1
121
2
2
)1(
2 22
1
1
)( ??
?
? ??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
???
dd
d
d
m
CA
m
abfcAI ?
2
1
)1(
1)(
?
? ?
?
??
?
?
?
?? dd
m m
abfcC ?
12
,
21 1
21
2)1(
12
?
???
?
??
??? dd
AAAIAAC
,
12 1
21
t o ld
AA ??
?
?
?
If A2 is acceptable
function numI = AdaptQNC(fname,a,b,m,tol)
A1 = CompQNC(fname,a,b,m,1);
A2 = CompQNC(fname,a,b,m,2);
d = 2*floor((m-1)/2)+1;
error = (A2-A1)/(2^(d+1)-1);
if abs(error) <= tol
numI = A2+error; % A2 is acceptable
else
mid = (a+b)/2;
numI = AdaptQNC(fname,a,mid,m,tol/2) + …
AdaptQNC(fname,mid,b,m,tol/2);
end
global FunEvals VecFunEvals
close all
x = linspace(0,1,100);y = humps(x);
plot(x,y);u=[];
for tol = [.01,001,0001,00001]
for m=3:2:9
figure; plot(x,y); hold on
FunEvals = 0; VecFunEvals = 0;
num0 = AdaptQNC('SpecHumps',0,1,m,tol);
xLabel(sprintf('Scalar Evals = %3.0f Vector Evals…
= %3.0f',FunEvals,VecFunEvals))
hold off ; u = [u; FunEvals VecFunEvals];
end;end
function y = SpecHumps(x)
global FunEvals VecFunEvals;
y = 1,/ ((x-.3).^2 +,01) + 1,/ ((x-.9).^2 +,04) - 6;
hold on
plot(x,y,'*')
hold off
FunEvals = FunEvals + length(x);
VecFunEvals = VecFunEvals + 1;
global
Quad,Quad8 in Matlab
quad,QNC(3)(Simpson's Rule) quad8,QNC(9)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
10
20
30
40
50
60
70
80
90
100
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
10
20
30
40
50
60
70
80
90
100
quad('specHumps',0,1) quad8('specHumps',0,1)
Special Topic
Gauss quadrature is useful in certain specialized
seetings,When the higher derivatives of a function
are smooth,fewer function evaluations are
required to attain a given level of accuracy.
In setting where the function evaluations are
experimentally determined,spline quadrature has a
certain appeal.
Gauss Quadrature
?? ??1 1 2211 )()()( xfwxfwdxxf
211 121 ?? ??? dxww
01 12211 ??? ?? x d xxwxw
2
31
1
22
22
2
11 ??? ?? dxxxwxw
01 1 3322311 ??? ?? dxxxwxw
,121 ?? ww
3/1,3/1 21 ??? xx
32)(f o r a cc u r a r y dxcxbxaxf ????
two-point Gauss-Legendre rule
?? ???1 1 )3/1()3/1()( ffdxxf
m-point Gauss-Legendre rule
)()( 11)( mmmGL xfwxfwQ ??? ?
where the wi and xi are chosen to make the rule exact for
polynomials of degree 2m-1
12:0,1)1(1
1
2211 ???
?????? ? mk
kxwxwxw
k
k
mm
kk ?
?? ???? 1 111 )(,)()()( kmm xxfdxxfxfwxfw ?
? ????ba dxxgabdxxf 1 1 )(2)(
?
?
??
?
? ???? xabbafxg
22
)(
???????? ?????? ??????????? ????
b
amm
dxxfxabbafwxabbafwab )(
22222 11
?
Integrals from,[a,b]
function numI = QGL(fname,a,b,m)
[w,x] = GLWeights(m);
fvals = feval(fname,((b-a)/2)*x + ((a+b)/2)*ones(m,1));
numI = ((b-a)/2)*w'*fvals;
Approximating the integral of sin from 0 to pi/2.
m NC(m) GL(m)
------------------------------------------------
2 0.7853981633974483 0.9984726134041149
3 1.0022798774922104 1.0000081215555008
4 1.0010049233142790 0.9999999771971151
5 0.9999915654729927 1.0000000000395670
6 0.9999952613861668 0.9999999999999533
function numI = QNC(fname,a,b,m)
w = NCweights(m);
x = linspace(a,b,m)';
f = feval(fname,x);
numI = (b-a)*(w'*f);
function numI = QGL(fname,a,b,m)
[w,x] = GLWeights(m);
fvals = feval(fname,((b-a)/2)*x + ((a+b)/2)*ones(m,1));
numI = ((b-a)/2)*w'*fvals;
Spline Quadrature
?? nxx dxxsI
1
)(
31,22,3,4,)()()()( iiiiiiii xxxxxxxq ??????? ????
41,32,23,
4,432)(
1
i
i
i
i
i
i
ii
x
x i
hhhhdxxqi
i
???? ????? ?
Suppose S(x) is a cubic spline interpolant of (xi,yi)
If the ith cubic is represented by
function numI = SplineQ(x,y)
S = spline(x,y);
[x,rho,L,k] = unmkpp(S);
sum = 0;
for i=1:L
h = x(i+1)-x(i);
subI = h*(((rho(i,1)*h/4 + rho(i,2)/3)*h + rho(i,3)/2)*h + rho(i,4));
sum = sum + subI;
end
numI = sum;
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
?? ba dxxfI )(
An m-point quadrature rule Q for the definite integral
is an approximation of the form
?
?
??
m
k
kk xfwabQ
1
)()(
abscissasweights
Efficiency essentially depends upon the
number of function evaluations.
The Newton-Cotes Rules
)()( xfxp ?
dxxfdxxp b
a
b
a ??
? )()(
M-point Newton-Cotes rule
dxxpQ b
a mmNC ? ?
? )(1)(
where pm-1(x) interpolates f(x) at
miab
m
iax
i,1),(1
1 ??
?
???
m=2 trapezoidal rule
?
?
?
?
?
?
???
?
?
?
?
?
?
?
?
?
?? ?
)(
2
1
)(
2
1
)(
)(
)()(
)(
)2(
bfafab
dxax
ab
afbf
afQ
b
a
NC
m=3 and c=(a+b)/2 Simpson rule
?
?
?
?
?
?
?
?
??
?
?
?
?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
??
?
?
?? ?
)(
2
4)(
6
))((
)()()()(
)(
)()(
)(
)3(
bf
ba
faf
ab
dxcxax
ab
ac
afcf
cb
cfbf
ax
ac
afcf
afQ
b
a
NC
? ?? ?
?
?
?
? ?
?
??
?
? ?? m
k
k
i
ikm xxcxp
1
1
1
1 )(
? ? ?? ???
?
???
? ??? ?
??
?
b
a
b
a
k
i
i
m
m
kmmNC dxxxcdxxpQ
1
11
1)( )()(
? ?
?
?
? ???
1
0
1
1)( )(
m m
k
mk
k
kmmNC ShcdsshaphQ
? ?? ??
?
?
???
?
???
? ??? 1
0
1
1
1
m k
i
mk dsisS
Equal spacing,x=a+sh
)!3/()23(
)2/()2(
,/)(,
3
12344
2
1233
12211
hffffc
hfffc
hffcfc
????
???
???
?
?
? ?
?
?
? ?
??????
?????
??????
1
0
22
4
1
0
2
3
1
0
1
0
2
11
4/)3()1()2)(1(
,3/)2/5()1()1(
2/)1(),1(1
m
m
m
m
m m
mm
mmdssssS
mmdsssS
ms d sSmdsS
8/)33)((
)33(
8
3
4321
4321)4(
ffffab
ffff
h
Q NC
?????
????
Weight [1 3 3 1]/8 for QNC(4)
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)
?
?
?
?
?
?
?
?
?
?
???? ?
?
)(
)(
],,)[()(
1
1
1
)(
m
m
m
i
iimNC
xf
xf
wwabfwabQ ??
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 )
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
If f(x) and its first four derivatives are
continuours on [a,b],then
4
5
)3( 2880
)()( MabQdxxfb
a NC
????
Where M4 is an upper bound on |f(4)(x)| on [a,b]
Note that if f(x) is a cubic polynomial,then f(4)(x)=0 and so
Simposn's rule is exact,This is somewhat surprising because
The rule is based on the integration of a quadratic interpolant.
2
)1(
)( 1)()(
?
? ?
?
?
?
?
?
?
???
?
d
d
mmNC
b
a m
abfcQdxxf ?
where cm is a small constant,and
?
?
? ?
?
e v e n is if
e v e n is if1
mm
mm
d
2
1)( 1)(
?
? ?
?
?
?
?
?
?
???
?
d
dmmNC
b
a m
abMcQdxxf
1
)1( )( If
?
? ?
d
d Mxf
function error = NCerror(a,b,m,M)
if m==2,d=1; c = -1/12;
elseif m==3,d=3; c = -1/90;
elseif m==4,d=3; c = -3/80;
elseif m==5,d=5; c = -8/945;
elseif m==6,d=5; c = -275/12096;
elseif m==7,d=7; c = -9/1400;
elseif m==8,d=7; c = -8183/518400;
elseif m==9,d=9; c = -2368/467775;
elseif m==10,d=9; c = -173/14620;
else d=11; c = -1346350/326918592;
end;
h = (b-a)/(m-1);
error = abs(c*M*h^(d+2));
m QNC(m) Error Error Bound
2 0.7853981633974483 2.146e-001 3.230e-001
3 1.0022798774922104 2.280e-003 3.321e-003
4 1.0010049233142790 1.005e-003 1.476e-003
5 0.9999915654729927 8.435e-006 1.219e-005
6 0.9999952613861668 4.739e-006 6.867e-006
7 1.0000000258372352 2.584e-008 3.714e-008
8 1.0000000158229039 1.582e-008 2.277e-008
9 0.9999999999408976 5.910e-011 8.466e-011
10 0.9999999999621676 3.783e-011 5.417e-011
11 1.0000000000001021 1.021e-013 1.460e-013
? 2/0 )s in (? dxx
Open Newton-Cotes Rules
Close Newton-Cotes Rules
For m=3,5,6,7,…,the open rules involve negative weights,
a feature that can unermine the numerical stability of the
rule,(Numerical cancellation can occur with positive
integrands,and this can jeopardize the relative error.)
Composite Rules
? ??
?
??
n
i
z
z
b
a
i
i
dxxfdxxf
1
1 )()(
Composite Simpson rule
?
?
?? ??
?? n
i
iii
i zfzfzfQ
1
12/1 ))()(4)((6
2/)(,12/11 ??? ????? iiiiii zzzzz
numI=0
For I=1:length(z)-1
numI=numI+QNC(fname,z(I),z(I+1),m)
end
n
abiaz
i
??????? )1(
numI=0; Delta=(b-a)/n
For I=1:n
numI=numI+QNC('f',a+(I-1)*Delta,a+I*Delta,m)
end
))17:13()13:9()9:5()5:1(()4( )5( fwfwfwfwQ TTTTNC ?????
f(1:n)
)( )(n mNCQ
? ? ? ?
function numI = CompQNC(fname,a,b,m,n)
Delta = (b-a)/n; h = Delta/(m-1);
x = a+h*(0:(n*(m-1)))'; w = NCWeights(m);
x = linspace(a,b,n*(m-1)+1)'; f = feval(fname,x);
numI = 0; first = 1;
last = m;
for i=1:n
numI = numI + w'*f(first:last);
first = last;
last = last+m-1;
end
numI = Delta*numI;
? ??? ??
?
?
??
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
???? ? n
i
d
ii
i
d
mi
n
i
z
z
b
a m
zzfcQdxxfdxxf i
i 1
2
1)1(
1 1
)()( 1 ?
?
?
?
n
i
i
n
mNC QQ
1
)(
)( n
abzz
ii
?????
? 1
?
?
?? ?
n
i
d
i
d ff
n 1
)1()1( )()(1 ??
)(
)1(
)( )1(
2
)(
)( ?
?
?
???
?
???
?
?
???? dd
m
n
mNC
b
a
nf
mn
abcQdxxf
1
2
1
)(
)(
1
)1(
)( ?
?
?
?
?
?
?
?
?
?
?
???
?
???
?
?
???
? d
d
dm
n
mNC
b
a nm
abMcQdxxf
2
1)( 1)(
?
? ?
?
?
?
?
?
?
???
?
d
dmmNC
b
a m
abMcQdxxf
1
2
1
)(
)(
1
)1(
)( ?
?
?
?
?
?
?
?
?
?
?
???
?
???
?
?
???
? d
d
dm
n
mNC
b
a nm
abMcQdxxf
t o l
n
abM ??
?
??
?
? ?
4
5
4
1
290
1
4 4
2 8 8 0
)()(
to l
abMabn
?
???
0 5 10 15 20 25 30 35 40
10-15
10-10
10-5
100
m = 3
m = 5
m = 7
Example 1,QNC(m,n) error for integral of sin from 0 to pi/2
n = number of subintervals.
Er
ro
r in
Q
NC
(m
,n)
Adaptive Quadrature
Adaptive quadrature methods,to discover where the
integrand is ill behaved and shortening the subintervals
accordingly.
If the magnitude of the derivative varies widely across
the interval of integration,then the error control process
may result in an unnecessary number of function
evaluations.
An Adaptive Newton-Cotes Procedure
function numI=AdaptQNC(fname,a,b,…)
Compute the integral from a to b in two ways,Call the values A1 and
A2 and assume that A2 is better.
Estimate the error in A2 base on |A1-A2|.
If the error is sufficiently small,then
numI=A2;
Else
mid=(a+b)/2;
numI=AdaptQNC(fname,a,mid,…)+AdaptQNC(fname,mid,b,….);
end
Details:
Set ）（）（, 2
)(21 )(1 mNCmNC QAQA ??
CAI ?? 1
121
2
2
)1(
2 22
1
1
)( ??
?
? ??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
???
dd
d
d
m
CA
m
abfcAI ?
2
1
)1(
1)(
?
? ?
?
??
?
?
?
?? dd
m m
abfcC ?
12
,
21 1
21
2)1(
12
?
???
?
??
??? dd
AAAIAAC
,
12 1
21
t o ld
AA ??
?
?
?
If A2 is acceptable
function numI = AdaptQNC(fname,a,b,m,tol)
A1 = CompQNC(fname,a,b,m,1);
A2 = CompQNC(fname,a,b,m,2);
d = 2*floor((m-1)/2)+1;
error = (A2-A1)/(2^(d+1)-1);
if abs(error) <= tol
numI = A2+error; % A2 is acceptable
else
mid = (a+b)/2;
numI = AdaptQNC(fname,a,mid,m,tol/2) + …
AdaptQNC(fname,mid,b,m,tol/2);
end
global FunEvals VecFunEvals
close all
x = linspace(0,1,100);y = humps(x);
plot(x,y);u=[];
for tol = [.01,001,0001,00001]
for m=3:2:9
figure; plot(x,y); hold on
FunEvals = 0; VecFunEvals = 0;
num0 = AdaptQNC('SpecHumps',0,1,m,tol);
xLabel(sprintf('Scalar Evals = %3.0f Vector Evals…
= %3.0f',FunEvals,VecFunEvals))
hold off ; u = [u; FunEvals VecFunEvals];
end;end
function y = SpecHumps(x)
global FunEvals VecFunEvals;
y = 1,/ ((x-.3).^2 +,01) + 1,/ ((x-.9).^2 +,04) - 6;
hold on
plot(x,y,'*')
hold off
FunEvals = FunEvals + length(x);
VecFunEvals = VecFunEvals + 1;
global
Quad,Quad8 in Matlab
quad,QNC(3)(Simpson's Rule) quad8,QNC(9)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
10
20
30
40
50
60
70
80
90
100
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
10
20
30
40
50
60
70
80
90
100
quad('specHumps',0,1) quad8('specHumps',0,1)
Special Topic
Gauss quadrature is useful in certain specialized
seetings,When the higher derivatives of a function
are smooth,fewer function evaluations are
required to attain a given level of accuracy.
In setting where the function evaluations are
experimentally determined,spline quadrature has a
certain appeal.
Gauss Quadrature
?? ??1 1 2211 )()()( xfwxfwdxxf
211 121 ?? ??? dxww
01 12211 ??? ?? x d xxwxw
2
31
1
22
22
2
11 ??? ?? dxxxwxw
01 1 3322311 ??? ?? dxxxwxw
,121 ?? ww
3/1,3/1 21 ??? xx
32)(f o r a cc u r a r y dxcxbxaxf ????
two-point Gauss-Legendre rule
?? ???1 1 )3/1()3/1()( ffdxxf
m-point Gauss-Legendre rule
)()( 11)( mmmGL xfwxfwQ ??? ?
where the wi and xi are chosen to make the rule exact for
polynomials of degree 2m-1
12:0,1)1(1
1
2211 ???
?????? ? mk
kxwxwxw
k
k
mm
kk ?
?? ???? 1 111 )(,)()()( kmm xxfdxxfxfwxfw ?
? ????ba dxxgabdxxf 1 1 )(2)(
?
?
??
?
? ???? xabbafxg
22
)(
???????? ?????? ??????????? ????
b
amm
dxxfxabbafwxabbafwab )(
22222 11
?
Integrals from,[a,b]
function numI = QGL(fname,a,b,m)
[w,x] = GLWeights(m);
fvals = feval(fname,((b-a)/2)*x + ((a+b)/2)*ones(m,1));
numI = ((b-a)/2)*w'*fvals;
Approximating the integral of sin from 0 to pi/2.
m NC(m) GL(m)
------------------------------------------------
2 0.7853981633974483 0.9984726134041149
3 1.0022798774922104 1.0000081215555008
4 1.0010049233142790 0.9999999771971151
5 0.9999915654729927 1.0000000000395670
6 0.9999952613861668 0.9999999999999533
function numI = QNC(fname,a,b,m)
w = NCweights(m);
x = linspace(a,b,m)';
f = feval(fname,x);
numI = (b-a)*(w'*f);
function numI = QGL(fname,a,b,m)
[w,x] = GLWeights(m);
fvals = feval(fname,((b-a)/2)*x + ((a+b)/2)*ones(m,1));
numI = ((b-a)/2)*w'*fvals;
Spline Quadrature
?? nxx dxxsI
1
)(
31,22,3,4,)()()()( iiiiiiii xxxxxxxq ??????? ????
41,32,23,
4,432)(
1
i
i
i
i
i
i
ii
x
x i
hhhhdxxqi
i
???? ????? ?
Suppose S(x) is a cubic spline interpolant of (xi,yi)
If the ith cubic is represented by
function numI = SplineQ(x,y)
S = spline(x,y);
[x,rho,L,k] = unmkpp(S);
sum = 0;
for i=1:L
h = x(i+1)-x(i);
subI = h*(((rho(i,1)*h/4 + rho(i,2)/3)*h + rho(i,3)/2)*h + rho(i,4));
sum = sum + subI;
end
numI = sum;
