孤立特征值情况的矩阵摄动法算例(matlab编程)

合集下载

matlab中eig函数用法

matlab中eig函数用法
Matlab的eig函数是一个多功能的命令，用于计算矩阵的特征值和特征向量。

它可以被用来求解一些振动类型的简小问题，解决斯坦纳网络的传播问题，判断系统的稳定性以及识别线性系统的状态。

eig函数是用来计算矩阵特征值和特征向量的函数。

简单地说，特征值是一个
复数，表示变换向量空间中的某一矢量。

特征向量是一个列向量，表示变换向量空间中的某一个方向。

eig函数可以求解一个n阶实矩阵A的n个特征值和n个特征
向量。

使用函数的方法如下：
[V,D] = eig(A)
其中A为一个实矩阵，V是A的特征向量矩阵，D是包含A的特征值的对角矩阵。

它们之间的关系是：
AV=VD
其中V是A的特征向量矩阵，D是A的特征值矩阵。

eig函数是一个有效的求解矩阵特征值和特征向量的函数，可以解决一些系统
分析和振动传播类型的问题，同时也可以用它来识别线性系统的状态、判断系统的稳定性等。

matlab 约束矩阵正定

matlab 约束矩阵正定
约束矩阵正定性在优化问题中是非常重要的，因为正定矩阵具有很多良好的性质，能够保证优化问题的最优解存在且唯一。

在MATLAB 中，判断一个矩阵是否正定可以使用eig函数，该函数计算矩阵的特征值和特征向量。

如果一个矩阵的所有特征值均大于零，则该矩阵为正定矩阵。

因此，我们可以通过eig函数计算出约束矩阵的特征值，判断它们是否全部大于零，从而判断约束矩阵是否正定。

下面是一个简单的MATLAB代码示例：
C = [1 2 3; 2 5 6; 3 6 9]; % 定义一个矩阵C
[eigVec, eigVal] = eig(C); % 计算矩阵C的特征值和特征向
量
if all(diag(eigVal) > 0) % 判断特征值是否全部大于零
disp('矩阵C为正定矩阵'); % 如果全部大于零，则矩阵C为正定矩阵
else
disp('矩阵C不是正定矩阵'); % 如果存在小于等于零的特征值，则矩阵C不是正定矩阵
end
通过这个简单的代码示例，我们可以看到MATLAB如何判断约束
矩阵的正定性。

在优化问题中，正确地判断约束矩阵是否正定，可以帮助我们更好地处理优化问题，得到更优的解。

matlab特征分解

matlab特征分解一、Matlab特征分解简介特征分解是矩阵理论中的一个重要概念，可以将一个方阵分解为对角矩阵和相应的特征向量矩阵。

在Matlab中，可以使用eig函数进行特征分解。

二、Matlab特征分解函数eig的使用方法1.基本语法：[V,D] = eig(A)，其中A为待分解的方阵，V为特征向量矩阵，D为对角矩阵。

2.示例代码：A = [1 2 3; 4 5 6; 7 8 9]; %定义一个3*3的方阵[V,D] = eig(A); %进行特征分解disp(V); %输出特征向量矩阵disp(D); %输出对角矩阵三、Matlab特征分解的应用场景1.求解线性微分方程组；2.求解差分方程组；3.求解最小二乘问题；4.求解最大化或最小化目标函数等。

四、Matlab特征值和特征向量的性质1.若A是实数矩阵，则其复共轭转置也是其伴随矩阵，即A* = A'；2.A与其伴随矩阵A'之积是一个Hermite（共轭对称）矩阵，即AA' = A'A；3.特征向量矩阵V的列向量是相互正交的，即V'V = I；4.特征向量矩阵V的逆矩阵等于其转置矩阵，即V^-1 = V'。

五、Matlab特征分解的相关函数1.eig函数：用于计算方阵的特征值和特征向量；2.svd函数：用于计算奇异值分解（SVD）；3.qr函数：用于计算QR分解。

六、Matlab特征分解的注意事项1.在进行特征分解时，需要注意输入的方阵必须是方阵，否则会出现错误；2.在使用eig函数进行特征分解时，需要注意输出结果中特征向量矩阵和对角矩阵并不是按顺序排列的，需要手动排序后才能正确使用。

