自适应梯形积分算法
引言:笔者为学习数值分析的苦逼大二学生,为了做展示需要用到自适应梯形积分算法,但是遍历我的书本和百度也没能找到,故自己根据自适应辛普森公式推导之。
一、梯形公式和辛普森公式
对于数值积分,我们最容易想到的就是梯形公式了,也就是取该段函数的起点终点的平均值作为整段函数的估计值。就像下图所示:
![image-20201204224416709](https://i-blog.csdnimg.cn/blog_migrate/2e0f503f5bb37ab1d1659bdab100cf00.png)
书本上的公式肯定是对的(\狗头)
二、复化梯形公式
有时候,为了提高精度,而不愿意改变算法,我们可以采用复化的形式,也就是把函数分成多段,每一段用梯形公式进行计算得到该段的估计值,再进行累加。(下图为书本上的公式)
三、自适应复合梯形公式推导
但是,我们往往碰到的函数可能一段上很均匀,但是在另一段上波动很大,例如我们这里测试用的函数(下图为在定义域为
(
1
,
3
)
(1,3)
(1,3)上的形状)。
F
=
(
100
/
x
2
)
∗
(
s
i
n
(
10
/
x
)
)
;
F = (100/x^2)*(sin(10/x));
F=(100/x2)∗(sin(10/x));
这时候我们的自适应梯形公式就可以派上用场了,它的实现手段就是在波动较大的地方多进行分段,而波动较小的地方则较少的分段足以达到精度要求。而如何检验波动大小就成了问题的痛点所在,下面我们重点放在如何分析现有估计的误差,并作出合理的判断是否继续细分。
1.推导
先利用梯形公式:
∫
a
b
f
(
x
)
d
x
=
S
(
a
,
b
)
−
h
3
12
f
′
′
(
ξ
)
(1)
\int_a^bf(x)dx = S(a,b)-\frac{h^3}{12}f''(\xi)\\ \tag{1}
∫abf(x)dx=S(a,b)−12h3f′′(ξ)(1)
其中,
S
(
a
,
b
)
=
h
2
[
f
(
a
)
+
f
(
b
)
]
(2)
\ S(a,b) = \frac{h}{2} [f(a)+f(b)]\\ \tag{2}
S(a,b)=2h[f(a)+f(b)](2)
T h e n , a p p l y t h e C o m p o s i t e T r a p e z o i d a l r u l e w i t h n = 3 a n d s t e p s i z e ( b − a ) = h Then\ ,apply \ the \ Composite \ Trapezoidal\ rule \ with\ n=3\ and\ step\ size\ (b-a)=h\\ Then ,apply the Composite Trapezoidal rule with n=3 and step size (b−a)=h
∫ a b f ( x ) d x = h 2 [ f ( a ) + 2 × f ( a + h 2 ) + f ( b ) ] − ( h 2 ) 2 b − a 12 f ′ ′ ( ξ ˉ ) = S ( a , a + b 2 ) + S ( a + b 2 , b ) − h 3 48 f ′ ′ ( ξ ˉ ) (3) \int_a^bf(x)dx = \frac{h}{2} [f(a)+2\times f(a+\frac{h}{2})+f(b)]-(\frac{h}{2})^2\frac{b-a}{12}f''(\bar \xi)\\ \tag{3} =S(a,\frac{a+b}{2})+S(\frac{a+b}{2},b)-\frac{h^3}{48}f''(\bar \xi)\\ ∫abf(x)dx=2h[f(a)+2×f(a+2h)+f(b)]−(2h)212b−af′′(ξˉ)=S(a,2a+b)+S(2a+b,b)−48h3f′′(ξˉ)(3)
为了计算误差,我们做出以下假设:
w
e
s
u
p
p
o
s
e
f
′
′
(
ξ
ˉ
)
≈
f
′
′
(
ξ
)
we\ suppose\ f''(\bar \xi) \approx f''(\xi)
we suppose f′′(ξˉ)≈f′′(ξ)
于是有,
S
(
a
,
b
)
−
h
3
12
f
′
′
(
ξ
)
=
S
(
a
,
a
+
b
2
)
+
S
(
a
+
b
2
,
b
)
−
h
3
48
f
′
′
(
ξ
)
(4)
S(a,b)-\frac{h^3}{12}f''(\xi)=S(a,\frac{a+b}{2})+S(\frac{a+b}{2},b)-\frac{h^3}{48}f''(\xi)\\ \tag{4}
S(a,b)−12h3f′′(ξ)=S(a,2a+b)+S(2a+b,b)−48h3f′′(ξ)(4)
S ( a , b ) − [ S ( a , a + b 2 ) + S ( a + b 2 , b ) ] = h 3 12 f ′ ′ ( ξ ) − h 3 48 f ′ ′ ( ξ ) (5) S(a,b)-[S(a,\frac{a+b}{2})+S(\frac{a+b}{2},b)]=\frac{h^3}{12}f''(\xi)-\frac{h^3}{48}f''(\xi)\\ \tag{5} S(a,b)−[S(a,2a+b)+S(2a+b,b)]=12h3f′′(ξ)−48h3f′′(ξ)(5)
S ( a , b ) − [ S ( a , a + b 2 ) + S ( a + b 2 , b ) ] ≈ h 3 16 f ′ ′ ( ξ ) (5) S(a,b)-[S(a,\frac{a+b}{2})+S(\frac{a+b}{2},b)]\approx \frac{h^3}{16}f''(\xi) \tag{5} S(a,b)−[S(a,2a+b)+S(2a+b,b)]≈16h3f′′(ξ)(5)
因此误差可以估计为,
∣
∫
a
b
f
(
x
)
d
x
−
S
(
a
,
a
+
b
2
)
−
S
(
a
+
b
2
,
b
)
∣
≈
1
3
∣
S
(
a
,
b
)
−
S
(
a
,
a
+
b
2
)
−
S
(
a
+
b
2
,
b
)
∣
(6)
|\int_a^bf(x)dx -S(a,\frac{a+b}{2})-S(\frac{a+b}{2},b)|\approx \frac{1}{3}|S(a,b)-S(a,\frac{a+b}{2})-S(\frac{a+b}{2},b)| \tag{6}
∣∫abf(x)dx−S(a,2a+b)−S(2a+b,b)∣≈31∣S(a,b)−S(a,2a+b)−S(2a+b,b)∣(6)
这个结果表明,我们用 S ( a , a + b 2 ) + S ( a + b 2 , b ) S(a,\frac{a+b}{2})+S(\frac{a+b}{2},b) S(a,2a+b)+S(2a+b,b)进行估计积分值的大小的误差为 1 3 ∣ S ( a , b ) − S ( a , a + b 2 ) − S ( a + b 2 , b ) ∣ \frac{1}{3}|S(a,b)-S(a,\frac{a+b}{2})-S(\frac{a+b}{2},b)| 31∣S(a,b)−S(a,2a+b)−S(2a+b,b)∣,我们就可以用函数本身的性质得到计算误差的大小了,故称之为自适应(adapted)。
2、自适应复合梯形公式的matlab 代码
定义精度要求为 1 × 1 0 − 3 1\times 10^{-3} 1×10−3
format long
f1 = @(x) (100/x^2)*(sin(10/x)); % 测试函数(the test function)
% 积分上下限和精度要求,注意这里是整体的精度要求
a=1;
b=3;
e=1e-3;
p=[a,b];
p0=p;
ep=[e];
m=0; %分点数
q=0;
I=0; %最后的结果
while(1)
n1=length(ep);
n=length(p0);
if n==1
break;
end
h=p0(2)-p0(1);
whole=h/2*(f1(p0(1))+f1(p0(1)+h));
left=h/4*(f1(p0(1))+f1(p0(1)+h/2));
right=h/4*(f1(p0(1)+h/2)+f1(p0(1)+h));
if abs(whole-left-right)<=3*ep(1) %精度满足要求,则该部分不再细分
I=I+left+right;
p0=p0(2:n);
if n1>=2
ep=ep(2:n1);
end
q=q+1;
else %精度不满足要求,继续细分
m=m+1; %分层数++
p0=[p0(1),p0(1)+h/2,p0(2:n)];
if n1==1
ep=[ep(1)/2,ep(1)/2]; %第一次分段,将原有的精度要求提高,两段分别不超过ep/2
else
ep=[ep(1)/2,ep(1)/2,ep(2:n1)];%分段,将原有的精度提高,两段分别不超过ep/2
end
if q==0
p=p0;
else
p=[p(1:q),p0];
end
end
end
fprintf('求得的解为:%f\n',I);
fprintf('分层数为%d\n',m);
p=p';
disp('分点为');p
![image-20201204233247643](https://i-blog.csdnimg.cn/blog_migrate/27858d2bdd57502341072cadd0cc4a6f.png)
3、与一般复化梯形公式的比较
clc;
clear;
f=@(x) (100/x^2)*(sin(10/x));
result = -1.426014;
a=input('Enter lower limit a: '); % exmple a=1
b=input('Enter upper limit b: '); % exmple b=3
e=input('Enter the eps: '); % exmple eps=10^-3
n=2;
MAXN = 3000; %define the maximum number of segments
while n<MAXN
sum=0;
h=(b-a)/n;
for k=1:1:n-1
x(k)=a+k*h;
y(k)=f(x(k));
sum=sum+y(k);
end
% Formula: (h/2)*[(y0+yn)+2*(y2+y3+..+yn-1)]
answer=h/2*(f(a)+f(b)+2*sum);
if (abs(answer-result)<=e) %the precision of the result is satisfied.
break;
end
n = n + 1;
end
fprintf('\n The value of integration is %f',answer);
fprintf('\n The number of intermediate points is %f',n);
结论:
可以看出,自适应梯形公式计算的点数要比一般复化梯形公式少,从这个角度上看,自适应梯形公式更加适合计算该函数的数值积分。