数值计算方法
5.1引言
本章将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程求解等。但与一般数值计算教科书不同，本章的讨论重点是：如何利用现有的世界顶级数值计算资源MATLAB。至于数学描述，本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和MATLAB计算指令之间的内在联系及区别。 对于那些熟悉其他高级语言（如FORTRAN，Pascal，C++）的读者来说，通过本章，MATLAB卓越的数组处理能力、浩瀚而灵活的M函数指令、丰富而友善的图形显示指令将使他们体验到解题视野的豁然开朗，感受到摆脱烦琐编程后的眉眼舒展。 对于那些经过大学基本数学教程的读者来说，通过本章，MATLAB精良完善的计算指令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。 对于那些熟悉MATLAB基本指令的读者来说，通过本章，围绕基本数值问题展开的内容将使他们体会到各别指令的运用场合和内在关系，获得综合运用不同指令解决具体问题的思路和借鉴。 由于MATLAB的基本运算单元是数组，所以本章内容将从矩阵分析、线性代数的数值计算开始。然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，拟合和插值，Fourier分析，和一般常微分方程初值问题。本章的最后讨论稀疏矩阵的处理，因为这只有在大型问题中，才须特别处理。 从总体上讲，本章各节之间没有依从关系，即读者没有必要从头到尾系统阅读本章内容。读者完全可以根据需要阅读有关节次。除特别说明外，每节中的例题指令是独立完整的，因此读者可以很容易地在自己机器上实践。
5.1 LU分解和恰定方程组的解
5.1.1 LU分解、行列式和逆
（1）LU分解 （2）行列式和逆
5.1.2 恰定方程组的解
【*例5.2.2-1】“求逆”法和“左除”法解恰定方程的性能对比
（1）为对比这两种方法的性能，先用以下指令构造一个条件数很大的高阶恰定方程。
rand('state',12); %选定随机种子，目的是可重复产生随机阵A。
A=rand(100,100)+1.e8; %rand(100,100)生成（100×100）均匀分布随机矩阵。 x=ones(100,1); b=A*x; cond(A) ans =
1.4426e+012
flops(0);tic
%每个随机阵元素加10的目的是使A阵条件数升高。 %令解向量 x 为全1的100元列向量。
%为使 Ax=b 方程一致，用A和 x 生成 b 向量。 %求A阵的条件数。
8（2）“求逆”法解恰定方程的误差、残差、运算次数和所用时间
%浮点运算计数器置0 ；启动计时器Stopwatch Timer
1
xi=inv(A)*b; ti=toc ci=flops eri=norm(x-xi) rei=norm(A*xi-b)/norm(b) ti =
0.9300 ci =
2070322 eri =
3.0708e-004 rei =
6.6280e-007
% xi 是用“求逆”法解恰定方程所得的解。 %关闭计时器，并显示解方程所用的时间。 %“求逆”法解方程所用的运算次数
%解向量 xi 与真解向量 x 的范-2误差。 %方程的范-2相对残差
（3）“左除”法解恰定方程的误差、残差、运算次数和所用时间
flops(0);tic;xd=A\\b; %是用“左除”法解恰定方程所得的解。 td=toc,cd=flops,erd=norm(x-xd),red=norm(A*xd-b)/norm(b) td =
0.2200 cd =
741872 erd =
3.2243e-004 red =
2.0095e-016
5.1.3 范数、条件数和方程解的精度
【*例5.2.3-1】Hilbert矩阵是著名的病态矩阵。MATLAB中有专门的Hilbert矩阵及其准确逆矩阵的生成函数。本例将对方程Hx?b近似解和准确解进行比较。所谓n阶Hilbert矩阵
?1??i?j的形式是：Hn??hij?n?n,hij??1???i?j?1i?j 。 i?jN=[6 8 10 12 14]; %本例计算的矩阵阶数 for k=1:length(N) n=N(k); %矩阵的阶 H=hilb(n); %产生n阶Hilbert矩阵 Hi=invhilb(n); %产生完全准确的n阶逆Hilbert矩阵 b=ones(n,1); %生成n阶全1向量 x_approx=H\\b; %利用左除H求近似解 x_exact=Hi*b; %利用准确逆Hilbert矩阵求准确解 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); %计算Hilbert矩阵的条件数 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('')
2
Hilbert矩阵阶数
6 8 10 12 14 实际误差 er_actual
5.0339e-011 8.5981e-008 2.2819e-004 1.3381e-001 3.9641e+000 近似的最大可能误差 er_approx
3.3198e-009 3.3879e-006 3.5583e-003 3.9259e+000 3.4573e+002 最大可能误差 er_max
6.0095e-007 2.4531e-002 1.4094e+003 2.9206e+007 2.4178e+010
5.2 矩阵特征值和矩阵函数
5.2.1 特征值和特征向量的求取
【例5.3.1-1】简单实阵的特征值问题。
A=[1,-3;2,2/3];[V,D]=eig(A) V =
-0.7728 + 0.0527i -0.7728 - 0.0527i 0 + 0.6325i 0 - 0.6325i D =
0.8333 + 2.4438i 0 0 0.8333 - 2.4438i
【*例5.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 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 1.1227 ER2 =
1.0e-015 *
0.4441 -0.2220 0.1471 -0.2220 0 0.0555 -0.3629 0.2776 -0.0172 -0.0015 0.0066 0 0 -0.2220 -0.1110 0.1388
【*例5.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)%指令eig的运作时间。 options.tol=1e-8; %为eigs设定计算精度。 options.disp=0; %使中间迭代结果不显示。
t0=clock;[v,d]=eigs(A,1,'lr',options);%计算最大实部特征值和特征向量。 T_part=etime(clock,t0) %指令eigs的运作时间。
[Dmr,k]=max(real(diag(D))); %在eig求得的全部特征值中找最大实部的那个。 d,D(1,1) T_full = 2.8000 T_part = 19.5500 d =
3.0140 - 0.2555i ans =
3.0140 + 0.2555i vk1=V(:,k+1); %与d相同的特征向量应是V的第k+1列。 vk1=vk1/norm(vk1);v=v/norm(v); %向量长度归一。
3
V_err=acos(norm(vk1'*v))*180/pi %求复数向量之间的夹角（度）。 D_err=abs(D(k+1,k+1)-d)/abs(d) %求两个特征值间的相对误差。 V_err =
8.5377e-007 D_err =
1. 7098e-010
5.2.2 特征值问题的条件数
【例5.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.500000e-018. > In E:\\MAT53\\toolbox\\matlab\\matfun\\condeig.m at line 30 c_eig =
1.0000e+000 1.0000e+000 3.3333e+017 3.3333e+017
【*例5.3.2-2】对亏损矩阵进行Jordan分解。
A=gallery(5) %MATLAB设置的特殊矩阵，它具有五重特征值。 [VJ,DJ]=jordan(A); %求出准确的特征值，使A*VJ=VJ*D成立。 [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.0328 + 0.0243i 0 0 0 0 -0.0328 - 0.0243i 0 0 0 0 0.0130 + 0.0379i 0 0 0 0 0.0130 - 0.0379i
0 0 0 0 Column 5
0 0 0 0 0.0396
c_eig =
1.0e+010 * 2.1016
4
2.1016 2.0251 2.0251 1.9796 c_equ =
5.2133e+017
5.2.3 复数特征值对角阵与实数块特征值对角阵的转化
【*例5.3.3-1】把例5.3.1-1中的复数特征值对角阵D转换成实数块对角阵，使VR*DR/VR=A。
[VR,DR]=cdf2rdf(V,D)
VR =
-0.7728 0.0527 0 0.6325 DR =
0.8333 2.4438 -2.4438 0.8333
5.2.4 矩阵的谱分解和矩阵函数
【*例5.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
【*例5.3.4- 2】标量的数组乘方和矩阵乘方的比较。（A取自例5.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
【*例5.3.4-3】sin的数组运算和矩阵运算比较。（A取自例5.3.4-1）
A_sinA=sin(A) %数组运算 A_sinM=funm(A,'sin') %矩阵运算 A_sinA =
0.8415 0.9093 0.1411 -0.7568 -0.9589 -0.2794 0.6570 0.9894 0.4121 A_sinM =
-0.6928 -0.2306 0.2316 -0.1724 -0.1434 -0.1143
0.3479 -0.0561 -0.4602
5.3 奇异值分解
5
数值计算方法精讲
