数值积分解法

介绍

这篇博客讨论的是一个数值积分问题，对于一些难于求积的函数，使用牛顿-莱布尼茨公式显然是不科学的。因此对于这类问题，我们可以从积分中值定理出发，使用矩形或梯形的面积去近似积分值。
这里主要讨论的是复合梯形公式和复合辛普森公式，统称为复合求积法。这种通过把积分区间细分成若干个子区间（通常是等分），再在每个子区间上使用低阶求积公式，从而提高了计算精度。

分析

1.复合梯形公式

该公式核心思想是对于细分后的每一个子区间，使用梯形公式求积。假设将区间$[a,b]$等分为$n$个子区间，分点$x_k = a + kh，h = \frac{b-a}{n}，k=0,1,…,n$，公式的原理就是对于每一块小梯形而言，上底为$f(x_k)$，下底为$f(x_{k+1})$，高为$h$。
复合梯形公式如下：
$$
I = \int^b_af(x)dx = \sum^{n-1}{k=0}\int^{x{k+1}}{x_k}f(x)dx = \frac{h}{2}\sum^{n-1}_{k=0}[f(x_k) + f(x_{k+1})]，
$$
而在计算机编程时，通常使用的是下面这个形式:
$$
T_n = \frac{h}{2}\sum^{n-1}_{k=0}[f(x_k) + f(x_{k+1})] = \frac{h}{2}[f(a)+2\sum^{n-1}{k=1}f(x_k) + f(b)]
$$

2.复合辛普森公式

复合辛普森公式的思路和复合梯形公式的思路一致，差别在于它的低阶求积公式使用了辛普森公式。辛普森公式是经过加权改造的梯形公式，集中在中间的点的权值更高，因此相比起简单的梯形公式，辛普森公式拥有更高的数值精度。
将区间&[a,b]&分为&n&等分，在每个字区间$[x_k, x_{k+1}]$上采用辛普森公式，记$x_{k+1/2} = x_k + \frac{1}{2}h$，得：
$$
I = \int^b_af(x)dx = \sum^{n-1}{k=0}\int^{x{k+1}}{x_k}f(x)dx \ = \frac{h}{6}\sum^{n-1}_{k=0}[f(x_k) + 4f(x_{k+1/2}) + f(x{k+1})]
$$
在计算机编程时，通常写成下列形式：
$$
S_n = \frac{h}{6}[f(a) + 4\sum^{n-1}{k=0}f(x{k+1/2}) + 2\sum^{n-1}_{k=1}f(x_k) + f(b)]
$$

代码实现

1.复合梯形公式

根据上述推到的复合梯形公式，编写matlab代码。

% 复合梯形公式
function [result] = CompositeTrapezoid(a, b, n)
    if (b < a)
        c = b;
        b = a;
        a =c;
    end
    h = (b - a) / n; % 计算步长
    result = myFun(a) + myFun(b);
    for k = 1:n-1
        x = a + k * h;
        result = result + 2 * myFun(x);
    end
    result = (h / 2) * result;
end

function [y] = myFun(x)
    if (x == 0)
        y = 1;
    else
        y = sin(x) / x;
    end
end

2.复合辛普森公式

根据上述推导的复合辛普森公式，编写matlab代码。

% 复合辛普森公式
function [result] = CompositeSimpson(a, b, n)
    if (b < a)
        c = b;
        b = a;
        a = c;
    end
    h = (b - a) / n; % 计算步长
    result = myFun(a) + myFun(b);
    for k = 1:n-1
        x = a + k * h;
        result = result + 2 * myFun(x);
    end
    for k = 1:n
        x = a + (k-1)*h + 1/2 * h;
        result = result + 4 * myFun(x);
    end
    result = result * (h / 6);
end

function [y] = myFun(x)
    if (x == 0)
        y = 1;
    else
        y = sin(x) / x;
    end
end