发布于  更新于 

数值分析作业 - 浮点运算和误差

数值分析

Problem 1

说明

分析单精度计算 fl(9.4)−fl(9)−fl(0.4) 的结果，并进行计算机实践。

Symbol Exponent mantissa
9.4f + 0 3 + 127 = 10000010 100101100110011001100110
9.0f + 0 3 + 127 = 10000010 100100000000000000000000
9.4f-9.0f + 0 3 + 127 = 10000010 000001100110011001100110
Lsh 5 + 0 -2 + 127 = 01111101 110011001100110011000000
0.4f + 0 -2 + 127 = 01111101 110011001100110011001101
9.4f-9.0f-0.4f - 1 -2 + 127 = 01111101 000000000000000000001101
Lsh 20 - 1 -22 + 127 = 01101001 110100000000000000000000

得到计算结果为

编写辅助函数用于输出浮点数的比特

1
2
3
4
5
6
7
8
9
10
11
12
13
#define print_fpbits(x) print_fpbits_impl(x, #x)
template <typename T>
void print_fpbits_impl(T value, std::string expr) {
    auto y = reinterpret_cast<std::uint8_t*>(&value);
    std::cout << expr << ": " << std::endl;
    for (int i = sizeof(T) - 1; i >= 0; i--) {
        char byte = y[i];
        for (int j = 7; j >= 0; j--) {
            std::cout << ((byte >> j) & 1);
        }
    }
    std::cout << std::endl;
}

验证计算过程

1
2
3
4
5
6
print_fpbits(9.4f);
print_fpbits(9.0f);
print_fpbits(0.4f);
print_fpbits(9.4f - 9.0f);
print_fpbits(9.4f - 9.0f - 0.4f);
std::cout << std::format("{:.12f}\n", 9.4f - 9.0f - 0.4f);

得到输出

1
2
3
4
5
6
7
8
9
10
11
9.4f:
01000001000101100110011001100110
9.0f:
01000001000100000000000000000000
0.4f:
00111110110011001100110011001101
9.4f - 9.0f:
00111110110011001100110011000000
9.4f - 9.0f - 0.4f:
10110100110100000000000000000000
-0.000000387430

Problem 2

说明

设计高效的多项式算法 , 计算机编程比较直接算法与优化算法的计算时间（采用双精度进行计算，分别循环 次）

容易写出直接的 多项式计算方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
const double coef[] = {1, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 4, 0, 0, 0, 5};

double pow_naive(double x, int n) {
    double res = 1.0;
    for (int i = 1; i <= n; i++) {
        res *= x;
    }
    return res;
}

double f_naive(double x) {
    double res = 0.;
    for (int i = 0; i <= 15; i++) {
        res += coef[i] * pow_naive(x, i);
    }
    return res;
}

和简单的 优化

1
2
3
4
5
6
7
8
9
double f_fast(double x) {
    double res = 0.;
    double pow = 1.;
    for (int i = 0; i <= 15; i++) {
        res += coef[i] * pow;
        pow *= x;
    }
    return res;
}

需要注意的是, 针对 C++ 一类编译器具有强大静态优化能力的编程语言, 将简单的算法在同样的输入下单纯地重复运行并不是合理的性能测量方式, 因为编译器很容易将整个循环删除掉. 简单起见, 可以在循环中要求从 std::vector 中读取并存储数据, 以防止无用的循环被删除. 下面报告程序在 MSVC 19.33.31629.0 下, 简单循环 次的用时

1
2
3
4
x = 2, naive: 2272.14ms
x = 2, fast: 418.473ms
x = 2.22... naive: 2245.62ms
x = 2.22..., fast: 426.92ms

启用 /O2 优化后, 性能有显著的提升

1
2
3
4
x = 2, naive: 686.996ms
x = 2, fast: 113.86ms
x = 2.22... naive: 670.233ms
x = 2.22..., fast: 108.462ms

Problem 3

说明

推导计算 的递推公式,用正向和逆向递推计算 ，比较并分析两种迭代的误差和稳定性。

容易推导公式

初始值可以取

这个正向递推公式每次递推将误差放大五倍

1
2
forward = (NestList[{E - #[[2]] #[[1]], #[[2]] + 1} &, {1., 2}, 20] // Transpose)[[1]];
forward // InputForm
1
2
3
4
5
6
7
{1., 0.7182818284590451, 0.5634363430819098,
 0.4645364561314058, 0.395599547802016, 0.34468454164694906,
 0.30549003693040166, 0.27436153301583177, 0.24902803131655915,
 0.22800151529345358, 0.21026516023105568, 0.19509990568637692,
 0.18198305453614516, 0.17051906495301283, 0.1604958541638526,
 0.15034816183740363, 0.16236307722318344, -0.2042535615582568,
 6.599099498065924, -129.26370813285942, 2717.256152618507}

注意到在 16 项处由于误差的累积导致积分有界性已被破坏, 之后的迭代过程明显不收敛.

考虑反向迭代公式

取 反向迭代

1
2
3
backward = (NestList[{1/#[[2]] (E - #[[1]]), #[[2]] - 1} &, {1., 30},
      29] // Transpose)[[1]] // Reverse;
backward // InputForm
1
2
3
4
5
6
7
8
9
10
{1., 0.7182818284590452, 0.5634363430819095,
 0.4645364561314071, 0.3955995478020096, 0.3446845416469873,
 0.30549003693013366, 0.27436153301797606, 0.24902803129726034,
 0.2280015154864418, 0.21026515810818536, 0.19509993116082064,
 0.18198272336837695, 0.17052370130176742, 0.16042630893253398,
 0.15146088553850115, 0.14344677430452524, 0.13623989097759062,
 0.12972389988482333, 0.12380383076257832, 0.11840138244490007,
 0.1134514146712435, 0.10889929102044421, 0.10469884396838404,
 0.10081072924944416, 0.09720286797349674, 0.09380439317463288,
 0.09175881956932448, 0.057276060948634834, 1.}

与 NIntegrate 直接积分 (保证 20 位有效数字) 的结果比较并绘制图像

1
2
3
4
direct =
 Table[NIntegrate[x^n E^x, {x, 0, 1}, AccuracyGoal -> 20,
   WorkingPrecision -> 20], {n, 1, 20}];
ListLinePlot[{forward, backward, direct}]

绘制两种算法和直接数值积分法的残差

1
2
ListLinePlot[
 Abs[{Take[forward, 20] - direct, Take[backward, 20] - direct}]]

由于反向迭代逐渐减小误差的特性, 表现远好于正向迭代

许可协议

本文采用 署名-相同方式共享 4.0 国际 许可协议，转载请注明出处。

分享文章
快来参与讨论吧