数值分析实验

实验一 1.3

实验内容

编程观察无穷级数$$\sum_{n=1}^\infty\frac{1}{n}$$的求和运算.

  • 采用IEEE单精度浮点数,观察当n为何值时求和结果不再变化,将它与理论分析的结论进行比较(注:在MATLAB中可用single命令将变量转成单精度浮点数)。
  • 用IEEE双精度浮点数计算(1)中前n项的和,评估IEEE单精度浮点数计算结果的误差。
  • 如果采用IEEE双精度浮点数,估计当n为何值时求和结果不再变化,这在当前做实验的计算机上大概需要多长的计算时间?

实现思路:

  • 第一问:采用IEEE单精度浮点数, 使用MATLAB的single命令即为IEEE单精度数值,使用single_sumlast_sum记录求和结果,当两个值相等时跳出while循环,保留n值,输出求和结果
  • 不使用single即为Double型数据,对其循环1中的n次,记录求和结果
  • 利用matlab的tic,toc组合命令记录双精度计算n次运行时间,进而推测求值不再变化时的运行时间

结果及分析

程序运行结果如下

1
2
3
4
双精度运算结果: 15.403683 n = 2097152
时间已过 0.011240 秒。
双精度运算结果: 15.133307
单精度运算误差: -0.270376
数据类型位数符号位指数位尾数有效位(10进制)
float3218237
double641115216

1、 根据定理 1.6(即“大数吃小数”定理):$$|\frac{x_2}{x_1}|\le0.5\epsilon, x_1+x_2==x_1$$其中$\epsilon$为机器精度。则当前的 1/n 与当前求和值之比小于等于 $0.5\epsilon$时退出循环。单精度浮点计数时的机器精度为$$2^{-24} \approx 5.96010^{-8},n -> \infty$$时有$$\sum_{n=1}^\infty\frac{1}{n} = ln(n)$$从而当$$\frac{1}{nln(n)} \le 0.5 5.96010^{-8}$$时计算结果不再变化,估算该n值为2000000,这与程序所得结果2097152基本符合.

2、单精度计算前2097152项和结果与双精度计算结果绝对误差为0.270376。

3、单精度浮点计数时的机器精度为$$2^{-53} \approx 1.11010^{-16}$$,有 $$\sum_{n=1}^\infty\frac{1}{n} = ln(n)$$从而当$$\frac{1}{nln(n)} \le 0.5 *1.11010^{-16} = 5.55 * 10^{-17}$$时计算结果不再变化,估算该n值为$n = 500000000000000 = 5.0*10^{14}$, 计算2097152次双精度加法时耗时0.011240秒,则计算n次耗时$$0.011240 * \frac{n}{2097152} \approx 2679824.8291015625 s$$

估算结果

实验心得

初步接触了MATLAB的部分编程方法,对机器精度的理解更加深入了,对单精度和双精度的计算准确度也有了一个直观的了解。

实验二 2.2

实验内容

编程实现阻尼牛顿法。要求:

  • 设定阻尼因子的初始值$\lambda_0$及解的误差阈值$\epsilon$;
  • 阻尼因子$\lambda$用逐次折半法更新;
  • 打印每个迭代步的最终值及近似解。

用所编程序求解:

  1. $x^3 - x - 1 = 0$,取$x_0 = 0.6$
  2. $-x^3 + 5x = 0$, 取$x_0 = 1.2$

实现思路

阻尼牛顿法中迭代新解为

$$x_{k+1} = x_k - \lambda_i \frac{f(x_k)}{f^\prime(x_k)}$$

单调性要求为$$|f(x_{k+1})|<|f(x_k)|, k = 0,1,2,\cdots$$

算法按照课本中算法2.5的伪代码实现

结果及分析

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
========阻尼牛顿法========
初始值x_0:0.6000000000 lambda:1.000000
x_1: 1.1406250000 lambda:0.015625
x_2: 1.3668136616 lambda:1.000000
x_3: 1.3262798040 lambda:1.000000
x_4: 1.3247202256 lambda:1.000000
x_5: 1.3247179572 lambda:1.000000
迭代解为: 1.32471795725 2.04494199352e-11

