第四章 数值计算
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版启用新积分指令quadl ；6.5版新增三重积
分指令triplequad。本章新版对此作了相应的改变。
4.2 LU分解和恰定方程组的解
1
4.2.1 LU分解、行列式和逆 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 = 0 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
2
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 0.0194 -0.0222 0.0222 0.0333
【例4.3.1-3】指令eig与eigs的比较。
rand('state',1),A=rand(100,100)-0.5;
t0=clock;[V,D]=eig(A);T_full=etime(clock,t0) options.tol=1e-8; options.disp=0;
t0=clock;[v,d]=eigs(A,1,'lr',options); T_part=etime(clock,t0)
[Dmr,k]=max(real(diag(D))); d,D(1,1) T_full =
3
0.2200 T_part = 3.1300 d =
3.0140 + 0.2555i ans =
3.0140 + 0.2555i
vk1=V(:,k);
vk1=vk1/norm(vk1);v=v/norm(v); V_err=acos(norm(vk1'*v))*180/pi D_err=abs(D(k,k)-d)/abs(d) V_err =
1.2074e-006 D_err =
4.2324e-010
4.3.2 特征值问题的条件数
【例4.3.2-1】矩阵的代数方程条件数和特征值条件数。
B=eye(4,4);B(3,4)=1;B
format short e,c_equ=cond(B),c_eig=condeig(B) B =
1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1 c_equ =
2.6180e+000
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.110223e-016. > In D:\\MATLAB6P1\\toolbox\\matlab\\matfun\\condeig.m at line 30 c_eig =
1.0000e+000 1.0000e+000 4.5036e+015 4.5036e+015
【例4.3.2-2】对亏损矩阵进行Jordan分解。
A=gallery(5)
[VJ,DJ]=jordan(A);
[V,D,c_eig]=condeig(A);c_equ=cond(A); DJ,D,c_eig,c_equ A =
-9 11 -21 63 -252 70 -69 141 -421 1684 -575 575 -1149 3451 -13801 3891 -3891 7782 -23345 93365 1024 -1024 2048 -6144 24572 DJ =
0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 D =
Columns 1 through 4
-0.0408 0 0 0 0 -0.0119 + 0.0386i 0 0 0 0 -0.0119 - 0.0386i 0 0 0 0 0.0323 + 0.0230i
4
0 0 0 0 Column 5
0 0 0 0 0.0323 - 0.0230i c_eig =
1.0e+010 * 2.1293 2.0796 2.0796 2.0020 2.0020 c_equ =
2.0253e+018
4.3.3 复数特征值对角阵与实数块特征值对角阵的转化
【例4.3.3-1】把例4.3.1-1中的复数特征值对角阵D转换成实数块对角阵，使VR*DR/VR=A。
A=[1,-3;2,2/3];[V,D]=eig(A); [VR,DR]=cdf2rdf(V,D) VR =
0.7746 0 0.0430 -0.6310 DR =
0.8333 2.4438 -2.4438 0.8333
4.3.4 矩阵的谱分解和矩阵函数
【例4.3.4-1】数组乘方与矩阵乘方的比较。
clear,A=[1 2 3;4 5 6;7 8 9]; A_Ap=A.^0.3 A_Mp=A^0.3 A_Ap =
1.0000 1.2311 1.3904 1.5157 1.6207 1.7118 1.7928 1.8661 1.9332 A_Mp =
0.6962 + 0.6032i 0.4358 + 0.1636i 0.1755 - 0.2759i 0.6325 + 0.0666i 0.7309 + 0.0181i 0.8292 - 0.0305i 0.5688 - 0.4700i 1.0259 - 0.1275i 1.4830 + 0.2150i
【例4.3.4- 2】标量的数组乘方和矩阵乘方的比较。（A取自例4.3.4-1）
pA_A=(0.3).^A pA_M=(0.3)^A pA_A =
0.3000 0.0900 0.0270 0.0081 0.0024 0.0007 0.0002 0.0001 0.0000 pA_M =
2.9342 0.4175 -1.0993 -0.0278 0.7495 -0.4731 -1.9898 -0.9184 1.1531
【例4.3.4-3】sin的数组运算和矩阵运算比较。（A取自例4.3.4-1）
A_sinA=sin(A)
A_sinM=funm(A,'sin')
5