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

数值分析实验 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.49871112 0.94608307 0.39081185 0.27219826
    梯形公式 0.49870129 0.94569086 0.39091099 0.27076864
    Simpson 公式 0.49871112 0.94608309 0.39081186 0.27219871
    Romberg 算法 0.49871112 0.94608307 0.39081185 0.27219672
  • 分段数 \(n=16\)

    2

    积分 1 积分 2 积分 3 积分 4
    参考值 0.49871112 0.94608307 0.39081185 0.27219826
    梯形公式 0.49870866 0.94598503 0.39083664 0.27184119
    Simpson 公式 0.49871112 0.94608307 0.39081185 0.27219829
    Romberg 算法 0.49871112 0.94608307 0.39081185 0.27219827

结果 - 自适应积分

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

    • 结果

      3

      积分 1 积分 2 积分 3 积分 4
      参考值 0.498711117575 0.946083070367 0.390811845564 0.272198261288
      梯形公式 0.498711117575 0.946083070367 0.390811845564 0.272198261288
      Simpson 公式 0.498711117574 0.946083070367 0.390811845562 0.272198261327
      Romberg 算法 0.498711117575 0.946083070367 0.390811845556 0.272198261288
    • 步长

      4

      积分 1 积分 2 积分 3 积分 4
      梯形公式 1/512 1/1024 1/512 1/2048
      Simpson 公式 1/8 1/8 1/8 1/16
      Romberg 算法 1/8 1/16 1/16 1/64
  • 精度要求 \(\epsilon=10^{-10}\)

    • 结果

      5

      积分 1 积分 2 积分 3 积分 4
      参考值 0.498711117575 0.946083070367 0.390811845564 0.272198261288
      梯形公式 0.498711117575 0.946083070367 0.390811845564 0.272198261288
      Simpson 公式 0.498711117575 0.946083070367 0.390811845564 0.272198261288
      Romberg 算法 0.498711117575 0.946083070367 0.390811845564 0.272198261288
    • 步长

      6

      积分 1 积分 2 积分 3 积分 4
      梯形公式 1/8192 1/16384 1/8192 1/32768
      Simpson 公式 1/32 1/16 1/32 1/64
      Romberg 算法 1/8 1/16 1/32 1/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↩︎