《MATLAB教案.ch05.数值计算》由会员分享，可在线阅读，更多相关《MATLAB教案.ch05.数值计算（50页珍藏版）》请在金锄头文库上搜索。
1、5 数值计算5.1 引言本章将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程求解等。但与一般数值计算教科书不同，本章的讨论重点是：如何利用现有的世界顶级数值计算资源MATLAB。至于数学描述，本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和MATLAB计算指令之间的内在联系及区别。对于那些熟悉其他高级语言（如FORTRAN，Pascal，C+）的读者来说，通过本章，MATLAB卓越的数组处理能力、浩瀚而灵活的M函数指令、丰富而友善的图形显示指令将使他们体验到解题视野的豁然开朗，感受到摆脱烦琐编程后的眉眼舒展。对
2、于那些经过大学基本数学教程的读者来说，通过本章，MATLAB精良完善的计算指令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。对于那些熟悉MATLAB基本指令的读者来说，通过本章，围绕基本数值问题展开的内容将使他们体会到各别指令的运用场合和内在关系，获得综合运用不同指令解决具体问题的思路和借鉴。由于MATLAB的基本运算单元是数组，所以本章内容将从矩阵分析、线性代数的数值计算开始。然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，拟合和插值，Fourier分析，和一般常微分方程初值问题。本章的最后讨论稀疏矩
3、阵的处理，因为这只有在大型问题中，才须特别处理。从总体上讲，本章各节之间没有依从关系，即读者没有必要从头到尾系统阅读本章内容。读者完全可以根据需要阅读有关节次。除特别说明外，每节中的例题指令是独立完整的，因此读者可以很容易地在自己机器上实践。5.2 LU分解和恰定方程组的解5.2.1 LU分解、行列式和逆（1）LU分解（2）行列式和逆5.2.2 恰定方程组的解【*例】“求逆”法和“左除”法解恰定方程的性能对比（1）为对比这两种方法的性能，先用以下指令构造一个条件数很大的高阶恰定方程。rand(state,12);%选定随机种子，目的是可重复产生随机阵A。A=rand(100,100)+1.e8
4、;%rand(100,100)生成（100100）均匀分布随机矩阵。%每个随机阵元素加的目的是使A阵条件数升高。x=ones(100,1);%令解向量 x 为全1的100元列向量。b=A*x;%为使 Ax=b 方程一致，用A和 x 生成 b 向量。cond(A) %求A阵的条件数。 ans = 1.4426e+012 （2）“求逆”法解恰定方程的误差、残差、运算次数和所用时间flops(0);tic%浮点运算计数器置0 ；启动计时器Stopwatch Timerxi=inv(A)*b;% xi 是用“求逆”法解恰定方程所得的解。ti=toc%关闭计时器，并显示解方程所用的时间。ci=flops
5、%“求逆”法解方程所用的运算次数eri=norm(x-xi)%解向量 xi 与真解向量 x 的范-2误差。rei=norm(A*xi-b)/norm(b)%方程的范-2相对残差 ti = 0.9300ci = 2070322eri = 3.0708e-004rei = 6.6280e-007 （3）“左除”法解恰定方程的误差、残差、运算次数和所用时间flops(0);tic;xd=Ab;%是用“左除”法解恰定方程所得的解。td=toc,cd=flops,erd=norm(x-xd),red=norm(A*xd-b)/norm(b) td = 0.2200cd = 741872erd = 3.2
6、243e-004red = 2.0095e-016 5.2.3 范数、条件数和方程解的精度【*例】Hilbert矩阵是著名的病态矩阵。MATLAB中有专门的Hilbert矩阵及其准确逆矩阵的生成函数。本例将对方程近似解和准确解进行比较。所谓n阶Hilbert矩阵的形式是： 。N=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=Hb; %利用左除H求近似解x_exac
7、t=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 edisp(实际误差 er_actual),disp(er_actual)
8、,disp()disp(近似的最大可能误差 er_approx),disp(er_approx),disp()disp(最大可能误差 er_max),disp(er_max),disp() 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.
9、4531e-002 1.4094e+003 2.9206e+007 2.4178e+010 5.3 矩阵特征值和矩阵函数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.6325iD = 0.8333 + 2.4438i 0 0 0.8333 - 2.4438i 【*2】本例演示：如矩阵中有元素与截断误差相当时的特性值问题。A=3 -2 -0.9 2*eps -2 4 -1 -eps -eps/4 eps/2 -1 0 -
10、0.5 -0.5 0.1 1 ;V1,D1=eig(A);ER1=A*V1-V1*D1V2,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.1227ER2 = 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.2
11、220 -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);%在e
12、ig求得的全部特征值中找最大实部的那个。d,D(1,1) T_full = 2.8000T_part = 19.5500d = 3.0140 - 0.2555ians = 3.0140 + 0.2555i vk1=V(:,k+1);%与d相同的特征向量应是V的第k+1列。vk1=vk1/norm(vk1);v=v/norm(v);%向量长度归一。V_err=acos(norm(vk1*v)*180/pi%求复数向量之间的夹角（度）。D_err=abs(D(k+1,k+1)-d)/abs(d)%求两个特征值间的相对误差。 V_err = 8.5377e-007D_err =1. 7098e-01
13、0 5.3.2 特征值问题的条件数【例1】矩阵的代数方程条件数和特征值条件数。B=eye(4,4);B(3,4)=1;Bformat 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 1c_equ = 2.6180e+000Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.500000e-018. In E:MAT53toolboxmatlabmatfuncondeig.m at line 30c_eig = 1.0000e+000 1.0000e+000 3.3333e+017 3.3333e+017 【*例】对亏损矩阵进行Jordan分解。A=gallery(5)
