1、1用程序实现节点导纳矩阵的形成：主程序：%节点导纳矩阵的形成%网络的节点数为 nn=input(网络的节点数);A=input(节点相关性矩阵 );%1 表示两节点相连，0 表示两节点不连Y=zeros(n,n);for i=n:-1:1for j=i-1:-1:1if (A(i,j)ijY(i,j)=input(i,j 支路之间的导纳值);Y(j,i)=Y(i,j);endendY(i,i)=sum(Y(i,:);endY%对节点导纳矩阵的修改%0 表示不修改；%1 表示给原有网络增加节点，且节点处的阻抗为 Zij;%2 表示在节点 i,j 之间增加一条阻抗为 Zij 支路；%3 表示在原有
2、网络 i,j 之间切除一条阻抗为 Zij 支路；%4 表示原有网络 i,j 之间支路阻抗由 Zij 变成 Z1ijm=1;for m=1:100m=input(导纳矩阵修改的情况分类);switch(m)case(0)breakcase(1) i=input(在第 i 点增加阻抗 ); Zij=input(增加节点的阻抗 );Y(n+1,:)=zeros(1,n);Y(:,n+1)=zeros(n+1,1);Y(n+1,n+1)= Y(i,i);Y(i,n+1)=-1/Zij;Y(n+1,i)= Y(i,n+1);Y(i,i)= Y(i,i)+ Y(i,i);Y2case(2)i=input(
3、请输入节点 i);j=input(请输入节点 j);Zij=input(请输入在节点 i,j 间增加的阻抗);Y(i,i)= Y(i,i)+1/Zij;Y(j,j)=Y(j,j)+1/Zij;Y(i,j)=Y(j,j)-1/Zij;Y(j,i)=Y(i,j);Ycase(3)i=input(请输入节点 i);j=input(请输入节点 j);Zij=input(请输入在节点 i,j 间增加的阻抗);Y(i,i)= Y(i,i)-1/Zij;Y(j,j)=Y(j,j)-1/Zij;Y(i,j)=Y(j,j)+1/Zij;Y(j,i)=Y(i,j);Ycase(4)i=input(请输入节点 i)
4、;j=input(请输入节点 j);Z1ij=input(请输入在节点 i,j 间变换的阻抗);Y(i,i)= Y(i,i)-1/Y(i,j)+1/Z1ij;Y(j,j)=Y(j,j)-1/Y(i,j)+1/Z1ij;Y(i,j)=Y(j,j)+1/Y(i,j)-1/Z1ij;Y(j,i)=Y(i,j);Y otherwisedisp(输入数据有误，请重新输入数据);endend运行结果：网络的节点数 3节点相关性矩阵1,1,1;1,1,0;1,0,1i =3j =1i,j 支路之间的导纳值 5i =2j =31i,j 支路之间的导纳值 2Y =7 2 52 2 05 0 5导纳矩阵修改的情况
5、分类 1在第 i 点增加阻抗 2增加节点的阻抗 4Y =7.0000 2.0000 5.0000 02.0000 4.0000 0 -0.25005.0000 0 5.0000 00 -0.2500 0 2.0000导纳矩阵修改的情况分类 2请输入节点 i3请输入节点 j1请输入在节点 i,j 间增加的阻抗 8Y =7.1250 2.0000 7.0000 02.0000 4.0000 0 -0.25007.0000 0 5.1250 00 -0.2500 0 2.0000导纳矩阵修改的情况分类 3请输入节点 i2请输入节点 j3请输入在节点 i,j 间增加的阻抗 1Y =7.1250 2.0
6、000 7.0000 02.0000 3.0000 5.1250 -0.25007.0000 5.1250 4.1250 00 -0.2500 0 2.0000导纳矩阵修改的情况分类 4请输入节点 i3请输入节点 j3请输入在节点 i,j 间变换的阻抗 6Y =7.1250 2.0000 7.0000 02.0000 3.0000 5.1250 -0.25007.0000 5.1250 4.0542 00 -0.2500 0 2.0000导纳矩阵修改的情况分类 7输入数据有误，请重新输入数据导纳矩阵修改的情况分类 04三角分解法求解线性方程组：主程序：%用三角分解法解线性组方程 AX=bn=i
7、nput(系数矩阵的阶数);A=input(线性方程组的系数矩阵);b=input(线性方程组的常数项);L,U,m=zhjLU(A,n);y=zeros(1,4);if (m=1)x=b/L/U;x end调用程序：function L,U,m=zhjLU(A,n)RA=rank(A); if RA=ndisp(请注意：因为 A 的 n 阶行列式 hl 等于零，所以 A 不能进行 LU 分解.A 的秩 RA 如下:), RA,hl=det(A);m=0;returnendif RA=nfor p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif
8、 h(1,i)=0disp(请注意：因为 A 的 r 阶主子式等于零，所以 A 不能进行 LU 分解.A 的秩 RA 和各阶顺序主子式值 hl 依次如下： ), hl;RAm=0;returnendendif h(1,i)=0 m=1;disp(请注意：因为 A 的各阶主子式都不等于零，所以 A 能进行 LU 分解.A 的秩 RA 和各阶顺序主子式值 hl 依次如下：)for j=1:nU(1,j)=A(1,j);end5for k=2:nfor i=2:nfor j=2:nL(1,1)=1;L(i,i)=1;if ijL(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)
9、=A(i,1)/U(1,1);L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k)/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendhl;RA,U,Lendend运行结果：系数矩阵的阶数 3线性方程组的系数矩阵1,2,3;1,7,9;3,6,8线性方程组的常数项1,2,7请注意：因为 A 的各阶主子式都不等于零，所以 A 能进行 LU 分解.A 的秩 RA 和各阶顺序主子式值 hl 依次如下：RA =3U =1 2 30 7 60 0 -1L =1 0 01 1 03 0 1x =-22.0000 6.5714 -33.5714