初始值x_0:1.2000000000 lambda:1.000000
x_1: -1.9411764706 lambda:0.250000
x_2: -2.3204623232 lambda:1.000000
x_3: -2.2404594360 lambda:1.000000
x_4: -2.2360808552 lambda:1.000000
x_5: -2.2360679776 lambda:1.000000
迭代解为: -2.23606797761 1.1124381416e-09

========牛顿法========
初始值x_0:0.6000000000
x_1: 17.9000000000
x_2: 11.9468023286
x_3: 7.9855203519
x_4: 5.3569093148
x_5: 3.6249960329
x_6: 2.5055891901
x_7: 1.8201294223
x_8: 1.4610441099
x_9: 1.3393232243
x_10: 1.3249128677
x_11: 1.3247179926
迭代解为: 1.32471799264 1.50938453736e-07

初始值x_0:1.2000000000
x_1: -5.0823529412
x_2: -3.6219359167
x_3: -2.7660437838
x_4: -2.3576006646
x_5: -2.2448622370
x_6: -2.2361193864
x_7: -2.2360679793
迭代解为: -2.23606797927 1.77279613212e-08

本实验中,$\lambda$均只在第一步迭代中有折半运算,其后都未进入折半循环,减少了总体计算次数,这是一个好结果。

实验三 3.6

实验内容

编程序生成Hilbert矩阵$\mathbf H_n$,以及n维向量$b = \mathbf H_nx$,其中x为所有分量都是1的向量。用Cholesky分解算法求解方程$\mathbf H_nx = b$,得到近似解$\hat x$,计算残差$r = b - H_n\hat x$ 和误差 $\Delta x = \hat x - x$ 的$\infty-$范数。

  1. 设n=10,计算$| r |_\infty、 | \Delta x |_\infty$。
  2. 在右端项上施加$10^{-7}$的扰动然后解方程组,观察残差和误差的变化情况。
  3. 改变n的值为8和12,求解相应方程,观察$| r |_\infty、 | \Delta x |_\infty$的变化情况。通过这个实验说明了什么问题?

实现思路

  1. 根据Hilbert矩阵的结构,生成n阶矩阵$\mathbf H$,继而得到向量$b = \mathbf H_nx$;
  2. 用Cholesky分解算法求解方程$\mathbf H_nx = b$,得到近似解$\hat x$,具体过程为:
    • 首先对$\mathbf H_n$进行Choklesky分解得到下三角矩阵$L_n$
    • $\mathbf H_nx = b, H_n = LL^T$, 则先计算$LY = b$得到Y, 然后计算$L^Tx = Y$,得到$x$
  3. 计算残差和误差的$\infty-$范数。

结果及分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
n =  8 扰动 =  0.0
误差x无穷范数: 6.27198750047e-07
残差r无穷范数: 2.22044604925e-16
n = 8 扰动 = 1e-07
误差x无穷范数: 0.0216211911554
残差r无穷范数: 4.4408920985e-16
n = 10 扰动 = 0.0
误差x无穷范数: 0.000643781555903
残差r无穷范数: 8.881784197e-16
n = 10 扰动 = 1e-07
误差x无穷范数: 0.699727076702
残差r无穷范数: 4.4408920985e-16
n = 12 扰动 = 0.0
误差x无穷范数: 0.355545598229
残差r无穷范数: 4.4408920985e-16
n = 12 扰动 = 1e-07
误差x无穷范数: 23.9972906485
残差r无穷范数: 1.11022302463e-15
  • n = 10, 无扰动时
    $$| r |_\infty = 8.881784197e-16$$
    $$| \Delta x |_\infty = 0.000643781555903$$

  • n = 10, 扰动$10^{-7}$时
    $$| r |_\infty = 4.4408920985e-16$$
    $$| \Delta x |_\infty = 0.699727076702$$

  • n = 8时
    $$| r |_\infty = 2.22044604925e-16$$
    $$| \Delta x |_\infty = 6.27198750047e-07$$

    n = 12时

    $$| r |_\infty = 4.4408920985e-16$$
    $$| \Delta x |_\infty = 0.355545598229$$

进一步实验,统计实验结果

N$| \Delta x |_\infty$引入$10^{-7}$扰动后$| r |_\infty$引入$10^{-7}$扰动后
86.27198750047e-070.02162119115542.22044604925e-164.4408920985e-16
100.0006437815559030.6997270767028.881784197e-164.4408920985e-16
120.35554559822923.99729064854.4408920985e-161.11022302463e-15

