《基于matlab的数值计算》由会员分享，可在线阅读，更多相关《基于matlab的数值计算（43页珍藏版）》请在金锄头文库上搜索。
1、 1 第四章 第四章 数值计算数值计算 4.1 引言引言 本章将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微 积分、数据分析、以及常微分方程（初值和边值问题）求解等。但与一般数值计算教科书不 同，本章的讨论重点是：如何利用现有的世界顶级数值计算资源 MATLAB。至于数学描述， 本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和 MATLAB 计算指令之间的内在联系及区别。 对于那些熟悉其他高级语言（如 FORTRAN，Pascal，C+）的读者来说，通过本章， MATLAB 卓越的数组处理能力、浩瀚而灵活的 M 函数指令、丰富而友善的图形
2、显示指令将 使他们体验到解题视野的豁然开朗，感受到摆脱烦琐编程后的眉眼舒展。 对于那些经过大学基本数学教程的读者来说，通过本章，MATLAB 精良完善的计算指 令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简 单算例通向科学研究和工程设计实际问题的一条途径。 对于那些熟悉 MATLAB 基本指令的读者来说，通过本章，围绕基本数值问题展开的内 容将使他们体会到各别指令的运用场合和内在关系， 获得综合运用不同指令解决具体问题的 思路和借鉴。 由于 MATLAB 的基本运算单元是数组，所以本章内容将从矩阵分析、线性代数的数值 计算开始。然后再介绍函数零点、极值的求取，
3、数值微积分，数理统计和分析，拟合和插值， Fourier 分析，和一般常微分方程初值、边值问题。本章的最后讨论稀疏矩阵的处理，因为 这只有在大型问题中，才须特别处理。 从总体上讲， 本章各节之间没有依从关系， 即读者没有必要从头到尾系统阅读本章内容。 读者完全可以根据需要阅读有关节次。除特别说明外，每节中的例题指令是独立完整的，因 此读者可以很容易地在自己机器上实践。 MATLAB 从 5.3 版升级到 6.x 版后，本章内容的变化如下： ? MATLAB 从 6.0 版起，其矩阵和特征值计算指令不再以 LINPACK和 EISPACK库为基 础，而建筑在计算速度更快、运行更可靠的 LAPAC
4、K和 ARPACK程序库的新基础上。 因此，虽然各种矩阵计算指令没有变化，但计算结果却可能有某些不同。这尤其突出地 表现在涉及矩阵分解、特征向量、奇异向量等的计算结果上。对此，用户不必诧异，因 为构成空间的基向量时不唯一的，且新版的更可信。本书新版全部算例结果是在 6.x 版 上给出的。 ? 在 5.3 版本中，泛函指令对被处理函数的调用是借助函数名字符串进行的。这种调用方 式在 6.x 版中已被宣布为“过渡期内允许使用但即将被淘汰的调用方式”；而新的调用 方式是借助“函数句柄”进行的。因此，关于述泛函指令，本章新版着重讲述如何使用 “函数句柄”，同时兼顾“函数名字符串”调用法。 ? MATL
5、AB 从 6.0 版起，提供了一组专门求微分方程“边值问题”数值解的指令。适应 这种变化，本章新增第 4.14.5 节，用 2 个算例阐述求解细节。 ? 5.3 版中的积分指令 quad8 已经废止；6.x 版启用新积分指令 quadl ；6.5 版新增三重积分 指令 triplequad。本章新版对此作了相应的改变。 4.2 LU 分解和恰定方程组的解分解和恰定方程组的解 2 4.2.1 LU 分解、行列式和逆分解、行列式和逆 4.2.2 恰定方程组的解恰定方程组的解 【例 4.2.2-1】“求逆”法和“左除”法解恰定方程的性能对比 （1） randn(state,0); A=gallery
6、(randsvd,100,2e13,2); x=ones(100,1); b=A*x; cond(A) ans = 1.9990e+013 （2） tic %启动计时器启动计时器 stopwatch timer xi=inv(A)*b; %xi 用求逆法解恰定方程所得的解用求逆法解恰定方程所得的解 ti=toc %关闭计时器，并显示方程所用的时间关闭计时器，并显示方程所用的时间 eri=norm(x-xi) %解向量解向量 xi 与真解向量与真解向量 x 的误差范数的误差范数 rei=norm(A*xi-b)/norm(b) %方程范数的相对残差方程范数的相对残差 ti = 0.7700 er
7、i = 0.0469 rei = 0.0047 （3） tic;xd=Ab; 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 矩阵及其准确逆 矩阵的生成函数。本例将对方程bHx =近似解和准确解进行比较。 N=6 8 10 12 14; for k=1:length(N) n=N(k); H=hilb(n);
8、 Hi=invhilb(n); b=ones(n,1); x_approx=Hb; 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 3 disp(Hilbert 矩阵阶数矩阵阶数),disp(N) format short e disp(实际误差实际误差 er_actual),disp(er_actu
9、al),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 最大可能
10、误差 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】本例演示：如矩阵中有元素与截断误差相当时的特性值
11、问题。 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.1
12、055 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(d
13、iag(D); 4 d,D(1,1) T_full = 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)
14、;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:MATLAB6P1toolboxmatlabmatfuncondeig.m at line 30 c_eig = 1.0000e+000 1.0000e+
15、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 5 -0.0408 0 0 0 0 -0.0119 + 0.0386i 0 0 0 0 -0.0119 -
