课后作业1

Q1

问题

考虑下面程序

x=1
delta=1/2^(53);
for j=1:2^(20)
    x=x+delta
end

根据等比数列求和公式，可知 $x=1+2^{-33}$ ，使用Matlab计算，实际结果是多少？

程序

x=1;
delta=1/2^(53);
for j=1:2^(20)
    x=x+delta;
end
disp(x);

输出

1
2

>> q1
     1

分析

一个64位的双精度浮点数在内存中被分为三个部分：

符号位 : 1 bit. 决定是正数还是负数。
指数位 : 11 bits. 决定了数的大小范围。
尾数位 : 52 bits. 决定了数的精度。

一个浮点数的值大致可以表示为：符号 * (1 + 尾数) * 2^(指数 - 偏置值)

对于一个浮点数系统，存在一个最小的数 ε，使得 1 + ε > 1。任何比 ε 更小的正数 d，在计算 1 + d 时，结果都会被舍入回 1。

这个 ε 在Matlab中可以用 eps 函数查看。对于数值 1，它的机器精度是：
eps(1) 的值为 2^-52 (大约是 2.2204e-16)。

>> eps(1)

ans =

     2.220446049250313e-16

这意味着，对于数字 1 来说，能够分辨的最小的增量就是 2^-52。任何小于这个值的增量（严格来说是小于 eps(1)/2）在与1相加时，都会因为舍入误差而被“吞掉”。

若要完成这次计算，则需要改变计算顺序

x=1;
delta=1/2^(53);
sum=0;
for j=1:2^(20)
    sum=sum+delta;
end
x=x+sum;
disp(x);

1 2	>> q1 1.000000000116415

‍

Q2

问题

当 $x>0$ 时，我们知道 $2x>x$ ，对浮点数该不等式是否成立？运行如下程序：

x=1
twox=2*x
k=0;
while(twox>x)
	x=twox;
	twox=2*x;
	k=k+1;
end

程序

clear;
x=1;
twox=2*x;
k=0;
while(twox>x)
	x=twox;
	twox=2*x;
	k=k+1;
end

输出

分析

在数学上，当 $x>0$ 时， $2x>x$ 永远成立，但从运行结果图上可以看到k为有限值1024。

这是由于在计算机中，浮点数是个有限值，大约为1.7977e+308，当计算超过这个值后，会发生溢出，Matlab就会将计算结果处理为一个特殊的值——无穷大，表示为Inf。

则当x和twox都溢出时$x=Inf=twox $，此时程序跳出循环。

由Q1中分析可知指数位: 11 bits. 决定了数的大小范围。而 $2^{11}=2048$ ，即双精度浮点数的最大指数可取0-2047。

当指数位存储值为 0 (全0) 时，用于表示0和非规格化数。
当指数位存储值为 2047 (全1) 时，用于表示无穷大 (Inf) 和 NaN (Not a Number)。

那么用于表示有限数的最大指数值为2046，而为了表示负指数，双精度浮点数规定了一个偏置值1023。

则能表示的最大正指数为 $2046-1023=1023$ 。所以循环能够进行1023次，最终k值为1024。

‍

Q3

问题

使用Matlab计算线性方程组

\left[ \begin{matrix} 2.00& 1.00\\ 1.99& 1.00\\ \end{matrix} \right] x=\left[ \begin{array}{c} 1.00\\ -1.00\\ \end{array} \right]

两方程在数学理论上完全等价，用Matlab计算 $x-z$ 的值， $x-z$ 是否为零？

程序

clear;
A=[2.00,1.00;1.99,1.00];
B=[1.00;-1.00];
x=A\B;
C=[0.02,0.01;1.99,1.00];
D=[0.01;-1.00];
z=C\D;
y=x-z;
disp(x);
disp(z);
disp(y);

输出

分析

在数学上x与z是完全想等的，即y的理想值为 $\left[ \begin{array}{c} 0\\ 0\\ \end{array} \right]$ 。而Matlab中实际的输出值为 $\left[ \begin{array}{c} -4.2633e-13\\ 8.5265e-13\\ \end{array} \right]$ 。