七、Matlab特征分解实例下面给出一个简单的实例来说明Matlab中如何使用eig函数进行特征分解：A = [4 2; 1 3]; %定义一个2*2的方阵[V,D] = eig(A); %进行特征分解disp(V); %输出特征向量矩阵disp(D); %输出对角矩阵运行以上代码后，输出结果为：V =0.8944 0.70710.4472 -0.7071D =2 00 5其中，V为特征向量矩阵，D为对角矩阵。

第7章矩阵特征值和特征向量的数值解法

n
i xi
i 1
n
i Ak xi
i 1
n
ii k xi
i 1
1k [1x1
n i ( i )k
i2
1
xi ]
(7.1.4)
主要用于求矩阵按模最大的特征值和相应的特征向量。设矩阵 A 的
n 个特征值 i (i 1,2,..., n) 满足：
| 1 || 2 || 3 | . . .| n | （7.1.1）
3 2.689 319 6.737 850 6.747 559 0.398 562 0.998 561 1.000 000
4 1.595 686 2.379 870 2.381 309 0.670 088 0.999 396 1.000 000
5 2.680 956 6.772 616 6.723 220 0.398 761 0.999 910 1.000 000
7 9.605 572,5.817 228,-3.778 139 1.000 000，0.605 777，-0.394 369 9.605 572
8 9.605 567,5.816 808,-3.788 717 1.000 000，0.605 566，-0.394 429 9.605 567
由表 7.1.1 知， m8 m8 105 ，故取 1 m8 9.605567 ，
u(k) ma x(u
(k
)
)
的最大分量为
1，即完成了规范化。
7.1.1 幂法原理及实用幂法
由于 v(k) 中最大分量为 1，即 max（ v(k) ）=1，故
v(k)
Ak u ( 0) max( Aku(0) )
(7.1.6)
由式（7.1.4）有

(完整版)关于matlab中的eig函数(求特征值和特征向量)

关于matlab中的eig函数（求特征值和特征向量）在MATLAB中，eig用途：Find eigenvalues（特征值）and eigenvectors（特征向量），常用的调用格式有5种：(1) E=eig(A)：求矩阵A的全部特征值，构成向量E。

（注意，第一列为对应第一个特征值的特征向量）(2) [V,D]=eig(A)：求矩阵A的全部特征值，构成对角阵D，并求A的特征向量构成V的全部列向量。

(3) [V,D]=eig(A,'nobalance')：与第2种格式类似，但第2种格式中先对A 作相似变换后求矩阵A的特征值和特征向量，而格式3直接求矩阵A的特征值和特征向量。

(4) E=eig(A,B)：由eig(A,B)返回N×N阶方阵A和B的N个广义特征值，构成向量E。

(5) [V,D]=eig(A,B)：由eig(A,B)返回方阵A和B的N个广义特征值，构成N×N 阶对角阵D，其对角线上的N个元素即为相应的广义特征值，同时将返回相应的特征向量构成N×N阶满秩矩阵，且满足AV=BVD。

Syntax（句法）如下：d = eig(A)d = eig(A,B)[V,D] = eig(A)[V,D] = eig(A,'nobalance')[V,D] = eig(A,B)[V,D] = eig(A,B,flag)d = eig(A)和 [V,D] = eig(A)最为常用注意，第一列为对应第一个特征值的特征向量，比如：B=rand(4)B =0.5653 0.7883 0.1365 0.97490.2034 0.5579 0.3574 0.65790.5070 0.1541 0.9648 0.08330.5373 0.7229 0.3223 0.3344>> [a,b]=eig(B) %求矩阵B的全部特征值，构成对角阵b，并求B的特征向量构成a的列向量。

关于matlab中的eig函数(求特征值和特征向量)

关于matlab中的eig函数（求特征值和特征向量）在MATLABxx，eig用途：Find eigenvalues（特征值）and eigenvectors（特征向量），常用的调用格式有5种：(1)E=eig(A)：求矩阵A的全部特征值，构成向量E。

（注意，第一列为对应第一个特征值的特征向量）(2)[V,D]=eig(A)：求矩阵A的全部特征值，构成对角阵D，并求A的特征向量构成V的全部列向量。

(3) [V,D]=eig(A,'nobalance')：与第2种格式类似，但第2种格式中先对A作相似变换后求矩阵A的特征值和特征向量，而格式3直接求矩阵A的特征值和特征向量。

