数值分析实验 - 数值积分与数值微分

发表于 更新于 分类于 ， ， ， 本文字数： 1.3k 阅读时长 ≈ 5 分钟

数值分析实验 4 - 数值积分与数值微分

目的和意义

实验目的

选用复合梯形公式，复合 Simpson 公式，Romberg 算法，计算

  1. \(\displaystyle I=\int_0^\frac{1}{4}\sqrt{4-\sin^2x}\mathrm{d}x\)
  2. \(\displaystyle I=\int_0^1\frac{\sin x}{x}\mathrm{d}x\)
  3. \(\displaystyle I=\int_0^1\frac{e^x}{4+x^2}\mathrm{d}x\)
  4. \(\displaystyle I=\int_0^1\frac{\ln(1+x)}{1+x^2}\mathrm{d}x\)

要求

  1. 编制数值积分算法的程序
  2. 分别用两种算法计算同一个积分，并比较其结果
  3. 分别取不同步长 \(n\), 试比较计算结果 (如 \(n=10,20\) 等)
  4. 给定精度要求 \(\epsilon\), 试用变步长算法，确定最佳步长

意义

  1. 深刻认识数值积分法的意义
  2. 明确数值积分精度与步长的关系
  3. 根据定积分的计算方法，可以考虑二重积分的计算问题

计算公式

梯形积分

\[ \int_a^bf(x)\mathrm{d}x\approx\frac{h}{2}\left(f(a)+2\sum_{i=1}^{n-1}f(x_k)+f(b)\right) \]

其中

  • 步长 \(h=\frac{b-a}{n}\)
  • 等分点 \(x_k=a+kh,~k=0,1,2,...,n\); \(x_0=a,x_n=b\)

Simpson 积分

\[ \int_a^bf(x)\mathrm{d}x\approx\frac{h}{6}\left(f(a)+4\sum_{i=1}^{n-1}f(x_{k+1/2})+2\sum_{i=1}^{n-1}f(x_k)+f(b)\right) \]

其中

  • 步长 \(h=\frac{b-a}{n}\)
  • 等分点 \(x_k=a+kh,~k=0,1,2,...,n\); \(x_0=a,x_n=b\)
  • \(x_{k+1/2}={x_k+x_{k+1}\over2},~k=0,1,...,n-1\)

Romberg 算法 (复化梯形)

  • Richardson 外推法

    对于

    \[ T_n=\frac{h}{2}\left(f(a)+2\sum_{i=1}^{n-1}f(x_k)+f(b)\right)\tag{1} \]

    有

    \[ T_n=I+\tau_1h^2+\tau_2h^4+O(h^6) \]

    其中

    • \(I=\displaystyle\int_a^bf(x)\mathrm{d}x\)
    • \(\tau_1,\tau_2\) 为与 \(h\) 无关的常数

    注意到

    \[ T_{2n}=I+\tau_1\left(\frac{h}{2}\right)^2+\tau_2\left(\frac{h}{2}\right)^4+O(h^6)\tag{2} \]

    \((1)-4\times (2)\), 得

    \[ S_n:={4T_{2n}-T_n\over3}=I-\frac{\tau_2}{4}h^4+O(h^6) \]

  • Romberg 算法

    记

    • \(h_{k+j}=\frac{1}{2}h_{k+j-1}\)
    • \(T_{0,0}=\frac{b-a}{2}(f(a)+f(b))\)
    • \(\displaystyle T_{0,k+1}=\frac{1}{2}T_{0,k}+\frac{h_k}{2}\sum_{i=0}^{n-1}f(x_{i+1/2})\)
    • \(\displaystyle T_{j,k}={4^jT_{j-1,k+1}-T_{j-1,k}\over4^j-1}\)

    则 \(T_{n,0}\) 即为答案

关于多重积分

只需嵌套进行数值积分即可

一些特殊处理

  • 由于积分 2 为瑕积分，故编写程序时将积分区间取为 [eps, 1]
  • 根据文献 1 的结果，编写自适应积分程序时，对自适应梯形积分和自适应 Simpson 积分进行了常数优化。如对自适应 Simpson 积分，终止条件改为 \(|S-(S_l+S_r)|<15\epsilon\), 返回值改为 \(S_l+S_r+\frac{S-(S_l+S_r)}{15}\)

程序设计

主程序

Show code

main.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
% Exp.4

% @Author: Tifa
% @LastEditTime: 2021-05-06 23:06:44

% Global arguments
now_type = 2;
round_digit = 12;
if now_type == 1
    now_integral = @integral_with_h;
    now_arg = 4;
    get_step = false;
elseif now_type == 2
    now_integral = @integral_with_eps;
    now_arg = 1e-10;
    get_step = true;