有两个因素共同导致了这一结果：

病态系统
1
2
3
4
5
6
7
8
9
10
11
12
13
>> A = [2.00 1.00; 1.99 1.00];
>> cond(A)

ans =

9.960089959930169e+02

>> C = [0.02 0.01; 1.99 1.00];
>> cond(C)

ans =

4.960599998033202e+04
Matlab中的cond()函数能够计算条件数，一般来说条件数大于100时，系统就属于病态了。条件数越大，系统就越病态，解对输入的微小扰动就越敏感。
浮点数精度限制

在求解 $x=A\backslash\mathrm{B}$ 的过程中，由于浮点运算，产生了一个无法避免的微小误差ε₁。

在求解 $z=C\backslash\mathrm{D}$ 的过程中，同样会有一个误差ε₂。

由于计算路径和初始值的不同，ε₁与ε₂是不相等的。

由于系统是病态的，ε₁与ε₂被系统本身放大了，最终得到 $x-z\ne 0$ 。

‍

Q4

问题

运行以下Matlab程序

x=0.1
sum=0
for i=1:100
	sum=sum+x
end

最终结果是否为10？如果不是10，为什么？

程序

clear;
x=0.1;
sum=0;
for i=1:100
	sum=sum+x;
end
disp(sum);

输出

1 2	>> q4 9.999999999999980

分析

计算机使用二进制来存储和处理数字，而二进制无法精确表示十进制的0.1，其在二进制下是一个无限循环小数：0.00011001100110011...，所以计算机只能存储这个无限小数的一个近似值，与原数值存在一个不可避免的微小误差。

而程序中这个存在微小误差的值被连续相加了100次，每一次相加都会导致误差的积累，最终导致结果不为准确的10。

‍

Q5

问题

编程实现Horner‘s算法（秦九韶算法）

（1）计算多项式 $P\left( x \right) =x^7-2x^6-3x^4+4x^3-x^2+6x-1$ 在 $x=2$ 处的值。

（2）计算多项式 $P\left( x \right) =1+x+...+x^{50}$ 在 $x=1.00001$ 的值，并与使用 $Q\left( x \right) =\frac{x^{51}-1}{x-1}$ 的结果进行比较。

（3）计算多项式 $P\left( x \right) =1-x+x^2-x^3+...+x^{98}-x^{99}$ 在 $x=1.0001$ 处的值。

程序

%1
clear;
d=7;
c=[1 -2 0 -3 4 -1 6 -1];
x=2;
y=nest(d,c,x);
fprintf('y=%.15f\n',y)
%%
%2
clear;
d=50;
c=ones(1,51);
x=1.00001;
y1=nest(d,c,x);
y2=(x^51-1)/(x-1);
fprintf('y1=%.15f\n',y1)
fprintf('y2=%.15f\n',y2)
%%
%3
clear;
d=99;
n=100;
idx=0:n-1;
c=2*mod(idx,2)-1;
x=1.0001;
y=nest(d,c,x);
fprintf('y=%.15f\n',y)

输出

（1） $y=261.000000000000000$

（2） $y1=51.012752082749991$ ， $y2=51.012752082745230$

（3） $y=0.005024579817447$

分析

虽然 P(x) 和 Q(x) 在数学上是完全等价的，但在计算机进行浮点数运算时，它们的数值稳定性有着显著的差异。

计算 P(x) 的方法 (霍纳法，得到 y1):
- 霍纳法的计算过程是一系列的乘法和加法：(...((c₅₀*x + c₄₉)*x + c₄₈)...)*x + c₀。
- 在这个问题中，所有的系数都是1，且 x = 1.00001 是一个正数。
- 整个计算过程中只涉及正数的乘法和加法。虽然每一步都有微小的舍入误差，但这些误差是逐渐累积的，不会被急剧放大。因此，霍纳法在计算多项式时通常是数值稳定的。
计算 Q(x) 的方法 (直接套用公式，得到 y2):
- 这个计算涉及到 x⁵¹ - 1 和 x - 1。
- 由于 x = 1.00001 非常接近 1，那么 x⁵¹ 的值也会非常接近 1。
- 计算 x⁵¹ - 1 时，就出现了**“相近数相减”**的情况。
- 这会导致灾难性抵消或称为有效数字损失。当两个非常接近的数相减时，它们前面相同的大部分有效数字都被抵消掉了，结果的相对误差会被急剧放大，从而导致精度严重下降。

