模拟物理PPT课件

1、引言 在描述物理过程时，线性化是一个常见的假设和近似。 因此线性方程组在计算物理中是到处都会遇到的。 与求本征值和解线性方程组相联系的矩阵计算常常是许多计算物理问题中计算工作量的一大部分。第1页/共56页稠密与稀疏 本章讨论几种“直接”方法 它们适用于维数不超过几百的“稠密”矩阵（大部分元素不是零）。 稀疏矩阵中大量的元素等于零。它们出现在对常微分方程和偏微分方程进行离散化的情况下。 随后两章讨论稀疏矩阵的迭代方法。第2页/共56页 同特殊函数的情况一样，处理矩阵的子程序已经包含在一些程序库中。在大型的计算机上都是现成的。 但是，对于基本概念有一些了解是必要的。第3页/共56页5.1 矩阵求逆

2、 我们考虑一个 的方阵A求逆的问题。 例如在求解线性方程组的过程需要求逆 上式中b是已知矢量，x是未知矢量。 解决这个问题的最自然的方法是，计算A-1，它的计算公式NN |*1AAAA*是A的伴随矩阵。第4页/共56页计算量|*1AAANNppppAAAAA2121) 1()det(|NNNNNNAAAAAAAAAA212222111211*ijjiijM) 1(A伴随矩阵代数余子式余子式Mij：在A中去掉i行j列以后剩下的行列式。第5页/共56页计算量 pi有N！种排列。 这意味着有N!个项相加，每一项是N个数的乘积。 当矩阵是20阶方阵，有2.4E18项。 假设一次乘法需要10个时钟周期（

3、具体数目依赖于cpu型号），计算一项需200个时钟周期，计算一个行列式需要4.8E20时钟周期。4.8Ghz的cpu需要E11秒，约E3年。NNppppAAAAA2121) 1()det(|第6页/共56页 行列式不能这么计算。 矩阵的逆不能这样求。第7页/共56页Gauss-Jordan方法 求A-1的一种最简单的实用方法是Gauss-Jordan方法。 它的做法是对矩阵A进行初等行变换第8页/共56页 定义：下面三种变换称为矩阵的初等行变换1.对调两行2.以数k（k不等于0）乘某一行中的所有元素3.把某一行所有元素的k倍加到另一行对应的元素上去。第9页/共56页借助于符号的推导 这三种初等

4、行变换中的每一种都可以用一个简单的矩阵T左乘A来表示。 例如，当N=3下面三个矩阵的作用是 使A的第三行乘以2，使第一行和第二行互换，从第三行减去第一行的一半。第10页/共56页 Gauss-Jordan方法是要找这种运算构成的一个序列 它作用于A左边把A化为单位矩阵，即 因此T就是A的逆矩阵。第11页/共56页实际的计算过程 用一个例子来说明。 考虑3阶方阵A和单位阵I 首先做一个变换，使A的第一列除第一个元素外其他元素都变成零。第12页/共56页 从第二行减去第一行的4倍。 然后，从第三行减去第一行的2倍 现在转而处理第二列第13页/共56页 把第二行乘以1/3加到第一行上。 第二行乘以-

5、1/6 最后处理第3列第14页/共56页 把第三行的1/3加到第一行和第二行上，并将第三行乘以-1 这样TA就化成了单位阵。从而求出来逆矩阵TI。第15页/共56页补充 可能发生这样的情况：在计算过程中的某处，TA中我们正在处理的那一列上的对角元素为零。 在这种情况下，两行交换将带给这个对角元素一个非零值。 如果没有行交换能够做到这一点，那么矩阵A是奇异的。 实际上，数值精确度要求行交换把正处理的那一列中绝对值最大的元素移到对角位置上。第16页/共56页 在求逆过程中，如果矩阵的各个元素的大小相差过于悬殊，还可能发生与数值舍入误差相联系的问题。 在此情况下，对各行或各列进行标度化（即应用第二种

6、初等变换），常常是有用的。第17页/共56页 在一些特殊情况下计算量可以减少。例如1.若A是一个对称矩阵，2.只对b求解x有兴趣而对A-1没兴趣。 对情况2，我们可以把T作用到b上，而不是作用到整个单位矩阵上。（利用结合律可证。）第18页/共56页 若我们只是感兴趣|A|，那么用行变换把TA变成下三角或者上三角形式就可以了。 此时行列式就是对角元素的乘积。第19页/共56页流程图第20页/共56页第21页/共56页第22页/共56页第23页/共56页作业 写一个程序包含调整最大对角元的过程，计算下面矩阵的逆矩阵。 结果是：3431223211115 . 235 . 1231第24页/共56页5

7、.2 三对角矩阵的本征值 现在讨论方阵A的本征值和本征向量的问题。 即：求满足下面本征方程的 和 等价地，本征值就是N次特征多项式的零点nn第25页/共56页 此处仅讨论实对称矩阵。 它的本征值全部为实数。 这也是在模拟物理系统时最常见的情况。第26页/共56页 此处考虑需要求出本征值（或许多个本征值），以及它们的本征向量的情况。 如果只需要求一个大矩阵的大的或最小的本征值，那么第7.4节的迭代方法效率可能更高。第27页/共56页 对A进行对角化的总方针是把这个问题化为求一个对称的三对角矩阵的本征值和本征向量。 这是一个容易处理的问题。 通过一系列正交变换，总可以把一个矩阵化为三对角矩阵。这将

8、在下节讨论。 本节先讨论如何求三对角矩阵的本征值。第28页/共56页 所谓三对角矩阵是：除对角线上的元素和相邻元素以外的一切元素都是零。665544332211aaaaaa66655655544544433433322322211211aaaaaaaaaaaaaaaa第29页/共56页 为求一个对称三对角矩阵的本征值，必须求特征方程的根。 这个特征方程的左边（特征多项式）由下面的行列式给出第30页/共56页 只要我们对于给定的 能够求出PA的数值，就可以使用第一章的求根方法求出特征多项式的零点，即A的本征值。 PA的数值可以通过递推方便地求出。 令 是前n行和前n列构成的 子行列式之值。 显然

9、)(nPnnANPP 第31页/共56页111)(AP2121222)()()(APAP)()()()(221,1nnnnnnnPAPAP第32页/共56页 需要搜索的范围：第33页/共56页 习题5.4 写一个程序，求一个三对角矩阵的全部本征值。计算一个特殊的三对角方阵。这个矩阵的元素为： 用解析的方法已知这个矩阵的本征值是第34页/共56页第35页/共56页第36页/共56页第37页/共56页5.3 化为三对角形式 为了应用上节讨论的求本征值的方法，必须把一般的实对称矩阵A化为三对角形式。 我们必须找出正交矩阵O使得 为三对角矩阵。 变换后的矩阵的本征值就是原来的矩阵的本征值。 问题就是，

10、对于一个特定的矩阵A，如何求O的精确形式。第38页/共56页Householder方法 Householder方法是把一个矩阵化为三对角形式的一种常用而且方便的方法。 它取矩阵O为N-2个正交矩阵的乘积： 其中每一个都依次把A的一行和一列变换成需要的形式。 最后的两行和两列不需要变换。 求出这种 的方法是简单的。第39页/共56页工作原理 求第一个正交变换 ，它要达到的效果是A11是原矩阵的相应元素，k(1)是一个可能不为零的元素，矩阵A(2)是把变换O1应用于A后得到的后N-1行和N-1列的结果。 然后对A(2)做相似的变换。第40页/共56页O1的形式 现在我们省略下标 取O的形式为 其中

11、 是一个N-1维的由0构成的列向量，P是一个N-1维对称矩阵。 为了使O是正交阵，P必须是正交阵。 一种取法是Householder矩阵第41页/共56页O1的形式 其中I是一个N-1维单位矩阵，u是一个N-1维单位向量， 上式中的 是u的外积，结果是一个N-1维矩阵。 矩阵P的元素为第42页/共56页n1时的On的形式 当n=2，取其形式为： 更大的n，以此类推。第43页/共56页求出On的精确形式 利用前述正交矩阵可得 其中 是A的第一列中除第一个元素外的全部元素 根据这个形式要求第44页/共56页 能把k和u都求出来吗？ 首先计算K和k 接下来计算u左右两边做内积，可得利用 值 求出u以

12、后我们就可以根据公式P=I-2uut计算出P。正负号如何取？处于数值稳定性的考虑，取使等号右边的值大的符号。第45页/共56页k和A(2) 刚才我们已经求出了k， 变换以后的子矩阵A(2)还没有求出来。 在计算A(2)时不必进行整个矩阵相乘 其中A是原来矩阵的末N-1行和N-1列的对称矩阵。 这里展开成只有矩阵向量相乘的形式。第46页/共56页n1时的On的形式 当n=2，取其形式为： 更大的n，以此类推。第47页/共56页 在接连变换的过程中，如果绝对值最大的对角元素处于待变换的子阵的左上角，那么数值稳定性会得到改进。 在每次实行正交变换之后进行适当的行变换和列变换，可以做到这一点。此时本征

13、值不变。第48页/共56页程序第49页/共56页求本征矢量 有了把A变换为三对角形式的算法，我们就完全知道了如何求一个实对称矩阵的本征值。 一旦求出了本征值，求本征向量就比较简单了。 这里选用的方法是“逆向量迭代法”。第50页/共56页逆向量迭代法 逆向量迭代法的工作过程如下： 令 是对与 相联系的本征向量的任意猜测。 通过迭代来改善这一猜测： 其中 是一个非零的小标量。 这个运算会加强 沿着真正本征矢量方向的分量。第51页/共56页逆向量迭代 考虑y是线性方程的解 b是一个随机矢量， 很靠近一个本征值 。 把y解出来，它靠近 的本征值。 用y的值代替b，解出新y，它更靠近真正的本征矢量。 用逆矩阵的表示方法，就得到逆向