(4)E=eig(A,B)：由eig(A,B)返回N×N阶方阵A和B的N个广义特征值，构成向量E。

(5)[V,D]=eig(A,B)：由eig(A,B)返回方阵A和B的N个广义特征值，构成N×N阶对角阵D，其对角线上的N个元素即为相应的广义特征值，同时将返回相应的特征向量构成N×N阶满秩矩阵，且满足AV=BVD。

Syntax（句法）如下：d = eig(A)d = eig(A,B)[V,D] = eig(A)[V,D] = eig(A,'nobalance')[V,D] = eig(A,B)[V,D] = eig(A,B,flag)d = eig(A)和[V,D] = eig(A)最为常用注意，第一列为对应第一个特征值的特征向量，比如：B=rand(4)B =0.5653 0.7883 0.1365 0.97490.2034 0.5579 0.3574 0.65790.5070 0.1541 0.9648 0.08330.5373 0.7229 0.3223 0.3344>> [a,b]=eig(B) %求矩阵B的全部特征值，构成对角阵b，并求B的特征向量构成a的列向量。

matlab如何求矩阵方程的基础解系

矩阵方程的基础解系是指由矩阵方程的特解和齐次方程的通解组成的一种特殊解系，它包括了矩阵方程的所有解。

在 MATLAB 中，我们可以利用一些基本的函数和方法来求解矩阵方程的基础解系。

接下来我们将介绍如何使用 MATLAB 求解矩阵方程的基础解系的基础方法。

1. 准备工作在使用 MATLAB 求解矩阵方程的基础解系之前，首先需要确保MATLAB 环境已经安装并且处于可用状态。

需要定义矩阵系数和常数向量。

