第四章数值计算
4.1引言
本章将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程（初值和边值问题）求解等。但与一般数值计算教科书不同，本章的讨论重点是：如何利用现有的世界顶级数值计算资源MATLAB。至于数学描述，本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和MATLAB计算指令之间的在联系及区别。
对于那些熟悉其他高级语言（如FORTRAN，Pascal，C++）的读者来说，通过本章，MATLAB 卓越的数组处理能力、浩瀚而灵活的M函数指令、丰富而友善的图形显示指令将使他们体验到解题视野的豁然开朗，感受到摆脱烦琐编程后的眉眼舒展。
对于那些经过大学基本数学教程的读者来说，通过本章，MATLAB精良完善的计算指令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。
对于那些熟悉MATLAB基本指令的读者来说，通过本章，围绕基本数值问题展开的容将使他们体会到各别指令的运用场合和在关系，获得综合运用不同指令解决具体问题的思路和借鉴。
由于MATLAB的基本运算单元是数组，所以本章容将从矩阵分析、线性代数的数值计算开始。然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，拟合和插值，Fourier分析，和一般常微分方程初值、边值问题。本章的最后讨论稀疏矩阵的处理，因为这只有在大型问题中，才须特别处理。
从总体上讲，本章各节之间没有依从关系，即读者没有必要从头到尾系统阅读本章容。读者完全可以根据需要阅读有关节次。除特别说明外，每节中的例题指令是独立完整的，因此读者可以很容易地在自己机器上实践。
MATLAB从5.3版升级到6.x版后，本章容的变化如下：
●MATLAB从6.0版起，其矩阵和特征值计算指令不再以LINPACK和EISPACK库为基础，
而建筑在计算速度更快、运行更可靠的LAPACK和ARPACK程序库的新基础上。因此，虽然各种矩阵计算指令没有变化，但计算结果却可能有某些不同。这尤其突出地表现在涉及矩阵分解、特征向量、奇异向量等的计算结果上。对此，用户不必诧异，因为构成空间的基向量时不唯一的，且新版的更可信。本书新版全部算例结果是在6.x版上给出的。
●在5.3版本中，泛函指令对被处理函数的调用是借助函数名字符串进行的。这种调用
方式在6.x版中已被宣布为“过渡期允许使用但即将被淘汰的调用方式”；而新的调用方式是借助“函数句柄”进行的。因此，关于述泛函指令，本章新版着重讲述如何使用“函数句柄”，同时兼顾“函数名字符串”调用法。
●MATLAB从6.0版起，提供了一组专门求微分方程“边值问题”数值解的指令。适应这
种变化，本章新增第4.14.5节，用2个算例阐述求解细节。
● 5.3版中的积分指令quad8已经废止；6.x版启用新积分指令quad l；6.5版新增三重
积分指令triplequad。本章新版对此作了相应的改变。
4.2LU分解和恰定方程组的解
4.2.1LU分解、行列式和逆
4.2.2恰定方程组的解
【例4.2.2-1】“求逆”法和“左除”法解恰定方程的性能对比
（1）
randn('state',0);
A=gallery('randsvd',100,2e13,2);
x=ones(100,1);
b=A*x;
cond(A)
ans =
1.9990e+013
（2）
tic
xi=inv(A)*b;
ti=toc
eri=norm(x-xi)
rei=norm(A*xi-b)/norm(b)
ti =
0.7700
eri =
0.0469
rei =
0.0047
（3）
tic;xd=A\b;
td=toc,erd=norm(x-xd),red=norm(A*xd-b)/norm(b)
td =
erd =
0.0078
red =
2.6829e-015
4.2.3数、条件数和方程解的精度
【例4.2.3-1】Hilbert矩阵是著名的病态矩阵。MATLAB中有专门的Hilbert矩阵及其准确
Hx 近似解和准确解进行比较。
逆矩阵的生成函数。本例将对方程b
N=[6 8 10 12 14];
for k=1:length(N)
n=N(k);
H=hilb(n);
Hi=invhilb(n);
b=ones(n,1);
x_approx=H\b;
x_exact=Hi*b;
ndb=norm(H*x_approx-b);nb=norm(b);
ndx=norm(x_approx - x_exact);nx=norm(x_approx);
er_actual(k)=ndx/nx;
K=cond(H);
er_approx(k)=K*eps;
er_max(k)=K*ndb/nb;
end
disp('Hilbert矩阵阶数'),disp(N)
format short e
disp('实际误差 er_actual'),disp(er_actual),disp('')
disp('近似的最大可能误差 er_approx'),disp(er_approx),disp('')
disp('最大可能误差 er_max'),disp(er_max),disp('')
Hilbert矩阵阶数
6 8 10 12 14
实际误差 er_actual
1.5410e-010 1.7310e-007 1.9489e-004 9.1251e-002
2.1257e+000
近似的最大可能误差 er_approx
3.3198e-009 3.3879e-006 3.5583e-003 3.9846e+000 9.0475e+001
最大可能误差 er_max
7.9498e-007 3.8709e-002 1.2703e+003 4.7791e+007 4.0622e+010
4.3矩阵特征值和矩阵函数
4.3.1特征值和特征向量的求取
【例4.3.1-1】简单实阵的特征值问题。
A=[1,-3;2,2/3];[V,D]=eig(A)
V =
0.7746 0.7746
0.0430 - 0.6310i 0.0430 + 0.6310i
D =
0.8333 + 2.4438i 0
0 0.8333 - 2.4438i
【例4.3.1-2】本例演示：如矩阵中有元素与截断误差相当时的特性值问题。A=[3 -2 -0.9 2*eps
-2 4 -1 -eps
-eps/4 eps/2 -1 0
-0.5 -0.5 0.1 1 ];
[V1,D1]=eig(A);ER1=A*V1-V1*D1
[V2,D2]=eig(A,'nobalance');ER2=A*V2-V2*D2
ER1 =
0.0000 0.0000 0.0000 0.0000
0 -0.0000 -0.0000 -0.0000
0.0000 -0.0000 -0.0000 0.0000
0.0000 0.0000 0.0000 -0.5216
ER2 =
1.0e-014 *
-0.2665 0.0111 -0.0559 -0.1055
0.4441 0.1221 0.0343 0.0833
0.0022 0.0002 0.0007 0