end

% Data
syms x
f = {@(x) sqrt(4 - sin(x).^2), @(x) sin(x) ./ x, @(x) exp(x) ./ (4 + x.^2), @(x) log(1 + x) ./ (1 + x.^2)};
ranges = [0, 0.25; 0 + eps, 1; 0, 1; 0, 1];
method_key = ['tpz'; 'sps'; 'rbg'];
method_name = {'Trapezoid'; 'Simpson'; 'Romberg'};

disp(['Now : ', func2str(now_integral)])

for i = 1:length(method_key)
    disp(['Method: ', method_name{i}])

    for j = 1:length(f)
        now_funcstr = regexp(func2str(f{j}), '^@\([a-z,]*[a-z]\)', 'split');
        fprintf(['\tFunction: ', now_funcstr{end}, '; from %.2f to %.2f\n'], ranges(j, 1), ranges(j, 2))

        if get_step
            [val, step] = now_integral(f{j}, method_key(i, :), ranges(j, 1), ranges(j, 2), now_arg);
            now_vals = [val, step];
        else
            now_vals = now_integral(f{j}, method_key(i, :), ranges(j, 1), ranges(j, 2), now_arg);
        end

        fprintf('\tValue: ')
        disp(vpa(now_vals, round_digit))
    end

end

输入数据检查

Show code

input_check.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
function input_check(input_function, integral_type, l, r, h_or_epsi)
    % Input check of Exp.4

    % @Author: Tifa
    % @LastEditTime: 2021-05-06 23:06:44

    if ~isnumeric(l) ||~isnumeric(r) ||~isnumeric(h_or_epsi)
        error('l, r, h and eps should be numerals')
    end

    if h_or_epsi <= 0
        error('h and eps should be positive')
    end

    if ~ischar(integral_type)
        error('Integral type should be a character array')
    end

    if ~isa(input_function, 'function_handle')
        error('Input function should be function handle')
    end

    if nargin(input_function) ~= 1
        error('Input function should be only one variant')
    end

end

按步长积分 (复化积分)

Show code

integral_with_h.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
function output = integral_with_h(input_function, integral_type, l, r, h)
    input_check(input_function, integral_type, l, r, h)

    % @Author: Tifa
    % @LastEditTime: 2021-05-06 23:06:44

    if ~isinteger(h) && h ~= round(h)
        h = round(h);
        fprintf(2, 'Warning: input h is not integer, using h = %d\n', h)
    end

    x = linspace(l, r, 2^h + 1);

    if l <= r
        output = 1;
    else
        output = -1;
        tmp = l; l = r; r = tmp;
    end

    switch integral_type
        case 'sps'
            output = output .* vpa(sum(simpson(input_function, x(1:end - 1), x(2:end))));
        case 'tpz'
            output = output .* vpa(sum(trapez(input_function, x(1:end - 1), x(2:end))));
        case 'rbg'
            output = output .* vpa(romberg(input_function, l, r, h));
        otherwise
            error('Unknown integral type')
    end

end

function output = romberg(inputf, l, r, h)
    n = 1; len = r - l;
    T = zeros(h + 1);

    T(1, 1) = len / 2 * (inputf(l) + inputf(r));

    for k = 2:h + 1
        len = len / 2;
        T(k, 1) = T(k - 1, 1) / 2 + len * sum(inputf(l + (2 * (1:n) - 1) .* len));

        for j = 1:k - 1
            T(k, j + 1) = T(k, j) + (T(k, j) - T(k - 1, j)) / (4^j - 1);
        end

        n = n * 2;
    end

    output = T(h + 1, h + 1);
end

按精度积分 (自适应积分)

Show code

integral_with_eps.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
function [val, min_step] = integral_with_eps(input_function, integral_type, l, r, epsi)
    input_check(input_function, integral_type, l, r, epsi)

    % @Author: Tifa
    % @LastEditTime: 2021-05-06 23:06:44

    if l <= r
        flag = 1;
    else
        flag = -1;
        tmp = l; l = r; r = tmp;
    end

    switch integral_type
        case 'sps'
            func = @asr;
        case 'tpz'
            func = @atr;
        case 'rbg'
            func = @romberg;
        otherwise
            error('Unknown integral type')
    end

    [val, min_step] = func(input_function, l, r, epsi);
    val = val * flag;

end

function [val, min_step] = asr(inputf, l, r, epsi)
    [val, min_step] = BASE(inputf, l, r, epsi, @simpson, 15);
end

function [val, min_step] = atr(inputf, l, r, epsi)
    [val, min_step] = BASE(inputf, l, r, epsi, @trapez, 3);
end