由上述结果分析可以得出如下结论:

  • 无论是施加扰动还是改变矩阵的阶数n, 对残差的无穷范数$| r |_\infty$影响不大,但误差的无穷范数$| \Delta x |_\infty$变化明显
  • 矩阵的阶数n增大,误差的无穷范数$| \Delta x |_\infty$增加明显,问题的条件数

\begin{align}
cond &= \frac{|\Delta x|/|x|}{|\Delta b|/|b|} \
&= \frac{|\Delta x|/|x|}{\frac{|\Delta b|}{|H|}/|x|} \
&= \frac{|\Delta x|}{|\Delta b|}|H| \
&= \frac{|\Delta x|}{|r|} |H|
\end{align}

希尔伯特矩阵是一个典型的病态矩阵,其条件数$cond(H_n)_\infty > 1$ 且随着矩阵的阶数n越大,其病态性越严重,因此本题中线性方程组求解问题也是敏感的。

实验四 4.1

实验内容

考虑10阶$Hilbert$矩阵作为系数阵的方程组$$Ax = b$$

其中,A的元素$a_{ij}=\frac{1}{i+j-1}$, $b = \begin{bmatrix}1&\frac{1}{2}&\cdots&\frac{1}{10}\end{bmatrix}^T$.取初始解$x^{(0)} = 0$,编写程序用$Jacobi$与$SOR$迭代法求解该方程组,将$|x^{(k+1)} - x^{(k)}|_\infty < 10^{-4}$作为终止迭代的判据。

  1. 分别用$Jacobi$与$SOR(w = 1.25)$迭代法求解,观察收敛情况;
  2. 改变$w$的值,试验$SOR$迭代法的效果,考察解的准确度。

实现思路

根据教程中的$Jacobi$迭代算法和$SOR$迭代算法的伪代码完成编码

因为A是$Hilbert$矩阵,是实对称矩阵,其谱半径$\rho(A) = |A|_2 > 1$,其$Jacobi$迭代法不收敛。

对于$SOR$迭代法的准确度,我们使用误差$\Delta x = \hat x - x$的无穷范数进行量化,易知线性方程组的正确解$x=\begin{bmatrix}1&0&\cdots&0\end{bmatrix}^T$