‍

Q6

问题

用单精度计算 $x^2-\left( 10^{18}+1 \right) x+10^{18}=0$ 的根，精确解为 $x_1=10^{18}$ ， $x_2=1$ 利用求根公式

x=\frac{-b\pm \sqrt{b^2-4ac}}{2a}

可算出

x_1=\frac{-b+\sqrt{b^2-4ac}}{2a}=10^{18},x_2=\frac{-b-\sqrt{b^2-4ac}}{2a}=0

为什么 $x_2=0$ ，请你用别的方法计算出接近真实解的近似解。

程序

% 单精度计算
a = single(1);
b = single(-(1e18 + 1));
c = single(1e18);

discriminant = b^2 - 4*a*c;
x1_direct = (-b + sqrt(discriminant)) / (2*a);
x2_direct = (-b - sqrt(discriminant)) / (2*a);

fprintf('x1 = %.15g\n', double(x1_direct));
fprintf('x2 = %.15g\n', double(x2_direct));

% 利用根与系数的关系
% x1 * x2 = c/a, 所以 x2 = c/(a*x1)
x1_stable = (-b + sqrt(discriminant)) / (2*a);
x2_stable = c / (a * x1_stable);

fprintf('x1 = %.15g\n', double(x1_stable));
fprintf('x2 = %.15g\n', double(x2_stable));

输出

>> q6
x1 = 9.99999984306749e+17
x2 = 0
x1 = 9.99999984306749e+17
x2 = 1

分析

直接求根公式计算 $x_2$ 会得到0，这是因为在单精度下 $b^2$ 和 $4ac$ 的值非常大且非常接近，导致相减时出现灾难性抵消。

因此要想计算结果更接近近似解，利用根与系数的关系 $x_1\times x_2=\frac{c}{a}$ 将 $b^2$ 与 $4ac$ 这两个接近值直接相减避开即可。

‍

Q7

问题

求积分 $I_n=\int_0^1{\frac{x^n}{x+4}}dx$ 。

不稳定算法

\begin{cases} I_n=\frac{1}{n}-4I_{n-1},& n=1,2,\cdots\\ I_0=0.2231(\approx ln1.25)& \\ \end{cases}

稳定算法

\begin{cases} I_{n-1} = \frac{1}{4n}(1 - nI_n), & n = 11, 10, \cdots, 1, 0. \\ I_{11} = \frac{3}{160} \end{cases}

请比较两种算法的结果，并进行误差分析。

程序

clear; clc; format long;

% 精确值计算函数（使用符号计算作为参考）
syms x;
I_exact = zeros(12,1);
for n = 0:11
    I_exact(n+1) = double(int(x^n/(x+4), x, 0, 1));
end

fprintf('精确值:\n');
for n = 0:11
    fprintf('I_%d = %.15f\n', n, I_exact(n+1));
end

%% 不稳定算法（正向递推）
I_unstable = zeros(12,1);
I_unstable(1) = 0.2231;

for n = 1:11
    I_unstable(n+1) = 1/n - 4 * I_unstable(n);
    fprintf('I_%d = %.15f\n', n, I_unstable(n+1));
end

%% 稳定算法（反向递推）
I_stable = zeros(12,1);
I_stable(12) = 3/160; % I_11 = 3/160

for n = 11:-1:1
    I_stable(n) = (1 - n * I_stable(n+1)) / (4 * n);
end

for n = 0:11
    fprintf('I_%d = %.15f\n', n, I_stable(n+1));
end

%% 误差分析
fprintf('\n=== 误差分析 ===\n');
fprintf('n\t不稳定算法误差\t\t稳定算法误差\t\n');
for n = 0:11
    error_unstable = abs(I_unstable(n+1) - I_exact(n+1));
    error_stable = abs(I_stable(n+1) - I_exact(n+1));
    fprintf('%d\t%.2e\t\t%.2e\t\n', n, error_unstable, error_stable);
