辛普森求积公式 和 复合辛普森求积公式 Matlab 实现

辛普森求积公式 和 复合辛普森求积公式 Matlab 实现

辛普森求积公式

利用区间[a,b]的端点及中点计算积分

$$\int_{x_1}^{x_2} f(x)dx \approx \frac{b-a}{6} \times (f(a) + 4 \times f(\frac{b-a}{2}) -+f(b))$$

复合辛普森求积公式

将求积区间[a,b]分为n个区间，每个区间步长为h($h = \frac{b-a}{n}$)在每个区间求辛普森积分

$$\int_{x_1}^{x_2} f(x)dx \approx \frac{h}{6} \int (f(a) + 4 \times\sum_{k = 0}^{n-1}f(x_\frac{k+1}{2}) + 2 \times\sum_{k = 1}^{n-1}f(x_k) + f(b))$$
其中 $x_\frac{k+1}{2} = x_k + \frac{h}{2}$

体会

$$\sum_{k = 0}^{n-1}f(x_\frac{k+1}{2})$$ 这项是复合辛普森求积公式比复合梯形求积公式多出的式子
这项是辛普森求积公式多余出计算量也正是 复合更精确的原因

说明

上述公式是我用mathjax写的，如有错误请联系我修正
敬请指正
概述省略了部分推导过程，请查阅详细推导资料

Matlab 实现代码

辛普森求积公式

将该函数存为m文件

function res = Simpson(f,a,b)
    format long;
    if b < a
        c = b;
        b = a;
        a = c;
    end
    res = (b-a)/6 * (f(a) + 4*f((a+b)/2) + f(b));
end

调用下面语句测试函数

f = inline('sin(x)','x')
Trapezium(f,0,pi/2)

复合辛普森求积公式

将该函数存为m文件

function res = ComSimpson(f,n,a,b)
    format long;
    if b < a
        c = b;
        b = a;
        a = c;
    end
    h = (b-a)/n;
    d = f(a);
    for i = a+h:h:b-h 
        d = d + (2 * f(i));
    end
    for i = a+h/2:h:b-h/2 
        d = d + (4 * f(i));
    end
    d = d + f(b);
    res = (d * h / 6);

调用下面语句测试函数

f = inline('sin(x)','x')
ComSimpson(f,4,0,pi/2)