微积分学习笔记-数值积分

梯型法

为逼近 $\int _a^bf(x)dx$ , 用
$T=\frac{h}{2}(y_0 + 2y_1 + 2y_2 + ... + 2y_{n-1}+y_n)$ .
各y是f在分点
$x_0 = a, x_1 = a + h, x_2 = a+2h,...x_{n-1}, a + (n-1)h, x_n = b$ 的函数值，其中 $h = \frac{(b - a)}{ n}$ .

例1. 使用n = 4时的梯型法估计 $\int _1^2 x^2dx$ , 比较估计值和精确值.
$h = \frac{(2 - 1)}{4} = 0.25$
x和y值如下:

序号	x	y
0	1	1
1	1.25	1.5625
2	1.5	2.25
3	1.75	3.0625
4	2	4

我们有 $T = \frac{h}{2}(y_0 + 2y_1 + 2y_2 + 2y_3 + y_4) = 0.125 * (1 + 2 * 1.5625 + 2 * 2.25 + 2 * 3.0625 + 4) = 2.34375$
积分精确值是:
$\int_1^2 x^2 dx = [\frac{x^3}{3}]_1^2 = \frac{7}{3} = 2.3333333$
两者相对误差为0.446%
梯形法用程序实现比较方便, C++实现如下：

#include <iostream>
using namespace std;
double integrate(std::function<double(double)> f, double a, double b, int n) {
    double h = (b - a) / n, sum = 0, x = a;
    for (int i = 0; i <= n; ++i) {
        x = a + i * h;
        if (i == 0 || i == n) sum += f(x);
        else sum += 2 * f(x);
    }
    return sum * h / 2;
}
int main(int argc, const char * argv[]) {
    auto f = [](double x) { return x * x; };
    cout << integrate(f, 1, 2, 4) << endl;
    return 0;
}

梯形法的误差估计

如果 $f''$ 连续且M是 $|f''|$ 的值在 $[a,b]$ 上的一个上界，则
$|E_T| \leq \frac{b-a}{12}h^2M$ , 其中 $h = \frac{b - a}{n}$ .

例2. 用n=10步长时的体形法估计 $\int _0^\pi x sin(x) dx$ 时带来的误差的上界。

对于 $a = 0, b = \pi$ 和 $h = \frac{b - a}{n} = \frac{\pi}{10}$
$|E_r| \leq \frac{b - a}{12}h^2 M = \frac{\pi^3}{1200} M$
$f''(x) = |2cosx - xsinx|$
$\leq 2 |cosx| + |x||sinx|$
$\leq 2 + \pi * 1 = 2 + \pi$ .
则
$|E_T| \leq \frac{\pi^3}{1200} * (2 + \pi) = 0.1328$ .

Simpson法

为逼近 $\int_a^bf(x)dx$ , 用
$S = \frac{h}{3} (y_0 + 4y_1 + 2y_2 + ... + 2y_{n-2} + 4y_{n-1} + y_n)$ .
各y是f在分点
$x_0 = a, x_1 = a + h, x_2 = a + 2h, ..., x_{n-1} = a + (n - 1)h, x_n = b的值. n是偶数，而h = \frac{(b -a)}{n}$ .

Simpson法的误差估计

若 $f^{(4)}$ 连续，而M是 $f^{(4)}$ 在[a, b] 上的一个上界, 则
$|E_S| \leq \frac{b - a}{180}h^4M$ .
其中h = (b - a) / n.

例4 用n = 4时的Simpson逼近 $\int_0^{2} 5x^4 dx$ , 并计算误差。

$h = \frac{2 - 0} {4} = \frac{1}{2}$

序号	x	y
0	0	0
1	0.5	0.3125
2	1	5
3	1.5	25.3125
4	2	80

则
$S = \frac{h}{3}(y_0 + 4y_1 + 2y_2 + 4y_3 +y_4) = \frac{1}{6}(0 + 4 * 0.3125 + 2 * 5 + 4 * 25.3125 + 80) = 32.083333$
$f^{(4)}(x) = 120$
则
$|E_S| \leq \frac{2 - 0} {180} * \frac{1}{2}^4 * 120 = 0.0833333$
用C++实现Simpson法的代码如下：

#include <iostream>
#include <math.h>
using namespace std;
double integrate(std::function<double(double)> f, double a, double b, int n) {
    double h = (b - a) / n, sum = 0, x = a;
    for (int i = 0; i <= n; ++i) {
        x = a + i * h;
        if (i == 0 || i == n) sum += f(x);
        else sum += (i % 2 ? 4 : 2) * f(x);
    }
    return sum * h / 3;
}
int main(int argc, const char * argv[]) {
    auto f = [](double x) { return 5 * pow(x, 4); };
    cout << integrate(f, 0, 2, 4) << endl;
    return 0;
}

例5 比较梯形法和Simpton法逼近

ln2可以由公式 $\int _1^2\frac{1}{x}dx$ 计算，分别用梯形法和Simpton法计算。
用C++代码计算方法如下:

#include <iostream>
#include <math.h>
using namespace std;
enum IntegrateMethod {
    Trapezoid,
    Simpton
};
double integrate(std::function<double(double)> f, double a, double b, int n, IntegrateMethod method) {
    double h = (b - a) / n, sum = 0, x = a;
    for (int i = 0; i <= n; ++i) {
        x = a + i * h;
        if (i == 0 || i == n) sum += f(x);
        else {
            if (method == Simpton) {
                sum += (i % 2 ? 4 : 2) * f(x);
            } else {
                sum += 2 * f(x);
            }
        }
    }
    return sum * h / (method == Simpton ? 3 : 2);
}
 
int main(int argc, const char * argv[]) {
    auto f = [](double x) { return 1 / x; };
    cout << integrate(f, 1, 2, 100, Trapezoid) << endl;
    cout << integrate(f, 1, 2, 100, Simpton) << endl;
    return 0;
}

最后编辑于：2019.12.02 22:42:57