end
%% 计算条件数
fprintf('\n条件数分析:\n');
for n = 1:10
    cond_unstable = 4^n;  % 不稳定算法的误差放大因子
    cond_stable = (1/4)^(11-n);  % 稳定算法的误差缩小因子
    fprintf('n=%d: 不稳定算法条件数=%.2e, 稳定算法条件数=%.2e\n', n, cond_unstable, cond_stable);
end

输出

精确值:
I_0 = 0.223143551314210
I_1 = 0.107425794743161
I_2 = 0.070296821027356
I_3 = 0.052146049223909
I_4 = 0.041415803104364
I_5 = 0.034336787582543
I_6 = 0.029319516336493
I_7 = 0.025579077511171
I_8 = 0.022683689955316
I_9 = 0.020376351289848
I_10 = 0.018494594840608
I_11 = 0.016930711546657
I_1 = 0.107600000000000
I_2 = 0.069600000000000
I_3 = 0.054933333333334
I_4 = 0.030266666666665
I_5 = 0.078933333333340
I_6 = -0.149066666666695
I_7 = 0.739123809523924
I_8 = -2.831495238095695
I_9 = 11.437092063493889
I_10 = -45.648368253975555
I_11 = 182.684382106811313
I_0 = 0.223143550880458
I_1 = 0.107425796478170
I_2 = 0.070296814087321
I_3 = 0.052146076984048
I_4 = 0.041415692063809
I_5 = 0.034337231744763
I_6 = 0.029317739687613
I_7 = 0.025586184106692
I_8 = 0.022655263573232
I_9 = 0.020490056818182
I_10 = 0.018039772727273
I_11 = 0.018750000000000

=== 误差分析 ===
n	不稳定算法误差		稳定算法误差	
0	4.36e-05		4.34e-10	
1	1.74e-04		1.74e-09	
2	6.97e-04		6.94e-09	
3	2.79e-03		2.78e-08	
4	1.11e-02		1.11e-07	
5	4.46e-02		4.44e-07	
6	1.78e-01		1.78e-06	
7	7.14e-01		7.11e-06	
8	2.85e+00		2.84e-05	
9	1.14e+01		1.14e-04	
10	4.57e+01		4.55e-04	
11	1.83e+02		1.82e-03	

条件数分析:
n=1: 不稳定算法条件数=4.00e+00, 稳定算法条件数=9.54e-07
n=2: 不稳定算法条件数=1.60e+01, 稳定算法条件数=3.81e-06
n=3: 不稳定算法条件数=6.40e+01, 稳定算法条件数=1.53e-05
n=4: 不稳定算法条件数=2.56e+02, 稳定算法条件数=6.10e-05
n=5: 不稳定算法条件数=1.02e+03, 稳定算法条件数=2.44e-04
n=6: 不稳定算法条件数=4.10e+03, 稳定算法条件数=9.77e-04
n=7: 不稳定算法条件数=1.64e+04, 稳定算法条件数=3.91e-03
n=8: 不稳定算法条件数=6.55e+04, 稳定算法条件数=1.56e-02
n=9: 不稳定算法条件数=2.62e+05, 稳定算法条件数=6.25e-02
n=10: 不稳定算法条件数=1.05e+06, 稳定算法条件数=2.50e-01

分析

从输出结果可以看到，不稳定算法的误差比稳定算法误差差了5个数量级，不稳定算法随着计算阶数得上升误差变得越来越大。

可以与用Q3“条件数越大，系统就越病态，解对输入的微小扰动就越敏感”来分析

首先算出误差因子：

不稳定算法的误差传播

递推关系：

I_n = \frac{1}{n} - 4I_{n-1}

误差分析：
设真实值为 $I_n$ ，计算值为

\tilde{I}_n = I_n + \varepsilon_n

代入递推关系：

\tilde{I}_n = \frac{1}{n} - 4\tilde{I}_{n-1}

I_n + \varepsilon_n = \frac{1}{n} - 4(I_{n-1} + \varepsilon_{n-1})

