matlab电力系统快速解耦法潮流计算及短路计算程序

电力系统快速解耦法潮流分析及短路计算
一．程序设计的基本思想：
（1）由于电力系统潮流分析中要利用到矩阵运算，复数运算，故采用matlab编程。

采用文件输入，将系统的各个参数以文件的形式输入，便于程序的通用化。

（2）本程序共有两个输入文件，分别为线路参数的文件，和已知的节点状态文件（PQ）
（3）为了使程序不仅仅局限于计算9节点网络，在形成节点导纳的函数Yn（）中，利用循环，找出线路首节点中的最大编号数，自动确定节点导纳矩阵的维数。

故对于任意n节点网络，均可以计算出节点导纳矩阵
（4）在（3）的前提下，为了使程序支持系统增加节点，增加负荷等造成的PQ参数改变，或者PQ表的加长。

对程序做了如下优化。

首先，程序执行的基础是PQ表中平衡节点在第一行，接下来是PV节点，最后是PQ节点，如果系统添加节点，或者删除节点，均在PQ表的末端操作，会造成PQ表的顺序不是平衡节点、PV节点、PQ节点的顺序。

故引入了seqencing（）函数，其作用就是不论输入的PQ表是什么顺序，在程序读入后均按平衡-PV-PQ的顺序排列。

其次，顺序打乱的PQ表必须与支路参数表对应，故在Yn（）函数中加入了两段循环体，使之对应（见相应函数体注释）
（5）在满足了上述4个条件后，程序便可以通用化了。

当然，由于水平有限，且程序未能由大量数据测试，故缺陷在所难免，这里仅是做了通用化的尝试。

在本文最末附加了该程序通用化的实例。

