辛普森积分公式求曲线长度

一、积分的概念

积分（integral）的几何意义是函数的曲线上 x 的一段区间与 x 轴围成的曲边梯形的面积：

x x x的区间为 [ a , b ] [a,b] [a,b]，那么上图阴影面积为： ∫ a b f ( x ) d x int_{a}^{b}f(x)dx ∫abf(x)dx

计算方法一：分割成无穷多个小区间。

∫ a b f ( x ) d x = lim ⁡ x → ∞ ∑ i = 1 n b − a n f ( a + b − a n i ) int_{a}^{b}f(x)dx=limlimits_{x ightarrowinfty}sum limits_{i=1}^{n}{frac{b-a}{n}f(a+frac{b-a}{n}i)} ∫abf(x)dx=x→∞limi=1∑nnb−af(a+nb−ai)

计算方法二：牛顿-莱布尼茨公式。 若 F ′ ( x ) = f ( x ) , 那 么 : 若F(x)=f(x),那么: 若F′(x)=f(x),那么: ∫ a b f ( x ) d x = F ( b ) − F ( a ) int_{a}^{b}f(x)dx=F(b)-F(a) ∫abf(x)dx=F(b)−F(a) 如果容易求出 n n n趋近于无穷大时 f f f的和，可以使用方法一，如 f ( x ) = x 2 f(x)=x^2 f(x)=x2,而这时 f ( x ) = 1 x f(x)=frac{1}{x} f(x)=x1就不适用。如果容易求得 F F F,可以使用方法二， f ( x ) = 1 x f(x)=frac{1}{x} f(x)=x1.但如果两个特点都不满足，那么两种方法都无法使用。于是，我们引入了 数值积分。最常用的就是自适应辛普森积分。

二、辛普森公式

基本思想就是把复杂的函数 f f f近似成二次函数。

∫ a b f ( x ) d x ≈ ∫ a b ( A x 2 + B x + C ) d x = 1 x ( b 3 − a 3 ) + B 2 ( b 2 − a 2 ) + C ( b − a ) = 2 A ( b 3 − a 3 ) + 3 B ( b 2 − a 2 ) + 6 C ( b − a ) 6 = 2 A ( b − a ) ( b 2 + a b + a 2 ) + 3 B ( b + a ) ( b − a ) + 6 C ( b − a ) 6 = ( b − a ) ( 2 A b 2 + 2 A a b + 2 A a 2 + 3 B b + 3 B a + 6 C 6 = ( b − a ) ( A a 2 + B a + C + A b 2 + B b + C + A a 2 + 2 A a b + A b 2 + 2 B b + 2 B a + 4 C 6 = ( b − a ) ( f ( a ) + f ( b ) + 4 f ( a + b 2 ) 6 egin{aligned} &int_{a}^{b}f(x)dxapproxint_{a}^{b}(Ax^2+Bx+C)dx\ &=frac{1}{x}(b^3-a^3)+frac{B}{2}(b^2-a^2)+C(b-a)\ &=frac{2A(b^3-a^3)+3B(b^2-a^2)+6C(b-a)}{6}\ &=frac{2A(b-a)(b^2+ab+a^2)+3B(b+a)(b-a)+6C(b-a)}{6}\ &=frac{(b-a)(2Ab^2+2Aab+2Aa^2+3Bb+3Ba+6C}{6}\ &=frac{(b-a)(Aa^2+Ba+C+Ab^2+Bb+C+Aa^2+2Aab+Ab^2+2Bb+2Ba+4C}{6}\ &=frac{(b-a)(f(a)+f(b)+4f(frac{a+b}{2})}{6} end{aligned} ∫abf(x)dx≈∫ab(Ax2+Bx+C)dx=x1(b3−a3)+2B(b2−a2)+C(b−a)=62A(b3−a3)+3B(b2−a2)+6C(b−a)=62A(b−a)(b2+ab+a2)+3B(b+a)(b−a)+6C(b−a)=6(b−a)(2Ab2+2Aab+2Aa2+3Bb+3Ba+6C=6(b−a)(Aa2+Ba+C+Ab2+Bb+C+Aa2+2Aab+Ab2+2Bb+2Ba+4C=6(b−a)(f(a)+f(b)+4f(2a+b)

b − a b−a b−a 的值越小，上式两边越接近。

三、处理精度问题

有了 S i m p s o n Simpson Simpson 公式，一个自然的想法是把积分区间拆成多个小区间后求和。但是分成区间的个数和长度因积分区间和精度要求甚至被积函数而异。换句话说，分的区间数太少不满足精度要求，太多了会超时。自适应：就是求积分时能够自动控制切割的区间大小和长度，使精度能满足要求。

四、代码实现

4.1 辛普森积分函数

% 辛普森积分求曲线长度相关函数
function sum=AdaptiveSimpson(L,R,err)
    sum=SimpsonBasisFun(L,R,err,Simpson(L,R));
end

function sum=SimpsonBasisFun(L,R,err,LRSum)
    Mid=(L+R)/2.0;
    LSum=Simpson(L,Mid);
    RSum=Simpson(Mid,R);
    if(abs(LSum+RSum-LRSum)<=15.0*err)
        sum=LSum+RSum+(LSum+RSum-LRSum)/15.0;
    else
        sum=SimpsonBasisFun(L,Mid,err/2.0,LSum)+SimpsonBasisFun(Mid,R,err/2.0,RSum);
    end
end

function sum=Simpson(L,R)
    sum=(R-L)*(integralFun(L)+integralFun(R)+4*integralFun((L+R)/2.0))/6.0;%integralFun 为自动生成
end

4.2 测试代码

clc;clear;close all;
%% 1.生成 integralFun 文件
fileL=fopen(integralFun.m,w);
fprintf(fileL, function y=integralFun(u)   

);
fprintf(fileL, y=sqrt(u)+u^5;

);%待求曲线长度的函数表达式
fprintf(fileL, end 

);
fclose(fileL);

%% 2.测试1
AdaptiveSimpson(5,10,0.0001)

%% 3.测试2
figure(1);grid on;hold on;
curLen=0;
for i=0:0.05:10
    if i>0
        curLen=curLen+AdaptiveSimpson(i-0.05,i,0.0001);
    end
    y=sqrt(i)+i^5;
    plot(i,y,k.);
    plot(i,curLen,r.);
    drawnow;
end

4.3效果展示

五、参考资料