2. 使用 `\` 求解特解在 MATLAB 中，可以使用 `\` 运算符来求解线性方程组的特解。

具体而言，如果矩阵方程为 `A * X = B`，那么可以使用 `X = A \ B` 来求解特解。

3. 求解齐次方程的通解同样地，在 MATLAB 中，可以使用 `null` 函数来求解齐次方程的通解。

具体而言，如果矩阵方程为 `A * X = 0`，那么可以使用 `N =null(A,'r')` 来求解齐次方程的通解。

4. 合并特解和通解将特解和通解合并起来，就得到了矩阵方程的基础解系。

具体而言，可以使用 `X = special + nullspace` 来完成合并操作。

总结通过以上步骤，我们可以使用 MATLAB 求解矩阵方程的基础解系。

首先求解特解，然后求解齐次方程的通解，最后合并特解和通解即可得到基础解系。

这样的方法可以帮助我们更好地理解和求解矩阵方程的基础解系，同时也展现了 MATLAB 在数学计算方面的强大功能和灵活性。

在实际问题中，矩阵方程的基础解系往往是非常重要的，它可以帮助我们求解各种复杂的线性方程组和矩阵方程，从而进一步应用到实际的工程和科学问题中。

掌握如何使用 MATLAB 求解矩阵方程的基础解系是非常有必要的。

希望本文的介绍能够对读者有所帮助，也欢迎大家多多交流和共享关于 MATLAB 求解矩阵方程的经验和技巧。

在实际工作和学习中，矩阵方程的基础解系是一个非常重要的概念。

基于Matlab的数值分析实验的设计r——以矩阵的奇异值分解为例

基于Matlab的数值分析实验的设计r——以矩阵的奇异值分解为例杨海涛【摘要】基于Matlab的数学实验,目的是为了提高普通高校学生的数学兴趣和动手能力,为培养解决实际数学问题的能力而设计的实验方法和步骤.利用矩阵的奇异值分解svd函数的数值算法设计了一种在实验课上用svd算法仿真的实验方案.【期刊名称】《内蒙古民族大学学报（自然科学版）》【年(卷),期】2016(031)006【总页数】4页(P465-468)【关键词】数学实验;数值计算;Matlab;奇异值分解【作者】杨海涛【作者单位】内蒙古民族大学数学学院,内蒙古通辽 028043【正文语种】中文【中图分类】O241.4利用计算机计算解决实际问题时需要建立数学模型，同时也需要相应数学软件的支持，为了更好的掌握数学软件这一有利的工具，本文通过Matlab中的矩阵的奇异值分解〔1〕svd函数设计了如何学习和掌握该函数的实验方法，从而可以进一步学习其它函数.目前，Matlab在计算仿真和图形图像的压缩处理〔2～4〕有着广泛的应用，希望这个实验方法能够给正在学习数学软件的学生起到抛砖引玉的作用.1.1 理论基础任意复长方矩阵的奇异值分解定理如下（Autonne-Eckart-Young定理）：〔5〕设A为一个m×n的实矩阵（m≥n），则存在正交矩阵V和U，使得：其中∑=diag(σ1,σ2,…,σr)，而且σ1≥σ2≥…≥σr≥0.奇异值分解在统计分析、信号与图像处理、系统理论和控制中有着广泛的应用. 1.2 SVD的数值算法〔6〕求解一般矩阵奇异值问题最常用的算法可分为两大类：一类是QR分解的方法；另一类是Jacobi旋转的方法.奇异值分解的QR分解法分为两个阶段：一是矩阵的二重对角化，即利用Householder变换把矩阵A变换成二重对角矩阵.这个过程就是追赶（chasing）过程〔7〕.二是QR分解，即保持二重对角矩阵形式不变，利用正交变换使对角线元素逐渐减小，使得矩阵接近对角矩阵.1.2.1 svd函数.对于matlab中的svd函数不是开源的，但通过调研得知该函数的算法就是采用propack中的lansvd函数的算法.其算法过程如下：对于任意的矩阵Am×n可以利用lanczos过程将其双对角化，即存在m阶正交矩阵Q和n阶正交矩阵W满足：其中为一双对角阵，且αi, βi∈R (i=1, 2, …, n)，这里假设了m＞n.上面的分解是逐步的利用Householder左右变换得到的.算法中利用了已知函数完成的双对角化过程.这里如果我们已经知道Bn的svd分解为其中Vn+1为n+1阶正交矩阵，的广义对角矩阵Un为n阶正交矩阵.则有令V=Qn+1Vn+1，其前r列为对应的左奇异向量，令U=WUn，其前r列为对应的右奇异向量.因此算法的关键在于计算矩阵Bn的svd分解.因为Bn是一个双对角矩阵，它的奇异值分解计算比较简单，运算量较小，利用Givens变换，算法复杂度仅为O(n).在我们这个程序中，直接调用了库函数svd 计算Bn的奇异值.1.2.2 svds函数.对于matlab中的svds函数，算法流程如下：（1）对于m×n阶矩阵A，构造(m+n)×(m+n)阶矩阵，其中O是稀疏矩阵. （2）利用eigs函数计算出B矩阵的特征值矩阵D，以及D对应的正交阵W. （3）根据输入参数，按要求输出最大特征值、最接近某参数的特征值或对特征值进行排序等.（4）根据设定的阀值，对矩阵进行筛选得到U矩阵，对矩阵筛选得到V矩阵，对D矩阵进行筛选得到S矩阵.（5）根据S中元素的大小对U、V、S进行排序.（6）输出结果.有如下J-W定理（Jordan-Wielandt定理）：若σ1≥σ2≥…≥σp-1≥σp是矩阵Am×n的奇异值（其中，p=min{m,n}），则上述B矩阵具有2p个非零特征值-σ1,…,-σp,σp,…,σ1和|m-n|个零特征值，与非零特征值（±σj）相对应的特征向量为若m≠n，另有特征向量分别对应于m＞n或m＜n.svds函数的基本思想就是构造出这样的矩阵B，并且求出其特征值和特征向量，然后根据函数重构的条件（所需奇异值的数量、最大最小等），求出所需的奇异值和对应的向量.2.1 矩阵生成生成10个不同的m×n的矩阵，这里取m=120 , n=50，对这10个矩阵进行svd分解.矩阵生成的方法是先做出矩阵的奇异值，然后用QR分解做出10个120阶正交阵和10个50阶正交阵，再利用奇异值分解公式生成10个矩阵A.10个矩阵的奇异值（均为25个）如表1所示.2.2 仿真结果2.2.1 计算时间svd和svds两种方法对十个矩阵进行奇异值分解消耗的时间如表2，单位是秒.可见，svd函数运行时间一般要小于svds函数运行时间.2.2.2 最小奇异值计算准确率.通过对比最小特征值的计算准确度，将每种算法特征值的最小值与矩阵最开始生成时的最小特征值进行比较.比较结果两种方法对10个矩阵的svd分解运算最小特征值的误差均几乎为0，都在10-15次方的量级左右，见表3.其中delta_min_svd和delta_min_svds是matlab程序中两个变量值.这个误差已经和截断误差非常接近了，因此可以认为两种方法都能较好的进行奇异值计算.2.2.3 左右奇异值向量计算准确率.由于svd分解后U、V矩阵的不唯一性，这里我们只看生成的U、V的正交性是不是依然比较完好.通过计算每一列和它的转置相乘减去单位阵也都在10-15这样的量级，因此svd和svds对于正交性的保持都是非常好的.综上所述，通过上述设计的方法和步骤，完全可以掌握matlab中的M-文件和.mat文件的使用方法，并且能熟练掌握奇异值分解的理论知识和具体计算方法，才能为以后svd应用打下基础，这样的实验方法才有实际的意义，才是真正的综合设计实验.本文的宗旨在于引导学生如何用数学软件完成一个实验设计，同时掌握这个实验的目的、相关内容实验的方法和手段.〔1〕尹芳黎，杨雁莹，王传栋，等.矩阵奇异值分解及其在高维数据处理中的应用〔J〕.数学的实践与认识，2011，45（15）: 171-176.〔2〕胡乡峰，卫金茂.基于奇异值分解（SVD）的图像压缩〔J〕.东北师大学报:自然科学版，2006，38（3）:36-39.〔3〕曾超，张陈，徐永利.基于奇异值分解的图像压缩及其matlab实现〔J〕.科技信息，2010，14:484.〔4〕向培素.一种自适应AP算法的matlab实现〔J〕.西南民族大学学报:自然科学版，2014，40（16）:877-882.〔5〕张贤达.矩阵分析与应用〔M〕.北京：清华大学出版社，2004.358-360. 〔6〕姜健飞，胡良剑，唐俭.数值分析及其matlab实验〔M〕.北京：科学出版社，2004.6，219-224.〔7〕山晓东，李园.一类不可约L-矩阵的预条件AOR迭代法〔J〕.湖北民族学院学报:自然科学版，2014，32（2）:128-132.。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

设原始特征值问题是 )~1(00000nixMxKiii (1)

相应的正交规范条件是 )~1,(000njixMxijjTi (2)

先考虑一个五个自由度的系统．其中，0M、0K分别是原系统的nn阶对称质量阵、刚度阵，0i是特征值，且200ii，

0i是固有频率，0ix是相应的特征向量。

011110.5

M

02112112112111K









系统结构参数改变后，相应的质量阵、刚度阵均有相应的变化，设系统修改后的质量阵、刚度阵分别为

1010KKKMMM



(3)

其中取

111110.5

M 110000K





显然，新系统（即结构修改后的系统或称摄动系统）的特征值问题及相应的正交规范条件是

iiixMMxKK)()(1010

(4)

ijjTixMMx)(10 (5) 其中，i为新系统的特征值，ix 是相应的特征向量。根据摄动理论，可将i、ix按小参数展开成幂级数（因此，胡海昌院士称其为小参数法），即 )(32210iiii (6)

)(32210iiiixxxx (7)

现在来确定1i，2i，1ix，2ix，将（6）、（7）式代入（4）式得 ))()(())((2210102210221010iiiiiiiiixxxMMxxxKK 展开上式，略去)(3，并比较的同次幂系数可得 000000:iiixMxK (8-1)

00101010001101:iiiiiiiixMxMxMxKxK (8-2)

00201110111020011202:iiiiiiiiiiiixMxMxMxMxMxKxK (8-3)

将（6）、（7）式代入（5）式得

ijjjjTiiixxxMMxxx))(()(2210102210 展开上式，并令ji（从后面的分析可得，ji的关系式在推导中没有被直接利用），比较的同次幂系数可得

1:0000iTixMx (9-1)

0:0100011001iTiiTiiTixMxxMxxMx (9-2)

0:0021010111102002iTiiTiiTiiTiiTixMxxMxxMxxMxxMx (9-3)

至此，我们已推得了进行摄动分析（即求解i，ix中的1i，2i，1ix，2ix）所需的全部基本方程。实际上，(8-1)、(9-1)即为（1）、（2）式，是显然满足的。于是，可在原系统的特征值问题的基础上，通过(8-2)、(9-2)求解一阶摄动1i，1ix，通过(8-3)、(9-3)求解二阶

摄动2i，2ix，代回（6）、（7）式，即求得i，ix，且具有)(3精度。

一阶摄动公式可将1ix表示成原系统模态（向量）02010,,,nxxx的线性组合，即 nkkkixx1011 (10)

其中，1k是n个待定系数。将（10）代入(8-2)，并左乘Tkx0，得

0001010010100001010100iTkiiTkinkkkTkiiTknkkkTkxMxxMxxMxxKxxKx 利用正交规范关系（2）式，上式可简化成

000101000101001iTkiiTkiikiTkkkxMxxMxxKx (11)

ki 时，00ik，1000iTkxMx，由（11）式可得特征值的一阶摄动为

010101)(iiTiixMKx (12)

ki时，0000iTkxMx，则由（11）式可得

)()(00010101kixMKxkiiiTkk

 (13)

至此，n个待定系数1k中，只有1i尚未确定，现在来求1i。用00MxTi左乘（10）式两边，得

110001100inkkTikiTixMxxMx (14)

（14）式转置，并考虑到0M对称，且1i是一个常数，于是得 1001iiTixMx (15)

（14）、（15）式代入(9-2)可得 010121iTiixMx (16)

由（13）、（16）两式即完全确定了)~1(1nkk，也就确定了（10）式的1ix，而1i由（12）式给出，于是可得一阶摄动公式为







0010,1000010101010101)(21)()(iiTinikkkiiiTkiiiTiixxMxxxMKxxxMKx



(17)

二阶摄动公式为了得到更精确的摄动解，需要用到二阶摄动。根据展开定理，将2ix按02010,,,nxxx展开 nkkkixx1022 (18)

将（18）式代入(8-3)，并左乘Tkx0，得

0002011101110010200011010200)(iTkiiiiiiiTknkkkTkiiTknkkkTkxMxxMxMxMxxMxxKxxKx







利用（2）式，上式可简化成 000201110111000211002)(iTkiiiiiiiTkikiTkkkxMxxMxMxMxxKx



(19)

ki 时，00ik，1000iTkxMx，由（19）式可得特征值的二阶摄动为

)(0111011101102iiiiiiiTiixMxMxMxKx (20)

ki时，0000iTkxMx，则由（19）式可得

)()(000111011101102kixMxMxMxKxkiiiiiiiiTkk

 (21)

现在来确定最后一个系数2i，为此，用00MxTi左乘（18）式两边，得 210200200inkkkTiiTixMxxMx (22)

（22）式转置后，注意到0M对称，2i是一个数，可得 2002iiTixMx (23)

（22）、（23）式代入(9-3)可得 )(210111011102iTiiTiiTiixMxxMxxMx (24)

于是就得到了二阶摄动公式







0011101110,100001110111011020111011101102)(21)()(iiTiiTiiTinikkkiiiiiiiiTkiiiiiiiiTiixxMxxMxxMxxxMxMxMxKxxxMxMxMxKx



编程如下 M0=diag([1,1,1,1,0.5]); K0=[2,-1,0,0,0; -1,2,-1,0,0; 0,-1,2,-1,0; 0,0,-1,2,-1; 0,0,0,-1,1]; M1=diag([1,1,1,1,0.5]); K1=diag([-1,0,0,0,0]); [x0,r0]=eig(K0,M0); %%[V,D] = eig(A,B) produces a diagonal matrix D of generalized eigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that A*V = B*V*D . %%x0为特征向量矩阵，r0为特征值对角阵 %%下面来计算一阶摄动 for m=1:5 r1(m,m)=x0(:,m)'*(K1-r0(m,m)*M1)*x0(:,m); x1(:,m)=-0.5*x0(:,m)'*M1*x0(:,m)*x0(:,m); for n=1:5 if n~=m

x1(:,m)=x1(:,m)+x0(:,n)'*(K1-r0(m,m)*M1)*x0(:,m)/(r0(m,m)-r0(n,n))*M0(:,n); end end end %%再来计算二阶摄动 for m=1:5 r2(m,m)=x0(:,m)'*(K1*x1(:,m)-r0(m,m)*M1*x1(:,m)-r1(m,m)*M0*x1(:,m)-r1(m,m)*M1*x0(:,m));