function [val, min_step] = BASE(inputf, l, r, epsi, integral_function, C, S)
    narginchk(5, 7)

    if nargin == 5
        C = 1;
        S = integral_function(inputf, l, r);
    elseif nargin == 6
        S = integral_function(inputf, l, r);
    end

    mid = l + (r - l) / 2;
    Sl = integral_function(inputf, l, mid); Sr = integral_function(inputf, mid, r);

    if abs(Sl + Sr - S) < C * epsi
        val = vpa(Sl + Sr + (Sl + Sr - S) / C);
        min_step = r - l;
    else
        [vall, stpl] = BASE(inputf, l, mid, epsi / 2, integral_function, C, Sl);
        [valr, stpr] = BASE(inputf, mid, r, epsi / 2, integral_function, C, Sr);
        val = vall + valr;
        min_step = min(stpl, stpr);
    end

end

function [val, min_step] = romberg(inputf, l, r, epsi)
    k = 0; len = r - l;
    T = len * (inputf(l) + inputf(r)) / 2;
    err = inf;

    while err >= epsi
        k = k + 1;
        len = len / 2;
        T(k + 1, 1) = T(k, 1) / 2 + len * sum(inputf(l + (2 * (1:(2^(k - 1))) - 1) .* len));

        for j = 1:k
            T(k + 1, j + 1) = T(k + 1, j) + (T(k + 1, j) - T(k, j)) / (4^j - 1);
        end

        err = abs(T(k + 1, k + 1) - T(k, k));
    end

    val = vpa(T(k + 1, k + 1));
    min_step = 1/2^k;
end

辅助函数

梯形积分

Show code

trapez.mview raw
1
2
3
4
5
6
7
function output = trapez(inputf, l, r)

    % @Author: Tifa
    % @LastEditTime: 2021-05-06 23:06:44

    output = (inputf(l) + inputf(r)) .* (r - l) / 2;
end

Simpson 积分

Show code

simpson.mview raw
1
2
3
4
5
6
7
function output = simpson(inputf, l, r)

    % @Author: Tifa
    % @LastEditTime: 2021-05-06 23:06:44

    output = (inputf(l) + 4 * inputf (l + (r - l) / 2) + inputf(r)) .* (r - l) / 6;
end

结果讨论和分析

结果 - 复化积分

  • 分段数 \(n=8\)

    表 1

    积分 1积分 2积分 3积分 4
    参考值0.498711120.946083070.390811850.27219826
    梯形公式0.498701290.945690860.390910990.27076864
    Simpson 公式0.498711120.946083090.390811860.27219871
    Romberg 算法0.498711120.946083070.390811850.27219672

  • 分段数 \(n=16\)

    表 2

    积分 1积分 2积分 3积分 4
    参考值0.498711120.946083070.390811850.27219826
    梯形公式0.498708660.945985030.390836640.27184119
    Simpson 公式0.498711120.946083070.390811850.27219829
    Romberg 算法0.498711120.946083070.390811850.27219827

结果 - 自适应积分

  • 精度要求 \(\epsilon=10^{-8}\)

    • 结果

      表 3

      积分 1积分 2积分 3积分 4
      参考值0.4987111175750.9460830703670.3908118455640.272198261288
      梯形公式0.4987111175750.9460830703670.3908118455640.272198261288
      Simpson 公式0.4987111175740.9460830703670.3908118455620.272198261327
      Romberg 算法0.4987111175750.9460830703670.3908118455560.272198261288

    • 步长

      表 4

      积分 1积分 2积分 3积分 4
      梯形公式1/5121/10241/5121/2048
      Simpson 公式1/81/81/81/16
      Romberg 算法1/81/161/161/64

  • 精度要求 \(\epsilon=10^{-10}\)

    • 结果

      表 5

      积分 1积分 2积分 3积分 4
      参考值0.4987111175750.9460830703670.3908118455640.272198261288
      梯形公式0.4987111175750.9460830703670.3908118455640.272198261288
      Simpson 公式0.4987111175750.9460830703670.3908118455640.272198261288
      Romberg 算法0.4987111175750.9460830703670.3908118455640.272198261288

    • 步长

      表 6

      积分 1积分 2积分 3积分 4
      梯形公式1/81921/163841/81921/32768
      Simpson 公式1/321/161/321/64
      Romberg 算法1/81/161/321/64

分析

  • 三种积分方法中，梯形积分精度最差，其余两种方法所差无几
  • 观察 表 4 与 表 6 发现，积分精度为 \(n\) 位时，三种方法的最小步长分别为: \(O(n^{-4})\), \(O(n^{-2})\), \(O(n^{-2})\)

  1. J N Lyness. Notes on the Adaptive Simpson Quadrature Routine[J]. Journal of the ACM, 1969, 16(3): 483–495↩︎