I_n + \varepsilon_n = \frac{1}{n} - 4I_{n-1} - 4\varepsilon_{n-1}

由于 $I_n = \frac{1}{n} - 4I_{n-1}$ ，两边相减得到：

\varepsilon_n = -4\varepsilon_{n-1}

因此误差传播为：

|\varepsilon_n| = 4|\varepsilon_{n-1}| = 4^n|\varepsilon_0|

稳定算法的误差传播

递推关系：

I_{n-1} = \frac{1 - nI_n}{4n}

误差分析：
设真实值为 $I_n$ ，计算值为 $\tilde{I}_n = I_n + \varepsilon_n$

代入递推关系：

\tilde{I}_{n-1} = \frac{1 - n\tilde{I}_n}{4n}

I_{n-1} + \varepsilon_{n-1} = \frac{1 - n(I_n + \varepsilon_n)}{4n}

I_{n-1} + \varepsilon_{n-1} = \frac{1 - nI_n}{4n} - \frac{n\varepsilon_n}{4n}

由于 $I_{n-1} = \frac{1 - nI_n}{4n}$ ，两边相减得到：

\varepsilon_{n-1} = -\frac{1}{4}\varepsilon_n

因此误差传播为：

|\varepsilon_{n-1}| = \frac{1}{4}|\varepsilon_n| = \left(\frac{1}{4}\right)^{11-n}|\varepsilon_{11}|

则可知不稳定算法的误差因子为4，稳定算法的误差因子为 $\frac{1}{4}$ 。最终算出的条件数如输出所示，不稳定算法的条件数随着递推阶数提高而变大，导致误差越来越大；而稳定算法在递推过程中，随着条件数变小，误差也越来越小。

‍

Q8

问题

请计算 $\sin x$ 的10阶、20阶和100阶的泰勒展开式并画出图形。

程序

clear; clc; close all;
% 定义符号变量
syms z;

% 10阶泰勒展开
taylor10 = taylor(sin(z), z, 'Order', 11);
% 20阶泰勒展开
taylor20 = taylor(sin(z), z, 'Order', 21);
% 100阶泰勒展开
taylor100 = taylor(sin(z), z, 'Order', 101);

% 转换为数值函数用于绘图
taylor10_func = matlabFunction(taylor10);
taylor20_func = matlabFunction(taylor20);
taylor100_func = matlabFunction(taylor100);

% 绘制实数轴上的比较
figure('Position', [100, 100, 1200, 500]);

x_real = linspace(-3*pi, 3*pi, 1000);

% 绘制函数值比较
subplot();
plot(x_real, sin(x_real), 'k-', 'LineWidth', 3, 'DisplayName', 'sin(x) 精确值');
hold on;
plot(x_real, taylor10_func(x_real), 'r--', 'LineWidth', 1.5, 'DisplayName', '10阶泰勒展开');
plot(x_real, taylor20_func(x_real), 'g--', 'LineWidth', 1.5, 'DisplayName', '20阶泰勒展开');
plot(x_real, taylor100_func(x_real), 'b--', 'LineWidth', 1.5, 'DisplayName', '100阶泰勒展开');
xlabel('x'); ylabel('f(x)');
title('sin(x)的泰勒展开式比较');
legend('Location', 'northwest');
grid on;

% 设置函数值图的y轴范围
ylim([-2, 2]);

% 添加网格线
set(gca, 'YMinorGrid', 'on');

输出

分析

sin(x) = \sum_{n=0}^{\infty} \frac{(-1)^n x^{2n+1}}{(2n+1)!} = z - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \cdots

由输出图可见，随着泰勒展开的阶数不断提高，拟合曲线在整个定义域上快速趋近于 $\sin x$ ，尤其在 $\left| x \right|$ 较大的区间，拟合效果的改善尤为明显。

数据分析课后作业1

课后作业1

Q1

问题

程序

输出

分析

Q2

问题

程序

输出

分析

Q3

问题

程序

输出

分析

Q4

问题

程序

输出

分析

Q5

问题

程序

输出

分析

Q6

问题

程序

输出

分析

Q7

问题

程序

输出

分析

Q8

问题

程序

输出

分析