二、潮流计算框图
三．定义相应的函数
1.形成节点导纳矩阵的函数Yn（）
function Y=Yn(x,y)
%定义一名为Yn的函数，其功能是自动识别输入表中节点的个数，形成相应的节点导纳矩阵
[fid,message]=fopen(x,'r') ; %从x文件中读入支路参数
if fid==-1; %判断文件是否正确打开
error(message);
end;
[HeadPoint,HeadNumber, EndPoint,EndNumber,R,X,B,k]=textread(x,'%s %d %s %d %f %f %f %f'); %将读入的参数处理为以列为向量的数组
fclose(fid);%关闭文件
L=length(HeadNumber); %确定输入表的行数
[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);
%调用seqencing函数，引入y文件中的PQ参数
A=PointNumber;
for i=1:L; %通过以下两循环体，实现PQ参数与支路参数的编号对应
for j=1:L;
if HeadNumber(i)==j;
HeadNumber(i)=A(j);
break;
end;
end;
end;
for i=1:L;
for j=1:L;
if EndNumber(i)==j;
EndNumber(i)=A(j);
break;
end;
end;
end;
Y=zeros(L,L); %根据txt文件中数据表的长度建立空的节点导纳矩阵
for i=1:L
m=HeadNumber(i);n=EndNumber(i);
if k(i)==0; %判断是否何种元件，为输电线元件
if n~=0;
Y(m,m)=Y(m,m)+1j*B(i)+1/(R(i)+1j*X(i));
Y(n,n)=Y(n,n)+1j*B(i)+1/(R(i)+1j*X(i));
Y(m,n)=Y(m,n)-1/(R(i)+1j*X(i));
Y(n,m)=Y(n,m)-1/(R(i)+1j*X(i));
else
Y(m,m)=Y(m,m)+R(i)+1j*X(i);
end;
else %为变压器元件
if n~=0;
Y(m,m)=Y(m,m)+1/(R(i)+1j*X(i));
Y(m,n)=Y(m,n)-1/(k(i)*(R(i)+1j*X(i)));
Y(n,n)=Y(n,n)+1/(k(i)*k(i)*(R(i)+1j*X(i)));
Y(n,m)=Y(n,m)-1/(k(i)*(R(i)+1j*X(i)));
else
Y(m,m)=Y(m,m)+R(i)+1j*X(i);
end;
end;
end;
maxm=HeadNumber(1);
%通过下面两个循环体，确定输入表中节点编号的最大值，及为节点导纳矩阵的维数
for i=1:L;
if maxm<=HeadNumber(i);
maxm=HeadNumber(i);
end;
end;
maxn=EndNumber(1);
for i=1:L;
if maxn<=EndNumber(i);
maxn=EndNumber(i);
end;
end;
Y=Y(1:max(maxm,maxn),1:max(maxm,maxn));%形成导纳矩阵
2．对不满足要求的PQ参数表进行排序的函数seqencing（）
function [Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y)
%定义名为seqencing的函数，其功能是在系统添加节点，或输入的PQ参数的顺序不满足要求时，对PQ参数表进行重新排序，保证平衡节点放在第一行，接下来是PV节点，最后是PQ节点
[fid,message]=fopen(y,'r'); %从y文件中读入PQ参数
if fid==-1; %判断文件是否正确打开
error(message);
end;
[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=textread(y,'%f %f %f %f %f %f');
fclose(fid);
L=length(PointNumber);
%通过以下两个循环体，完成对PQ输入表的重新排序，其思想是，在PQ参数之前加入一列Pointstyle用于标识节点类型，平衡节点为0,PV节点为1，PQ节点为2,以Pointstyle列为基准进行排序
for i=1:L;
for j=1:L-i;
if Pointstyle(j)>Pointstyle(j+1);
t=Pointstyle(j+1);
Pointstyle(j+1)=Pointstyle(j);
Pointstyle(j)=t;
t=PointNumber(j+1);
PointNumber(j+1)=PointNumber(j);
PointNumber(j)=t;
t=Ps(j+1);
Ps(j+1)=Ps(j);
Ps(j)=t;
t=Qs(j+1);
Qs(j+1)=Qs(j);
Qs(j)=t;
t=Uk(j+1);
Uk(j+1)=Uk(j);
Uk(j)=t;
end;
end;
end;
3、形成解耦算法B’矩阵的函数 formB1（）
function B1=formB1(x,y)
%定义名为B1的函数形成解耦算法中的B’矩阵，得到的B’矩阵用B1表示
[fid,message]=fopen(x,'r') ; %从x文件中读入支路参数
if fid==-1; %判断文件是否正确打开
error(message);
end;
[HeadPoint,HeadNumber,EndPoint,EndNumber,R,X,B,k]=textread(x,'%s %d %s %d %f %f %f %f'); %将读入的参数处理为以列为向量的数组
fclose(fid);
L=length(HeadNumber);
[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);
%调用seqencing函数，引入y文件中的PQ参数
A=PointNumber;
%通过以下两循环体，实现PQ参数与支路参数的编号对应
for i=1:L;
for j=1:L;
if HeadNumber(i)==j;
HeadNumber(i)=A(j);
break;
end;
end;
end;
for i=1:L;
for j=1:L;
if EndNumber(i)==j;
EndNumber(i)=A(j);
break;
end;
end;
end;
B1=zeros(L,L);
for i=1:L %以行为单位，通过循环，用支路参数对B1进行修改，形成B’矩阵 m=HeadNumber(i);n=EndNumber(i);
B1(m,m)=B1(m,m)-1/X(i);
B1(n,n)=B1(n,n)-1/X(i);
B1(m,n)=B1(m,n)+1/X(i);
B1(n,m)=B1(n,m)+1/X(i);
end
maxm=HeadNumber(1);
for i=1:L;
if maxm<=HeadNumber(i);
maxm=HeadNumber(i);
end;
end;
maxn=EndNumber(1);
for i=1:L;
if maxn<=EndNumber(i);
maxn=EndNumber(i);
end;
end;
B1=B1(2:max(maxm,maxn),2:max(maxm,maxn)); %形成B’矩阵
4、形成解耦算法B’’矩阵的函数 formB11（）
function B11=formB11(x,y)
%定义名为B11的函数形成解耦算法中B’'矩阵，用B11表示从x文件中读入支路参数确定Y，从y文件中读入PQ参数确定B11的维数，即除去平衡节点和pv节点，此处要求PQ参数录入时，将平衡节点和PQ节点放在前排，这一要求在Yn函数中通过seqencing函数已经满足
Y=Yn(x,y);
B=imag(Y);
[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);
i=1;j=1;
while Pointstyle(i)<=1;
i=i+1;
j=j+1;
end;
B11=B(j:end,j:end); %形成B’’矩阵
5、计算正常情况下系统节点电压的函数 powerflow（）
function [U0,O0]=powerflow(x,y)
%定义名为powerflow的函数，利用快速解耦算法来计算正常情况下系统内各个节点的电压和相角
[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);% 调用seqencing函数对PQ参数表进行排序
Y=Yn(x,y); %形成节点导纳矩阵，Yn为n维
B1=formB1(x,y); %形成解耦算法中的B矩阵，B1为n-1维
B11=formB11(x,y); %形成解耦算法中的B'矩阵，B'为m维
G=real(Y); %取Y的实部
B=imag(Y); %取Y的虚部
U0=Uk;
O0=Ok;
L=length(PointNumber);
P=zeros(L,1);
Q=zeros(L,1);
dP=zeros(L,1);
dQ=zeros(L,1);
number=1;
i=1;k=1;
while Pointstyle(i)<=1; %通过k值确定系统中PQ节点的个数
i=i+1;
k=k+1;
end;
while number<100 %定义迭代次数上限为100次
for i=2:L;
sum1=0;
for j=1:L;
sum1=sum1+U0(j)*(G(i,j)*cos(O0(i)-O0(j))+B(i,j)*sin(O0(i)-O0(j))); %潮流方程，n-1维
end;
dP(i)=Ps(i)-U0(i)*sum1;
end
for i=k:L;
sum2=0;
for j=1:L;
sum2=sum2+U0(j)*(G(i,j)*sin(O0(i)-O0(j))-B(i,j)*cos(O0(i)-O0(j))); %潮流方程，m 维
end;
dQ(i)=Qs(i)-U0(i)*sum2;
end
dP1=dP(2:L)./U0(2:L);
dQ1=dQ(k:L)./U0(k:L);
a=max(norm(dP1,inf));
b=max(norm(dQ1,inf));
if max(a,b)<0.00001 %判断是否收敛
break;
disp(‘迭代’)
disp(k);
disp(‘次后收敛’);
else %如不收敛，
dO=-inv(B1)*dP1; %dO为n-1维
dU=-inv(B11)*dQ1; %dU为m维
zero1=zeros(k-1,1);
zero2=[0];
DU=[zero1;dU];
DO=[zero2;dO];
U0=U0+DU;
O0=O0+DO;
number=number+1;
end;
if number==100;
disp('迭代100次后不收敛，迭代结束');
end;
end;
四．对相应的系统进行潮流分析和短路计算
定义完上述函数之后，可直接调用函数形成导纳矩阵，计算正常情况下的节点电压，进行短路计算计算短路电流，短路后各个节点电压以及支路潮流分布。

程序的输入表共有三个，为‘network.txt’，‘network2.txt’，‘PQ参数.txt’。

其中‘network.txt’存放不含发电机节点导纳和负荷导纳的系统支路参数，‘network2.txt’存放含发电机节点导纳和负荷导纳的系统支路参数，‘PQ参数.txt’存放系统的PQ参数，现列表如下
第一列：首节点名称第二列：首节点编号第三列：尾节点名称第四列：尾节点编号第五列：支路参数R 第六列：支路参数X 第七列：支路参数B/2 第八列：变比K
第一列：节点类型，0-平衡节点，1-PV节点，2-PQ节点第二列：节点编号
第三列：节点的有功P 第四列：节点的无功Q 第五列：节点电压，PQ节点置为0 第六列：节点电压的相角，初始值为0
程序代码如下：
%正常情况下节点导纳矩阵
Y=Yn(‘network.txt’,PQ参数.txt’);
disp(‘节点导纳矩阵’)
disp(sparse(Y));
%正常情况下潮流计算结果兵绘制潮流计算框图
[U0,O0]=powerflow('network.txt','PQ参数.txt') %各个节点电压为U0，相角为O0 disp('正常情况下系统各个节点电压’);
disp(U0);
disp(‘正常情况下系统各个节点相角’);
disp(O0);
O0=O0*180/pi
disp(‘正常情况下系统各个节点相角角度值’)
%计算平衡节点功率和PV节点无功功率：
for m=1:9
p1=p1+U0(1)*U0(m)*(G(1,m)*cos(O0(1)-O0(m))+B(1,m)*sin(O0(1)-O0(m)));
q1=q1+U0(1)*U0(m)*(G(1,m)*sin(O0(1)-O0(m))-B(1,m)*cos(O0(1)-O0(m)));
q2=q2+U0(2)*U0(m)*(G(2,m)*sin(O0(2)-O0(m))-B(2,m)*cos(O0(2)-O0(m)));
q3=q3+U0(3)*U0(m)*(G(3,m)*sin(O0(3)-O0(m))-B(3,m)*cos(O0(3)-O0(m)));
end;
disp(‘平衡节点有功、无功，两个电压节点有功分别为’)
disp(p1); %有功功率和无功功率
disp(q1); %平衡节点无功功率
disp(q2); %PV节点2无功功率
disp(q3); %PV节点3无功功率
%增加发电机导纳yi和负荷导纳yLDi后的发电机节点和负荷节点的自导纳
Y=Yn('network2.txt','PQ参数.txt') %network2文件中存放的是经修改的支路参数
disp(‘加发电机导纳yi和负荷导纳yLDi后的发电机节点和负荷节点的自导纳’);
diag(Y)
%精确法计算短路电流，短路后各节点电压以及网络中各支路的电流分布
Y=Yn('network2.txt','PQ参数.txt'); %形成包括发电机内阻抗和符合阻抗节点导纳矩阵 Z=inv(Y); %形成节点阻抗矩阵 %题中给出的是4节点短路
Zf=Z(:,4);
[U0,O0]=powerflow('network.txt','PQ参数.txt');
U01=U0.*cos(O0)+1j*U0.*sin(O0);
If=U01(4)/Z(4,4); %计算短路电流If
U=U01-Zf*If; %计算短路后的各个节点电压
[fid,message]=fopen('network2.txt','r') ; %从x文件中读入支路参数
if fid==-1; %判断文件是否正确打开
error(message);
end;
[HeadPoint,HeadNumber,EndPoint,EndNumber,R,X,B,k]=textread('network2.txt','%s %d %s %d %f %f %f %f');
row=[HeadNumber,EndNumber];
for i=1:9 %计算系统中各个支路的潮流分布
I(row(i,1),row(i,2))=(U(row(i,1))-U(row(i,2)))/Z(row(i,1),row(i,2));
end;
disp(‘阻抗矩阵中的第f列’);
disp(‘精确算法短路电流If’);
disp(‘If模值’);
disp(abs(If));
disp(‘相角为’);
disp(angle(If));
disp(‘短路后各个节点的电压’);
disp(U);
disp(‘模值为’);
disp(‘短路后各个支路的电流值 I‘) ;
disp(‘模值为 I’);
disp(‘相角为（度）’);
%近似计算短路电流，短路后各节点电压以及网络中各支路的电流分布
If2=1/Z(4,4); %计算短路电流If
U2=1-Zf*If; %计算短路后的各个节点电压
disp(‘近似算法短路电流 If2’);
disp(‘If2模值’);
disp(abs(If));
disp(‘相角为（度）’);
disp(angle(If2));
disp(‘短路后各个节点的电压’);
disp(U2);
disp(‘模值为’);
程序运行结果如下：
节点导纳矩阵
Y =
(1,1) 0 -17.3611i
(4,1) 0 +17.3611i
(2,2) 0 -16.0000i
(7,2) 0 +16.0000i
(3,3) 0 -17.0648i
(9,3) 0 +17.0648i
(1,4) 0 +17.3611i
(4,4) 3.3074 -39.3089i
(5,4) -1.3652 +11.6041i
(6,4) -1.9422 +10.5107i
(4,5) -1.3652 +11.6041i
(5,5) 2.5528 -17.3382i
(7,5) -1.1876 + 5.9751i
(4,6) -1.9422 +10.5107i
(6,6) 3.2242 -15.8409i
(9,6) -1.2820 + 5.5882i
(2,7) 0 +16.0000i
(5,7) -1.1876 + 5.9751i
(7,7) 2.8047 -35.4456i
(8,7) -1.6171 +13.6980i
(7,8) -1.6171 +13.6980i
(8,8) 2.7722 -23.3032i
(9,8) -1.1551 + 9.7843i
(3,9) 0 +17.0648i
(6,9) -1.2820 + 5.5882i
(8,9) -1.1551 + 9.7843i (9,9) 2.4371 -32.1539i
迭代 k=9
次后收敛
正常情况下系统各个节点电压
U0 =
1.0400
1.0250
1.0250
1.0258
0.9956
1.0127
1.0258
1.0159
1.0324
正常情况下系统各个节点相角
O0 =
0.1620
0.0814
-0.0387
-0.0696
-0.0644
0.0649
0.0127
0.0343
正常情况下系统各个节点相角角度值
O0 =
9.2800
4.6647
-2.2168
-3.9888
-3.6875
3.7197
0.7275
1.9667
平衡节点有功、无功，两个电压节点有功分别为 p1=
0.7164
q1=
0.2705
q2=
0.0665
q3=
-0.1086
加发电机导纳yi和负荷导纳yLDi后的发电机节点和负荷节点的自导纳
ans =
0 -20.6944i
0 -19.3333i
0 -20.3982i
3.3074 -39.3089i
1.2917 -16.8338i
2.3466 -15.5484i
2.8047 -35.4456i
1.8033 -2
2.9641i
2.4371-32.1539i
阻抗矩阵中的第f列
Zf =
-0.0630 + 0.1475i
-0.0645 + 0.0860i
-0.0634 + 0.0864i
-0.0751 + 0.1758i
-0.0892 + 0.1534i
-0.0854 + 0.1532i
-0.0779 + 0.1040i
-0.0832 + 0.1032i
-0.0758 + 0.1032i
精确算法短路电流
If =
-2.2963 - 4.8486i
If模值为
5.3649
相角为
-115.3422
短路后各个节点的电压
U =
0.1801 + 0.0333i
0.4463 + 0.0502i
0.4573 - 0.0257i
-0.0000 - 0.0000i
0.0446 - 0.1498i
0.0716 - 0.1271i
0.3406 - 0.0725i
0.3243 - 0.1533i
0.3572 - 0.0950i
模值为
0.1831
0.4491
0.4580
0.0000
0.1563
0.1459
0.3482
0.3587
0.3696
短路后各个支路的电流值 I
(1,4) -0.2501 - 1.1140i (4,5) 0.8560 - 0.2072i (4,6) 0.8320 + 0.0041i (2,7) 0.5086 - 0.9694i (5,7) 0.6818 + 1.8688i (7,8) 0.3965 - 0.3175i (3,9) 0.2141 - 0.7887i (6,9) 0.8241 + 1.7182i (8,9) -0.2355 + 0.3676i
模值为 I
(1,4) 1.1417
(4,5) 0.8807
(4,6) 0.8320
(2,7) 1.0947
(5,7) 1.9893
(7,8) 0.5079
(3,9) 0.8172
(6,9) 1.9056
(8,9) 0.4365
相角为（度）
(1,4) -102.6530
(4,5) -13.6077
(4,6) 0.2855
(2,7) -62.3148
(5,7) 69.9555
(7,8) -38.6844
(3,9) -74.8141
(6,9) 64.3752
(8,9) 122.6426
近似算法短路电流
If2 =
0.2105 - 5.2357i
If2模值 5.2399
相角为（度） -87.6981
短路后各个节点电压值
U2 =
0.1611 + 0.0000i
0.4785 - 0.0351i
0.4773 - 0.0384i
0.0000
0.1155 - 0.0063i
0.1178 - 0.0098i
0.3698 - 0.0424i
0.3673 - 0.0436i
0.3752 - 0.0459i
模值为
0.1611
0.4797
0.4788
0.0000
0.1157
0.1182
0.3722
0.3698
0.3780
精确算法与近似算法的比较分析（1）短路电流误差
（If精-If近）/If精=2.33%
\φ精-φ近\=27.6441（度）
两电流相差进2.33%，可见在精度要求不是很高的前提下，近似算法是完全可取的。

（2）电压模值误差比较
（U精-U近）/U精=
0.1205
-0.0681
-0.0455
0.0000
0.2597
0.1900
-0.0690
-0.0311
-0.0227
可见电压模值的误差均在25%以内。

至此，上述程序已完成题目的所有要求。

为了实现程序的通用行，在编写程序的时候对程序进行了优化，使得上述程序不但能够满足本题要求，而且对于任意多节点的系统，只要将支路参数，PQ参数列表输入txt 文件中，有上述程序均可完成正常情况下潮流计算和短路计算。

而且，上述程序还支持支路参数的修改，添加或删除支路、节点。

四．附加：程序通用化实例
若在原题的基础上，在8母线和9母线之间增加一条母线10,，相应的支路参数如下表
7~8
0.0085 0.072
0.0745
8~10 0.0059 0.0508 0.1045 9~10 0.0050 0.0050 0.1045
2
9
83
6
4
1
5
(3)
f
7
新增10节点为PQ节点 P=0，Q=0；
由上述条件可以写出添加节点后的支路参数文件为PQ参数文件为
10
以这两个表作为输入，执行上述程序，计算潮流，得到各个节点的电压和相角为
U0 =
1.0400
1.0250
1.0250
1.0278
0.9988
1.0168
1.0304
1.0265
1.0397
1.0385
O0 =
0.1601
0.0795
-0.0386
-0.0695
-0.0643
0.0635
0.0113
0.0327
0.0215
即可得到十个节点的电压和相角，完成了程序的通用化。

五．结语
学习是一个理论与实践相结合的过程，学习电力系统分析，对相应的电力系统进行简单的潮流计算和短路分析是我们必须掌握的基本知识。

在用计算机编程是手段实现系统分析的功能时，我们的问题是只着眼于既定的题目和系统。

所以，怎样将程序通用化是一个值得考虑的问题，由于电力系统的复杂性以及个人知
识水平的限制，谨在此方面做了一些改进，真正的知识体系还亟待在今后的学习中得以建立和提高。

附录资料：matlab画二次曲面
一、螺旋线
a=0:0.1:20*pi;
h=plot3(a.*cos(a),a.*sin(a),2.*a,'b','linewidth',2); axis([-50,50,-50,50,0,150]);
grid on
set(h,'erasemode','none','markersize',22);
xlabel('x轴');ylabel('y轴');zlabel('z轴');
title('静态螺旋线');
2.动态螺旋线
t=0:0.1:10*pi;
i=1;
h=plot3(sin(t(i)),cos(t(i)),t(i),'*','erasemode','none'); grid on
axis([-2 2 -2 2 0 35])
for i=2:length(t)
set(h,'xdata',sin(t(i)),'ydata',cos(t(i)),'zdata',t(i));
drawnow
pause(0.01)
end
title('动态螺旋线');
(图略)
t=0:0.1:10*pi;
x=r.*cos(t);
y=r.*sin(t);
z=t;
plot3(x,y,z,'h','linewidth',2);
grid on
axis('square')
xlabel('x轴');ylabel('y轴');zlabel('z轴'); title('圆柱螺旋线')
二、旋转抛物面
b=0:0.2:2*pi;
[X,Y]=meshgrid(-6:0.1:6);
Z=(X.^2+Y.^2)./4;
meshc(X,Y,Z);
axis('square')
xlabel('x轴');ylabel('y轴');zlabel('z轴'); shading flat;
title('旋转抛物面')
或直接用：ezsurfc('(X.^2+Y.^2)./4')
三、椭圆柱面
load clown
ezsurf('(2*cos(u))','4*sin(u)','v',[0,2*pi,0,2*pi]) view(-105,40) %视角处理
shading interp %灯光处理
colormap(map) %颜色处理
grid on %添加网格线
axis equal %使x,y轴比例一致xlabel('x轴');ylabel('y轴');zlabel('z轴'); shading flat;
title('椭圆柱面') %添加标题
四、椭圆抛物面
b=0:0.2:2*pi;
[X,Y]=meshgrid(-6:0.1:6);
Z=X.^2./9+Y.^2./4;
meshc(X,Y,Z);
axis('square')
xlabel('x轴');ylabel('y轴');zlabel('z轴'); shading flat;
title('椭圆抛物面')
或直接用：ezsurfc('X.^2./9+Y.^2./4')
五、'双叶双曲面
ezsurf('8*tan(u)*cos(v)','8.*tan(u)*sin(v)','2.*sec(u)',[-pi./2,3*pi./2,0,2*pi]) axis equal
grid on
axis square
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shading flat;
title('双叶双曲面')
六、双曲柱面
load clown
ezsurf('2*sec(u)','2*tan(u)','v',[-pi/2,pi/2,-3*pi,3*pi]) hold on %在原来的图上继续作图
ezsurf('2*sec(u)','2*tan(u)','v',[pi/2,3*pi/2,-3*pi,3*pi]) colormap(map)
shading interp
view(-15,30)
axis equal
grid on
axis equal
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shading flat;
title('双曲柱面')
七、双曲抛物面（马鞍面）
[X,Y]=meshgrid(-7:0.1:7);
Z=X.^2./8-Y.^2./6;
meshc(X,Y,Z);
view(85,20)
axis('square')
xlabel('x轴');ylabel('y轴');zlabel('z轴'); shading flat;
title('双曲抛物面')
或直接用：ezsurfc('X.^2./8-Y.^2./6')
八、抛物柱面
[X,Y]=meshgrid(-7:0.1:7);
Z=Y.^2./8;
h=mesh(Z);
rotate(h,[1 0 1],180) %旋转处理%axis([-8,8,-8,8,-2,6]);
axis('square')
xlabel('x轴');ylabel('y轴');zlabel('z轴'); shading flat;
title('抛物柱面')
或直接用：ezsurfc('Y.^2./8')
九、环面
ezmesh('(5+2*cos(u))*cos(v)','(5+2*cos(u))*sin(v)','2*sin(u)',[0,2*pi,0,2*pi]) axis equal
grid on
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shading flat;
title('环面')
十、椭球
ezsurfc('(5*cos(u))*sin(v)','(3*sin(u))*sin(v)','4*cos(v)',[0,2*pi,0,2*pi]) axis equal
grid on
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shading flat;
title('椭球')
十一、单叶双曲面
ezsurf('4*sec(u)*cos(v)','2.*sec(u)*sin(v)','3.*tan(u)',[-pi./2,pi./2,0,2*pi]) axis equal
grid on
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shading flat;
title('单叶双曲面')
十二、旋转单叶双曲面
load clown
ezsurf('8*sec(u)*cos(v)','8.*sec(u)*sin(v)','2.*tan(u)',[-pi./2,pi./2,0,2*pi]) colormap(map)
view(-175,30)
%alpha(.2) %透明处理
axis equal
grid on
axis square
xlabel('x轴');ylabel('y轴');zlabel('z轴');
shading flat;
title('旋转单叶双曲面')
十三、圆柱面
subplot(1,2,1)
ezsurf('(2*cos(u))','2*sin(u)','v',[0,2*pi,0,2*pi]) grid on
shading interp
axis equal
xlabel('x轴');ylabel('y轴');zlabel('z轴');
title('圆柱面')
subplot(1,2,2)
cylinder(30)
shading interp
axis square
title('调用cylinder函数所得圆柱面')
十四、二次锥面
clc,clear;
P=[1,0,0;
0,cos(45*pi/180),sin(45*pi/180);
0,-sin(45*pi/180),cos(45*pi/180)];
for k2 = 1:31
for k1 = 1:31
x(k1,k2) = (k2-1)*cos ( (k1-1)*12*pi/180);
y(k1,k2) = (k2-1)*sin ( (k1-1)*12*pi/180);
z(k1,k2) = sqrt(x(k1,k2)^2+y(k1,k2)^2);
Vxyz = P*[x(k1,k2),y(k1,k2),z(k1,k2)]';
x1(k1,k2)=Vxyz(1);
y1(k1,k2)=Vxyz(2);
z1(k1,k2)=Vxyz(3);
end
end
surf(x,y,z)
hold on;
surf(x1,y1,z1);
shading flat;。