结果及分析

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
雅克比迭代法
jacobi迭代不收敛
w = 1.25
迭代次数: 187
迭代结果: [1.002022594409028, -0.021609987149109246, 0.05149585866739306, -0.028416470227697905, -0.006259608843587831, -0.006192548765741887, -0.0005373689012403876, 0.0018347675968743396, 0.003636620501503592, 0.004567285710065993]
w = 1.6
迭代次数: 493
迭代结果: [1.0021469657102071, -0.0283570320963052, 0.08880081215699734, -0.09535773222087951, 0.04336165034147876, -0.0355509984828663, 0.02513571167799267, -0.010543151501989416, 0.013356482959673942, -0.0023146920643427986]
w = 1.4
迭代次数: 270
迭代结果: [1.0022700223737249, -0.027450984607899714, 0.07485459487277796, -0.05746975449051693, 0.005911675449980705, -0.012463346972930323, 0.0036081977048833262, 0.0015499835906313381, 0.00513317184635749, 0.004901511994768535]
w = 1.2
迭代次数: 164
迭代结果: [1.0019214723944498, -0.019188626689019883, 0.04270174478597167, -0.02008673411613389, -0.007120242694937775, -0.005157615844817916, -0.0009805683321413268, 0.0013857829440947498, 0.0029820023152041636, 0.0039040390400999197]
w = 1.0
迭代次数: 2
迭代结果: [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
w = 0.8
迭代次数: 178
迭代结果: [0.9973420757985196, 0.022189810579507762, -0.03692216721730894, 0.0002172210934328124, 0.012508220958753955, 0.010340953138668948, 0.0047473091369648925, -0.00029255227286110113, -0.003947813463373296, -0.006348095947258598]
w = 0.6
迭代次数: 271
迭代结果: [0.9955833379103364, 0.038088293718681644, -0.06191668688342534, -0.008191119893787442, 0.022657640869212896, 0.024947505162692972, 0.014412072529676016, 0.0013699258800427523, -0.009816061664535391, -0.018012206700100845]
w = 0.4
迭代次数: 397
迭代结果: [0.9932920658677347, 0.05517010926044667, -0.08070115313478526, -0.023914639520845007, 0.024971861498459945, 0.03919911225602069, 0.029505998466182924, 0.008987379128371036, -0.013671954578940389, -0.03395765285066289]
w = 0.2
迭代次数: 657
迭代结果: [0.9886002689546933, 0.081327354825667, -0.09358751789761771, -0.048223519275488166, 0.012311655912579827, 0.04317339038553021, 0.043997018466545454, 0.024122149074273676, -0.007160625602977868, -0.0428600869164585]
w$| \Delta x |_\infty$迭代次数
1.60.0953577322209493
1.40.0748545948728270
1.250.0514958586674187
1.20.042701744786164
1.00.02
0.80.0369221672173178
0.60.0619166868834271
0.40.0807011531348397
0.20.0935875178976657
  1. $Jacobi$迭代法不收敛,$SOR$迭代法($w = 1.25$)迭代187次后收敛,得到近似解$\hat x$ = [1.002022594409028, -0.021609987149109246, 0.05149585866739306, -0.028416470227697905, -0.006259608843587831, -0.006192548765741887, -0.0005373689012403876, 0.0018347675968743396, 0.003636620501503592, 0.004567285710065993], 误差的无穷范数为$| \Delta x |_\infty = 0.0514958586674$

  2. 选取不同的w值,得到的迭代次数和误差的无穷范数变化如上表所示,可以得到的结论是w值越接近1,其收敛速度越快,误差越小。

实验心得

对$Jacobi$迭代法和$SOR$迭代法的收敛条件有了更深入的认识,对$SOR$迭代法中松弛因子$w$的大小对迭代算法收敛速度及计算准确度的影响有了直观的感受,针对一类应用问题选择合适的松弛因子是一个重要问题。

实验五 5.1

实验内容

用幂法求下列矩阵按模最大特征值$\lambda_1$及其对应的特征向量$x_1$,使$|(\lambda_1)_{k+1} - (\lambda_1)_k| < 10^{-5}$。

  1. $A = \begin{bmatrix}5&-4&1\-4&6&-4\1&-4&7\end{bmatrix}$
  2. $B = \begin{bmatrix}25&-41&10&-6\-41&68&-17&10\10&-17&5&-3\-6&10&-3&2\end{bmatrix}$

实现思路

按照课本中算法5.1实现算法

结果及分析

  1. 迭代次数 16

    主特征向量: [ 0.67401996 -1. 0.8895594 ]

    主特征值: 12.2543197032

  2. 迭代次数 6

    主特征向量: [-0.60397234 1. -0.25113513 0.14895345]

    主特征值: 98.5216977228

实验六 6.3

实验内容

对物理实验中所得的下列数据

$t_i$11.522.53.03.544.5
$y_i$33.4079.50122.65159.05189.15214.15238.65252.2
$t_i$5.05.56.06.57.07.58.0
$y_i$267.55280.50296.65301.65310.40318.15325.15
  1. 用公式$y = a + bt + ct^2$做曲线拟合.
  2. 用指数函数$y = ae^{bt}$做曲线拟合.
  3. 比较上述两条拟合曲线,哪条更好?

实现思路

  1. 参考课本算法6.2(用法方程方法求解曲线拟合的最小二乘问题)。
  2. 首先根据式(6.28)形成矩阵A,继而计算出$G = A^TA$和向量$b = A^Tf$,得到方程$Gx = b$。
  3. 解上述所得方程。可以采用$Cholesky$分解然后执行前代和回代过程,得到的向量即为多项式函数的各项系数(低次到高次)。
  4. 用指数函数拟合时,可以对函数式两边同时取对数,转化为一次多项式函数的曲线拟合问题。

结果及分析

  1. 拟合结果$y=-45.29423077+94.19429218t-6.12682612t^2$, 均方误差: 5.68393182348
  2. 拟合结果$y=67.39379285e^{0.23898344t}$, 均方误差: 63.1840873401

通过比较其均方误差可知:用2次多项式函数拟合曲线比指数函数拟合曲线好,其拟合图像如下所示

拟合图像

打赏