计算方法大作业

1、XXXX XXXX XXXX计算方法大作业班级XXXXXX学号XXXXXXX任课老师贺力平姓名XXX2013年 12月实验一 幂法与矩阵特征值1. 幂法求主特征值思路幂法的主要思想就是对假设的任意初始列向量作用n次A矩阵（左乘A矩阵）后，初始向量就接近A矩阵的主特征值对应的特征向量。由于左乘n次A矩阵有可能会造成计算量溢出，所以每次都对列向量作归一化处理。开始输入矩阵A,v(0)v(1)=Av(0)，i=0u(i)=v(i)/max(v(i)i=i+1v(i+1)=Au(i)max(v(i+1)-max(v(i)<?noyesu=v(i)/max(v(i)结束输出矩阵u,=uTAu/uT

2、u2. 程序代码function = mifa( A,v )%UNTITLED Summary of this function goes here% Detailed explanation goes hereA = 1 21 2 3 34 5 2 1 54 2 6 4 2 2 59 0; v=1 1 1 1'u(:,1)=v(:,1);fori=1:100v(:,i+1)=A*u(:,i);if abs(max(v(:,i+1)-max(v(:,i)<10-4breakend u(:,i+1)=v(:,i+1)/max(v(:,i+1);enddisp(u(:,i);disp

3、(max(v(:,i);k=u(:,i)'*A*u(:,i)/(u(:,i)'*u(:,i);disp(k);x,c=eig(A);disp(c);disp(x);disp(i);end3. 结果比较和结论初始矩阵 A = 1 21 2 334 5 2 154 2 6 4 2 2 59 0;初始向量v=1 1 1 1'幂法求得特征向量12.5147 14.9140 25.6935 39.6330归一化后特征向量 0.3158 0.3763 0.6483 1.0000列向量最大值近似主特征值39.6330Rayleigh商求出主特征值39.6330用eig()函数算出的特

4、征值和特征向量 39.6331 0 0 0 0 -18.2401 + 7.4985i 0 0 0 0 -18.2401 - 7.4985i 0 0 0 0 8.8471 -0.2450 -0.0471 + 0.0965i -0.0471 - 0.0965i -0.0574 -0.2919 0.1147 - 0.0939i 0.1147 + 0.0939i -0.1747 -0.5029 0.2862 - 0.1187i 0.2862 + 0.1187i 0.1535 -0.7758 -0.9330 -0.9330 0.9709 达到精度要求所需次数：17结论：可以看出初始列向量经过多次迭代后，

5、用幂法求出的特征值和用eig()函数求出的A的特征值，满足计算精度在，并且特征向量也具有数乘关系。另外发现取特征向量v的列最大值近似的特征值和用Rayleigh商求出的特征值相同。实验二 QR分解与特征值开始1. 幂法求主特征值思路重复直到a->上三角i=i+1<long(a)结束输出aa=R*QQ=h(i+1)R=aHousehold镜像矩阵h(i+1)=h* h(i)Household镜像矩阵a=h* aa的第i列作household变换输入矩阵a2. 源程序Household函数代码function Q R = hshd( a)%UNTITLED Summary of thi

6、s function goes here% Detailed explanation goes heren=length(a);h=eye(size(a);fori=1:n-1if(a(i,i)>0)sgn=1;elsesgn=-1;endx=sgn*(sum(a(i:n,i).2)0.5;u=a(i:n,i)+x*eye(n-i+1,1);ro=norm(u,2)2/2;H=eye(n-i+1)-u*u'/ro;H1=eye(size(a);H1(i:n,i:n)=H;a=H1*a;h=H1*h;endQ=inv(h);R=a;returnendQR分解代码function =

7、 qrfenjie(a)%UNTITLED Summary of this function goes here% Detailed explanation goes herefori=1:50Q,R=hshd(a);a=R*Q;enddisp(a)end3. 结果与结论假设a矩阵为如下3*3矩阵a = 2 2 10 1 43 0 2计算出a的特征值为5，+2.2361i，-2.2361i三个用QR分解逆序相乘法后得到矩阵 5.0000 -1.0690 1.2344 -0.0000 0.3571 2.7630-0.0000 -1.8558 -0.3571三对角线上要么是实特征值，要么是两两共轭

8、的虚特征值。对于上述矩阵，有实特征5和2*2矩阵构成的虚特征值，求得虚特征值正好为2.2361i。如果最后求得的矩阵为上三角矩阵，那么特征值就是对角线上元素，而且特征值均为实数。用eig()函数求得矩阵特征值如下5.0000 0.0000 + 2.2361i0.0000 - 2.2361i结论：可以看到两种算法在10-4次精度上特征值是相等的，这也验证了QR算法的正确性。实验三 数值分析在潮流计算方面的应用1 潮流计算的计算机算法简介潮流计算的计算机算法是以电网络理论为基础的，应用数值计算方法求解一组描述电力系统稳态特性的方程。从数学上讲是一组多元的非线性方程式的求解问题，这类方程的求解过程都

9、离不开迭代。由于电力系统结构及参数的一些特点，同时随着电力系统不断扩大，潮流问题的方程式的阶数也越来越高，这样的非线性方程式并不是任何数学方法都能保证给出正确答案的。这种情况就成为促使电力系统计算人员不断寻求新的且更可靠方法的一个重要因素。电网潮流计算的性能优劣一般依据的是能否可靠收敛，计算速度的快慢，内存占有多少，使用是否方便灵活，调整和修改是否容易，是否满足工程需要等来判别，其中以是否可靠收敛作为评价的主要标准。常用的分析法包括高斯-塞德尔法、牛顿-拉夫逊潮流算法、快速解耦算法（PQ 分解法）等。2 潮流计算的约束条件电力系统运行必须满足一定技术和经济上的要求。这些要求够成了潮流问题中某些

10、变量的约束条件，常用的约束条件如下：2.1节点电压应满足: （2.1） 2.2节点的有功功率和无功功率应满足： （2.2） 2.3节点之间电压的相位差应满足： （2.3） 3 节点导纳矩阵的形成与修改3.1 节点电压方程 （1）自、互导纳的物理意义自导纳在数值上等于与该节点I直接连接的所有支路导纳的总和。如。互导纳在数值上等于连接节点、支路导纳的负值，即。如。（2）节点导纳矩阵YB为对称方阵。（3）节点导纳矩阵YB为稀疏矩阵。（4）节点导纳矩阵具有对角优势。3.2 节点导纳矩阵的形成 用直接形成法形成节点导纳矩阵YB。节点导纳矩阵即可根据自导纳和互导纳的定义直接形成，也可用支路节点关联矩阵计算

11、。 3.3 节点导纳矩阵的修改（1）从原有网络引出一支路，同时增加一节点，节点导纳矩阵将增加一阶。新增的对角元，；新增的非对角元，；原有矩阵中的对角元将增加 ，。（2）在原有网络的节点、之间增加一支路。 ，（3）在原有网络的节点，之间切除一支路，（4）原有网络的节点、之间的导纳由改变为：，（5）原有网络节点i、j之间变压器的变比由改变为 ；4 牛顿-拉夫逊法（直角坐标）4.1概述1. 牛顿-拉夫逊法的意义和推导过程把按泰勒级数在点展开 修正方程 （4.1）2牛顿拉夫逊法的特点(1)牛顿-拉夫逊法是迭代法，逐渐逼近的方法；(2)修正方程是线性化方程，它的线性化过程体现在把非线性方程在按泰勒级数展

12、开，并略去高阶小量；(3)用牛顿拉夫逊法解题时，其初始值要求严格(较接近真解)，否则迭代不收敛。3多变量非线性方程的解牛顿拉夫逊法修正方程： （4.2）缩写为 (4.3) 4.2潮流计算时的修正方程（直角坐标）PQ节点 (4.4)PV节点 (4.5) 平衡节点只设一个,电压为已知,不参见迭代,其电压为 (4.6) 修正方程 (4.7) (4.8) 4.3雅可比矩阵各元素当时, 雅可比矩阵中非对角元素为 (4.9)当时,雅可比矩阵中对角元素为 (4.10)图1 牛顿-拉夫逊法计算步骤5 用MATLAB进行编程 牛顿-拉夫逊法（直角坐标）5.1 某电网接线图及给定的参数其中，1,2,3,4为PQ节

13、点，5为平衡节点各支路阻抗：Z12=Z21=0.06+j0.18 Z13=Z31=0.06+j0.18 Z14=Z41=0.04+j0.12Z15=Z51=0.02+j0.06 Z23=Z32=0.01+j0.03 Z25=Z52=0.08+j0.24 Z34=Z43=0.08+j0.24 各节点输出功率1：-0.2-j0.22: 0.45+0.15 3: 0.4+j0.054: 0.6+j0.1 5: 05.2 潮流计算计算机算法流程图开始形成节点导纳矩阵输入原始数据设节点电压，i=1,2,n,is置迭代次数置节点号i=1按式（4.9），（4.10）计算雅克比矩阵元素按式（4.4）计算节点的，节点的，求解修正方程式，得，雅克比矩阵是否已全部形成？计算平衡节点及PV节点功率求，迭代次数 k=k+1i=i+1？潮流计算完成计算各节点电压的新值：小结数值分析主要的思路就是用数值解代替精确解，从而寻找在一定精度内可以接受的数值算法。这是一种重要的简化