LU分解和恰定方程组的解

作者：未知信息来源:未知 2006-1-26

字体大小:小 中 大网友评论条 进入论坛

前言接下来将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程求解等。然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，拟合和插值， Fourier 分析，和一般常微分方程初值问题。2 LU分解和恰定方程组的解5。2 恰定方程组的解【 * 例 5。...

前言

接下来将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程求解等。但与一般数值计算教科书不同，本章的讨论重点：如何利用现有的世界顶级数值计算资源 MATLAB 。至于数学描述，本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和 MATLAB 计算指令之间的内在联系及区别。

对于那些熟悉其他高级语言（如 FORTRAN ， Pascal ， C++ ）的读者来说，通过本章， MATLAB 卓越的数组处理能力、浩瀚而灵活的 M 函数指令、丰富而友善的图形显示指令将使他们体验到解题视野的豁然开朗，感受到摆脱烦琐编程后的眉眼舒展。

对于那些经过大学基本数学教程的读者来说，通过本章， MATLAB 精良完善的计算指令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。

对于那些熟悉 MATLAB 基本指令的读者来说，通过本章，围绕基本数值问题展开的内容将使他们体会到各别指令的运用场合和内在关系，获得综合运用不同指令解决具体问题的思路和借鉴。

由于 MATLAB 的基本运算单元是数组，所以本章内容将从矩阵分析、线性代数的数值计算开始。然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，拟合和插值， Fourier 分析，和一般常微分方程初值问题。本章的最后讨论稀疏矩阵的处理，因为只有在大型问题中，才须特别处理。

从总体上讲，本章各节之间没有依从关系，即读者没有必要从头到尾系统阅读本章内容。读者完全可以根据需要阅读有关节次。除特别说明外，每节中的例题指令是独立完整的，因此读者可以很容易地在自己机器上实践。

5.2 LU分解和恰定方程组的解

5.2.2 恰定方程组的解

【 * 例 5.2.2 -1 】“求逆”法和“左除”法解恰定方程的性能对比

（1）为对比这两种方法的性能，先用以下指令构造一个条件数很大的高阶恰定方程。
rand('state',12); % 选定随机种子，目的是可重复产生随机阵 A 。
A=rand(100,100)+1.e8; %rand(100,100) 生成（ 100 × 100 ）均匀分布随机矩阵。
% 每个随机阵元素加的目的是使 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 Timer
xi=inv(A)*b; % xi 是用“求逆”法解恰定方程所得的解。
ti=toc % 关闭计时器，并显示解方程所用的时间。
ci=flops % “求逆”法解方程所用的运算次数
eri=norm(x-xi) % 解向量 xi 与真解向量 x 的范 -2 误差。
rei=norm(A*xi-b)/norm(b) % 方程的范 -2 相对残差
ti =
0.9300
ci =
2070322
eri =
3.0708e-004
rei =
6.6280e-007

（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.2.3 范数、条件数和方程解的精度

【 * 例 5.2.3 -1 】 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=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('')
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