MATLAB 期末复习
1．在MATLAB 帮助文件中查找有关title 的使用方法，并为y1对应的图形添加标题 y1=sin(t) ；为y2对应图形添加标题 y2=e -at cos(3t)，其中a 根据输入显示具体值（使用num2str 函数）。请写出有关指令。
答：title('y1=sin(t)') title(['y2=exp^{-',num2str(a),'t},cos(3t)'])
3. 运行以下指令
x=-3*pi:pi/15:3*pi;
y=x;
[X,Y]=meshgrid(x,y); warning off;
Z=sin(X).*sin(Y)./X./Y;
1）继续编写指令求出矩阵Z 中有多少个“非数”数据。
答： 指令：m=sum(sum(isnan(Z)))
运行结果： m =181
4）写出绘制相应“无裂缝”图形的全部指令。(参考课本p147)
答：x=-3*pi:pi/15:3*pi;
x=x+(x==0)*eps;%用逻辑数组参与计算，使0元素被一个“机器零”小数替代
y=x; [X,Y]=meshgrid(x,y); warning off;
Z=sin(X).*sin(Y)./X./Y; sum(sum(isnan(Z))) %统计Z 中非数的数目，数组要用两个sum surf(X,Y,Z); shading interp xlabel('x') ylabel('y') zlabel('z')
4.运行clear, rand('state',1),A=rand(3,3)
1）编写指令，获得一个对角矩阵B ，其相应元素由A 的对角元素构成
答：B=diag(A)
2）编写指令，获得一个矩阵C ，其对角元素全部为零，其余元素与A 对应元素相同。（使用diag 函数）
答：C=diag(diag(A));
5. 运行clear, rand('state',1),A=rand(3,5)生成二维数组A ，编写指令对A 矩阵进行以下操作，记录相应的结果。
1）对A 矩阵旋转180度，保存为B ；
答： B=rot90(A,2)
2）将A 矩阵排成5行3列，保存为C ；
答： C=reshape(A,5,3)
3）将A 矩阵删除第2及4列，保存为D ；
答： D=A;
D(:,[2,4])=[]
4）求该数组中所有大于0.5的元素的位置，分别求出它们的“全下标”
和“单下标”。（使用find 和sub2ind 函数）
答： n=(A>0.5) ; [m,n]=find(A>0.5; L=sub2ind(size(A),m,n)
6.求以下方程组的解1） ???????=--+=-+=-++=++184/102/2/d n p q p n d q p q d n q p n d 2）???????-=--+=--+=-++=-++-5
.385.04/31d n p q p n d q p q d n q p n d 答：（参考课本p173)
clear;
A=[1 -1 1/2 1/2;1 1 1 -1;1 1 -1/4 1;-8 1 -1 1]; b=[0 10 0 1]';
rank(A) %求A 的秩
rank([A,b]) %求增广矩阵的秩
B=A\b;
clear;
A=[-1 1 1 -1;1 1 -1 1;1 -1/4 -1 1;-8 -1 1 1]; b=[1 3 0.5 -3.5]';
rank(A) rank([A,b])
B=A\b;
C=null(A); A*(rand(1)*C+B) b
7.给定一个函数2
)(t e t f -=，利用diff 及gradient 函数近似求)(t f 的导函数)()(t f dt
t df '=，并将函数)(t f 和其导函数曲线绘制在一个图中，同时给出两种方法得到的导函数曲线的细节图，并进行对比。
答：clear;t=0:0.1:10; dif=diff(y)/0.1;
grad=gradient(y)/0.1;
subplot(1,2,1); plot(t(1:end-1),dif,t,grad);
legend('dfdt_{diff}','dfdt_{gradient}') %对图形按指令顺序进行标注
xlabel('t'); subplot(1,2,2);
plot(t(end-10:end-1),dif(end-9:end),'.k',t(end-9:end),grad(end-9:end),'or')
legend('dfdt_{diff}','dfdt_{gradient}'); xlabel('t'); 8.利用sum 、trapz 、及quad 函数求?=50)(dt t f S ，其中2)(t e t f -=。记录结果，
要求quad 的计算结果达到绝对精度为10-9。观察各函数计算结果之间的差异。 答：clear;format long; t=0:0.1:5; y=exp(-t.^2); s1=sum(y)*0.1 s2=trapz(t,y)
quad('exp(-t.^2)',0,5,1e-9) %1e-9代表计算结果的绝对精度为10-9
9.利用函数cumsum 和cumtrapz 求?=x dt t f x y 0)()(，其中2
)(t e t f -=。将两者计算结果及它们的误差绘制在同一个图中，同时给出误差的细节图。
在以上题目中要求在图中为曲线添加相应的图例，原函数曲线用细线，其他曲线用粗线。
答：clear;format long; t=0:0.1:5; y=exp(-t.^2);
s1=cumsum(y)*0.1; s2=cumtrapz(t,y);
subplot(1,2,1); hold on; plot(t,y)
plot(t,s1,'-k','LineWidth',2); %绘制t-s1关系图像，线条大小为2号
plot(t,s2,'-r','LineWidth',2);
plot(t,s1-s2,'-g','LineWidth',2);
xlabel('t');
legend('原函数','cumsum 积分','cumtrapz 积分','误差');
subplot(1,2,2);
plot(t,s1-s2,'LineWidth',2);
xlabel('t');
axis([0,5,0.04,0.11])
10.求函数5.08.1)2cos(5.1))5(sin()(2
06.02++-=t t t e t t f t 在区间[-5,5]中的最小值点及相应值。
答：clear; t=-5:0.1:5;
f=(sin(5*t)).^2.*exp(0.06*t.^2)-1.5*t.*cos(2*t)+1.8*(abs(t+0.5));
plot(t,f)
f1=@(t)((sin(5*t)).^2.*exp(0.06*t.^2)-1.5*t.*cos(2*t)+1.8*(abs(t+0.5)));
[sx,sfval,sexit,soutput]=fminbnd(f1,-5,5) 11.求0))sin(sin(105.02.0=-+--t e t t 的实数解。
答：clear;
t=-5:0.1:5;
f=-0.5+t-10*exp(-0.2*t).*abs(sin(sin(t)));
plot(t,f)
f1=@(t)(-0.5+t-10*exp(-0.2*t).*abs(sin(sin(t))));
%[x,fval]=fsolve(f1,1) %没有找到
[x,fval]=fsolve(f1,2) %找到
12.使用字符串、内联对象、匿名函数或函数文件的句柄四种方式求二元函数方
程组?
??=+=-0)cos(0)sin(y x y x 的数值解 答：clear all
x0=[0.5 0.5]
[x,fval]=fsolve('[sin(x(1)-x(2)),cos(x(1)+x(2))]',x0)
字符串解法
S=solve('sin(x-y)=0','cos(x+y)=0','x','y')
S.x, S.y
内联对象
f=inline('[sin(x(1)-x(2)),cos(x(1)+x(2))]');
[x,fval]=fsolve(f,x0)
匿名函数
f=@(x)([sin(x(1)-x(2)),cos(x(1)+x(2))]);
[x,fval]=fsolve(f,x0)
函数文件的句柄
（1）将以下指令创建M 文件，命名为Myfsove 保存在matlab 搜索路径下
function F=Myfsove(x)
F=[sin(x(1)-x(2));cos((x(1)+x(2)))];
（2）在matlab 指令窗口运行以下程序
[x,fval]=fsolve(@Myfsove,x0)
13. 产生均值为4，标准差为2的（10000*1）正态分布随机数组a ，分别用hist 和histfit 绘制该数组的频数直方图，观察两张图的差异。除histfit 上的拟合红线外，如何使这两个指令绘出相同的频率直方图。
答：clear;
N=10000;Mu=4;Sigma=2;
%x=randn(N,1)
%y=Mu+Sigma*x;
y=normrnd(Mu,Sigma,1,N);
figure();
hist(y,100);
xlabel('x')
ylabel('P')
figure();
histfit(y);
xlabel('x')
ylabel('P')
14.根据文件prob_data418.mat 的数据，求出这组数据的5阶拟合多项式。根据求得的多项式函数计算y 值，并与实测y 值比较，统计两者之间误差的均值及方差。 答：s=importdata('prob_data418.mat')
x=s.x;
y=s.y;
P=polyfit(x,y,5)
yy=polyval(P,x);
plot(x,y,'.r',x,yy,'k')
xlabel('x')
ylabel('y')
err=abs(y-yy);
em=mean(err)
ev=var(err)
15.求有理分式)
()()(x D x N x R =的商多项式)(x Q 和余多项式)(x r ，其中，)5.0)(3()(33++=x x x x N ，)125)(22()(232++-+=x x x x x D 。
答：N=conv([3 0 1 0],[1 0 0 0.5]);D=conv([1 2 -1],[5 2 0 1]);
[Q r]=deconv(N,D)
16.已知函数[]
111021)sin(1)()()(>=<≤?????????-=+-----ξξξββββξβξξt
t t
t e e te t e t y ，其中21ξβ-=。试在同一张图上绘出4.1:2.0:2.0=ξ不同取值时的各条曲线，时间区间为]18,0[∈t 。要求1<ξ的各条曲线用蓝色，1>ξ的用红色，1=ξ的用黑色，且线较粗。2.0=ξ及4.1=ξ的两条曲线添加text 标志。
答：clear;
t=0:0.01:18;
figure(1);hold on;
for z=0.2:0.2:1.4
b=sqrt(abs(1-z^2));
if (z<1)&(z>=0)
y=1/b*exp(-z*t).*sin(b*t);
plot(t,y,'-b');
elseif z>1
y=1/2/b*(exp(-(z-b)*t)-exp(-(z+b)*t));
plot(t,y,'-r');
elseif z==1
y=t.*exp(-t);
plot(t,y,'k','LineWidth',3);
end
end
text(2.5,0.6,'\zeta=0.2')
text(0.5,0.15,'\zeta=1.4')
text(10,0.7,'\Delta\zeta=0.2')
grid on;
box on;
axis([0 18 -0.4 0.8])
xlabel('x')
ylabel('y')
17.使用极坐标绘图指令polar 绘制表达式θρcos 1-=对应的心形图。
答：the=0:pi/100:2*pi;
P=1-cos(the);
h=polar(the,P);%返回线的句柄
set(h,'LineWidth',3)
title('P=1-cos(\theta)')
18.在区域]3,3[,-∈y x ，绘制224y x xe z --=的三维（透视）网格曲面，如图所示。要求，透明度设置为0.5，坐标轴范围与图示一致。
答：x=-3:0.1:3; y=x;
[X,Y]=meshgrid(x,y);
Z=4*X.*exp(-X.^2-Y .^2);
surfc(X,Y ,Z); %surfc 为绘制带等高线的三维表面图
colormap(cool);
box on;
%xlim([-3 3]);
%ylim([-3 3]);
axis([-3 3 -3 3 -2 2]); %将xyz 坐标的范围限制在（-3 3 -3 3 -2 2）
set(gca,'FontSize',14,'FontAngle','italic','XTick',-3:3);%字体设置对后面的label 起作用
xlabel('x'); ylabel('y'); zlabel('z')
alpha(0.5) %透明度设置为0.5
19.用指令fminbnd求)
y x-
=在0
x
-
(x
sin(cos
)
e
x附近的最小值。要求分别写出
=
fminbnd的第一个输入量使用匿名函数及函数文件句柄的相关程序，观察两种方式所得结果是否一致。
答：结果一致
>> %匿名函数法
y0=@(x)(-exp(-x)*abs(sin(cos(x)))); fplot(y0,[-6,4])
[x,fval]=fminbnd(y0,-1,1)
hold on, plot(x,fval,'*r')
%函数句柄法
hh=@ff; ezplot(hh,[-6,4]) , [x,fval]=fminbnd(hh,-1,1); hold on, plot(x,fval,'*r')
x =
-0.8634
fval =
-1.4348
20..创建一个名为showfig的函数文件，该函数文件根据输入宗量不同情况显示出不同的图形。其中绘图部分由函数文件中定义的子函数figplot完成。具体要求为：没有输入量时，画出单位圆；输入量为大于2的自然数N时，绘制正N 边形，图名反映多边形的边数；输入其它值时，提示出错信息“输入错误，请输入多边形边长数目”。新建一个M文件，在M文件中调用函数文件showfig，测试无输入、输入为2.5及输入为6的图形或提示结果。
答：
%showfig函数文件注：该题解法绘图部分没有定义子函数figplot完成
function showfig(n)
if nargin~=0 & (n<0 | fix(n)~=n)
disp('Input parameter error! Please input natural number n larger than 2.') else
switch nargin
case 0
m=100;t=(0:m)/m*2*pi;
x=sin(t);y=cos(t);
plot(x,y,'r','LineWidth',2.5);title('Circle')
axis equal, axis off
otherwise
if n==1 | n==2,disp('please input n>=3')
else
t=(0:n)/n*2*pi;
x=sin(t);y=cos(t);
plot(x,y,'r','LineWidth',2.5); title(['Poly gon nit',int2str(n),' edges'])
axis equal, axis off
end
end
end
%M 文件调用showfig
figure
showfig()
figure
showfig(2.5)
figure
showfig(6)
21.在matlab 的\toolbox\matlab\elmat\private 文件夹下有一个“烟圈矩阵”发生函数smoke.m 。运行指令smoke(3,0)，将生成一个三阶伪特征根矩阵：
????
??????--+-=0000.100000.10000.18660.05000.0000000.18660.05000.0i i A 在当前目录为“E:\MATLAB 实验”时的指令窗口运行smoke(3,0)，提示找不到该函数。利用函数句柄可实现在当前目录下正常调用smoke.m 函数，产生3阶伪特征根矩阵。请写出相应的操作步骤。
答：%在当前目录设置为X ：\matlab\matlab\work 时的程序 X 代表安装matlab 程序所在的分区盘 clear
b=pwd %获取当前目录名字符串
which('smoke') %检查在当前目录下能否看到smoke.m
b_d=b
b_d(end-4:end)=[]; %在b 字符串中去除最后的四个字符，即\work 。
str=[b_d,'\toolbox\matlab\elmat\private'];
cd(str) %把smoke.m 所在目录设置为当前目
SM=@smoke; %创建smoke.m 的函数句柄
cd(b) %使work 恢复为当前目录
disp(' ')
A=feval(SM,3,0) %feval 为函数宏指令
disp(' ')
pwd
which('smoke')
22.写出使用for 和while 循环语句计算∑==1000
02.0i i K 的程序，同时给出避免循环的
数值计算程序。
答：>> %for 循环程序
k1=0.2;sum1=1;
for i=0:1000000
sum1=sum1+k1;k1=k1*0.2;
end
sum1
%while 循环程序
k2=0.2;sum2=1; i=0;
while i<=1000000
sum2=sum2+k2; k2=k2*0.2; i=i+1;
end
sum2
%数值计算法
N=1000000; S=zeros(1,N);
k=0:N;
S=0.2.^k;
sum(S)
23. 设1)()()(2)(2233=+-+t y dt t dy dt t y d dt t y d ，t 范围为[10，20]，初始条件为0)10(=y ，1)10(=dt dy ，0)10(22=dt y d ，采用数值计算法求解)(t y 。绘制)(t y 的曲线，绘出)(t y 与)(t y '之间关系的相平面图，并在图上添加适当的标题。
答：%函数文件DyDt.m
function ydot=DyDt1(t,y)
ydot=[y(2);y(3);-2*y(3)+y(2)-y(1)+1];
%运行指令
tspan=[10,20];
y0=[0;1;0];
[tt,yy]=ode45(@DyDt1,tspan,y0);
plot(tt,yy(:,1)) %第1列为求解结果x
xlabel('t'),title('y(t)')
figure();
plot(yy(:,1),yy(:,2)) %相轨迹
title('y —dy/dt')
课本参考例题
p129 p143 p147 p167 p179 p213（可视化中标识指令字符的精细控制） P217
华中科技大学 信息工程学院课程设计报告书题目: 电力系统潮流计算 专业：电气工程及其自动化 班级： 学号： 学生姓名： 指导教师： 2015年 11 月 10 日
2015年11月12日
信息工程学院课程设计成绩评定表
摘要 电力系统稳态分析包括潮流计算和静态安全分析。本文主要运用的事潮流计算，潮流计算是电力网络设计与运行中最基本的运算，对电力网络的各种设计方案及各种运行方式进行潮流计算，可以得到各种电网各节点的电压，并求得网络的潮流及网络中的各元件的电力损耗，进而求得电能损耗。本位就是运用潮流计算具体分析，并有MATLAB仿真。 关键词：电力系统潮流计算 MATLAB仿真
Abstract Electric power system steady flow calculation and analysis of the static safety analysis. This paper, by means of the calculation, flow calculation is the trend of the power network design and operation of the most basic operations of electric power network, various design scheme and the operation ways to tide computation, can get all kinds of each node of the power grid voltage and seek the trend of the network and the network of the components of the power loss, and getting electric power. The standard is to use the power flow calculation and analysis, the specific have MATLAB simulation. Key words: Power system; Flow calculation; MATLAB simulation
潮流例题：根据给定的参数或工程具体要求（如图），收集和查阅资料；学习相关软件（软件自选：本设计选择Matlab进行设计）。 2.在给定的电力网络上画出等值电路图。 3.运用计算机进行潮流计算。 4.编写设计说明书。 一、设计原理 1．牛顿-拉夫逊原理 牛顿迭代法是取x0 之后，在这个基础上，找到比x0 更接近的方程的跟，一步一步迭代，从而找到更接近方程根的近似跟。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新
的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。 牛顿—拉夫逊迭代法的一般步骤： （1）形成各节点导纳矩阵Y。 （2）设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。 （3）计算各个节点的功率不平衡量。 （4）根据收敛条件判断是否满足，若不满足则向下进行。 （5）计算雅可比矩阵中的各元素。 （6）修正方程式个节点电压 （7）利用新值自第（3）步开始进入下一次迭代，直至达到精度退出循环。 （8）计算平衡节点输出功率和各线路功率 2．网络节点的优化 １）静态地按最少出线支路数编号 这种方法由称为静态优化法。在编号以前。首先统计电力网络个节点的出线支路数，然后，按出线支路数有少到多的节点顺序编号。当由n 个节点的出线支路相同时，则可以按任意次序对这n 个节点进行编号。这种编号方法的根据是导纳矩阵中，出线支路数最少的节点所对应的行中非零元素也２）动态地按增加出线支路数最少编号在上述的方法中，各节点的出线支路数是按原始网络统计出来的，在编号过程中认为固定不变的，事实上，在节点消去过程中，每消去一个节点以后，与该节点相连的各节点的出线支路数将发生变化(增加，减少或保持不变)。因此，如果每消去一个节点后，立即修正尚未编号节点的出线支路数，然后选其中支路数最少的一个节点进行编号，就可以预期得到更好的效果，动态按最少出线支路数编号方法的特点就是按出线最少原则编号时考虑了消去过程中各节点出线支路数目的变动情况。 3．MATLAB编程应用 Matlab 是“Matrix Laboratory”的缩写，主要包括：一般数值分析，矩阵运算、数字信号处理、建模、系统控制、优化和图形显示等应用程序。由于使用Matlab 编程运算与人进行科学计算的思路和表达方式完全一致，所以不像学习高级语言那样难于掌握，而且编程效率和计算效率极高，还可在计算机上直接输出结果和精美的图形拷贝，所以它的确为一高效的科研助手。 二、设计内容 1.设计流程图
2．一单位反馈控制系统，若要求： （1）跟踪单位斜坡输入时系统的稳态误差为2。 （2）设该系统为三阶，其中一对复数闭环极点为-1±j1 求满足上述要求的开环传递函数。 利用MATALAB 完成以下操作： （1）求控制系统的单位阶跃响应、单位斜坡响应、单位加速度响应； （2）绘制控制系统的根轨迹，并判断系统的稳定性； （3）绘制控制系统的奈奎斯特图和Bode 图。 解：由于控制系统为单位反馈系统 则 H （s ）=1 由条件一可知，输入信号R(s)= 21s 时，稳态误差为： 2) (10)()(1)(==→+=s sG s s H s G s sR Ess 由条件二可知，系统为三阶且有。两个极点s1=-1+j1 ；s2=-1-j1 则可得，开环传递函数为 ) 222(1)(++=s s s s G （1）求控制系统的单位阶跃响应： G=tf([1],[1, 2, 2, 0]); Step(G) title('单位阶跃响应'); grid
(2)控制系统的单位斜坡响应：z0=0; p0=[0, -1+j,-1-j]; k0=1; [num0,den0]=zp2tf(z0,p0,k0); g0=tf(num0,den0) g=feedback(g0,1); num=g.num{1}; den=[g.den{1},0]; sys1=tf(num,den) step(sys1)
(3)控制系统的单位加速度响应：z0=0; p0=[0, -1+j,-1-j]; k0=1; [num0,den0]=zp2tf(z0,p0,k0); g0=tf(num0,den0) g=feedback(g0,1); num=g.num{1}; den=[g.den{1} 0 0]; sys2=tf(num,den) step(sys2)
第1章MATLAB概论 1.1与其他计算机语言相比较，MATLAB语言突出的特点是什么？ 答：MATLAB语言有一下特点： 1、起点高。 2、人机界面适合科技人员。 3、强大而简易的作图功能。 4、智能化程度高。 5、功能丰富，可扩展性强。 1.2 MATLAB系统由那些部分组成？ 答：MATLAB系统主要由开发环境、MATLAB数学函数库、MMATLAB语言、图形功能和应用程序接口五个部分组成。 1.3 安装MATLAB时，在选择组件窗口中哪些部分必须勾选，没有勾选的部分以后如何补安装？ 答：在安装MATLAB时，安装内容由选择组件窗口中各复选框是否被勾选来决定，可以根据自己的需要选择安装内容，但基本平台（即MATLAB选项）必须安装。第一次安装没有选择的内容在补安装时只需按照安装的过程进行，只是在选择组件时只勾选要补装的组件或工具箱即可。 1.4 MATLAB操作桌面有几个窗口？如何使某个窗口脱离桌面成为独立窗口？又如何将脱离出去的窗口重新放置到桌面上？ 答：在MATLAB操作桌面上有五个窗口。在每个窗口的右上角有两个小按钮，一个是关闭窗口的Close 按钮，一个是可以使窗口成为独立窗口的Undock按钮，点击Undock按钮就可以使该窗口脱离桌面成为独立窗口。在独立窗口的Desktop菜单中选择Dock...项就可以将独立的窗口重新放置到桌面上。
1.5 如何启动M文件编辑/调试器？ 答：在操作桌面上选择“建立新文件”或“打开文件”操作时，M文件编辑/调试器将被启动。在命令窗口中键入edit命令时也可以启动M文件编辑/调试器。 1.6 存储在工作空间中的数组能编辑吗？如何操作？ 答：存储在工作空间的数组可以通过数组编辑器进行编辑：在工作空间浏览器中双击要编辑的数组名打开数组编辑器，再选中要修改的数据单元，输入修改内容即可。 1.7 命令历史窗口除了可以观察前面键入的命令外，还有什么用途？ 答：命令历史窗口除了用于查询以前键入的命令外，还可以直接执行命令历史窗口中选定的内容、将选定的内容拷贝到剪贴板中、将选定内容直接拷贝到M文件中。 1.8 如何设置当前目录和搜索路径，在当前目录上的文件和在搜索路径上的文件有什么区别？ 答：当前目录可以在当前目录浏览器窗口左上方的输入栏中设置，搜索路径可以通过选择操作桌面的file 菜单中的Set Path菜单项来完成。在没有特别说明的情况下，只有当前目录和搜索路径上的函数和文件能够被MATLAB运行和调用，如果在当前目录上有与搜索路径上相同文件名的文件时则优先执行当前目录上的文件，如果没有特别说明，数据文件将存储在当前目录上。 1.9 在MATLAB中有几种获得帮助的途径？ 答：（1）帮助浏览器：选择view菜单中的Help菜单项或选择Help菜单中的MATLAB Help菜单项可以打开帮助浏览器。 （2）help命令：在命令窗口键入“help”命令可以列出帮助主题，键入“help 函数名”可以得到
附录1 使用牛顿拉夫逊法进行潮流计算的Matlab程序代码 % 牛拉法计算潮流程序 %----------------------------------------------------------------------- % B1矩阵:1、支路首端号;2、末端号;3、支路阻抗;4、支路对地电纳 % 5、支路的变比;6、支路首端处于K侧为1,1侧为0 % B2矩阵:1、该节点发电机功率;2、该节点负荷功率;3、节点电压初始值 % 4、PV节点电压V的给定值;5、节点所接的无功补偿设备的容量 % 6、节点分类标号:1为平衡节点(应为1号节点);2为PQ节点;3为PV 节点; %------------------------------------------------------------------------ clear all; format long; n=input('请输入节点数:nodes='); nl=input('请输入支路数:lines='); isb=input('请输入平衡母线节点号:balance='); pr=input('请输入误差精度:precision='); B1=input('请输入由各支路参数形成的矩阵:B1='); B2=input('请输入各节点参数形成的矩阵:B2='); Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl); %------------------------------------------------------------------ for i=1:nl %支路数 if B1(i,6)==0 %左节点处于1侧 p=B1(i,1);q=B1(i,2); else %左节点处于K侧 p=B1(i,2);q=B1(i,1); end Y(p,q)=Y(p,q)-1、/(B1(i,3)*B1(i,5)); %非对角元 Y(q,p)=Y(p,q); %非对角元 Y(q,q)=Y(q,q)+1、/(B1(i,3)*B1(i,5)^2)+B1(i,4); %对角元K侧 Y(p,p)=Y(p,p)+1、/B1(i,3)+B1(i,4); %对角元1侧 end %求导纳矩阵 disp('导纳矩阵Y='); disp(Y) %------------------------------------------------------------------- G=real(Y);B=imag(Y); %分解出导纳阵的实部与虚部 for i=1:n %给定各节点初始电压的实部与虚部 e(i)=real(B2(i,3)); f(i)=imag(B2(i,3));
matlab数值计算功能 1,基础运算 （1）多项式的创建与表达 将多项式(x-6)(x-3)(x-8)表示为系数形式 a=[6 3 8] % 写成根矢量 pa=poly(a)% 求出系数矢量 ppa=poly2sym(pa,'x') % 表示成符号形式 ezplot(ppa,[-50,50]) 求3介方阵A的特征多项式 a=[6 2 4;7 5 6;1 3 6 ]; pa=poly(a)% 写出系数矢量 ppa=poly2sym(pa) %表示成符号形式 ezplot(ppa,[-50,50]) % 绘图 求x^3-6x^2-72x-27的根。 a=[1,-6,-72,-85]; % 写出多项式系数矢量 r=roots(a) % 求多项式的根 （2）多项式的乘除运算 c=conv(a,b) %乘法 [q,r]=deconv(c,a)% 除法 求a(s)=s^2+2s+3乘以b(s)=4s^2+5s+6的乘积 a=[1 2 3] b=[4 5 6] % 写出系数矢量 c=conv(a,b) c=poly2sym(c,'s') % 写成符号形式的多项式 展开(s^2+2s+2)(s+4)(s+1)并验证结果被(s+4),(s+3)除后的结果。c=conv([1,2,2],conv([1,4],[1,1])) cs=poly2sym(c,'s') c=[1 7 16 18 8] [q1,r1]=deconv(c,[1,4]) [q2,r2]=deconv(c,[1,3]) cc=conv(q2,[1,3]) test=((c-r2)==cc)
其他常用的多项式运算命令 pa=polyval(p,s) % 按数组规则计算给定s时多项式的值 pm=polyvalm(p,s)% 按矩阵规则计算给定s时多项式的值 [r,p,k]=residue(b,a) % 部分分式展开，b,a分别是分子，分母多项式系数矢量。r,p,k分别是留数，极点和值项矢量。 p=poly(x,y,n) % 用n介多项式拟合给定的数据 polyder(p) %多项式微分 x=[1 2 3 4 5]; y=[5.5 43.1 128 290.7 498.4]; p=polyfit(x,y,3) x2=1:0.1:5; y2=polyval(p,x2); plot(x,y,'o',x2,y2) 2，线性代数 1，求解方程的根 t=[0 0.3 0.8 1.1 1.6 2.3]'; y=[0.82 0.72 0.63 0.60 0.55 0.50]'; e=[ones(size(t)) exp(-t)] c=e\y t1=[0:0.1:2.5]'; y1=[ones(size(t1)),exp(-t1)]*c; plot(t1,y1,'b',t,y,'ro') 2，逆矩阵及行列式 inv(a) det(a) pinv(a) 3，矩阵分解（略） 4，数据分析 max(x)求x各列的最大元素 mean(x)求x各列的平均值 median(x)找出x各列的中位元素 min(x)求出x各列的最小元素 S=cumsum(x)求x各列元素累计和 sort(x)使x的各列元素按递增排序。 std(x)求x各列的标准差。 sum(x)求x各列元素之和 prod(x)求x各列元素之积
2. 用MATLAB 语句输入矩阵A 和B 3．假设已知矩阵A ，试给出相应的MATLAB 命令，将其全部偶数行提取出来， 赋给B 矩阵，用magic(8)A =命令生成A 矩阵，用上述命令检验一下结果是不是正确。 4．用数值方法可以求出∑=++++++==63 63622284212i i S Λ，试不采用循环的 形式求出和式的数值解。由于数值方法是采用double 形式进行计算的，难以保证有效位数字，所以结果不一定精确。试采用运算的方法求该和式的精确值。 5．选择合适的步距绘制出下面的图形。 （1）)/1sin(t ，其中)1,1(-∈t ； （2）)tan(sin )sin(tan t t -，其中),(ππ-∈t 6. 试绘制出二元函数2 2 2 2 )1(1)1(1),(y x y x y x f z +++ +-= =的三维图和三 视图 7. 试求出如下极限。 （1）x x x x 1)93(lim +∞ →； （2）1 1lim 0-+→→xy xy y x ； （3）2 2)()cos(1lim 2 2 220 0y x y x e y x y x +→→++- 8. 已知参数方程? ??-==t t t y t x sin cos cos ln ，试求出x y d d 和3 /2 2d d π=t x y 9. 假设?-=xy t t e y x f 0 d ),(2 ，试求2222 22y f y x f x f y x ??+???-?? 10. 试求出下面的极限。 （1）??????-++-+-+-∞→1)2(1 161141121lim 2222n n Λ； （2）)131211( lim 2 222π πππn n n n n n n ++++++++∞ →Λ 11. 试求出以下的曲线积分。 （1）?+l s y x d )(22，l 为曲线)sin (cos t t t a x +=，)cos (sin t t t a y -=， )20(π≤≤t 。
基于MATLAB的电力系统潮流计算 %简单潮流计算的小程序，相关的原始数据数据数据输入格式如下： %B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点%编号，将变压器的串联阻抗置于低压侧处理。 %第三列为支路的串列阻抗参数。 %第四列为支路的对地导纳参数。 %第五烈为含变压器支路的变压器的变比 %第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，%“0”为不含有变压器。 %B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点%负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中 %“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。 %X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地%参数。 n=input('请输入节点数:n='); n1=input('请输入支路数:n1='); isb=input('请输入平衡节点号:isb='); pr=input('请输入误差精度:pr='); B1=input('请输入支路参数:B1='); B2=input('请输入节点参数:B2='); X=input('节点号和对地参数:X='); Y=zeros(n); Times=1; %置迭代次数为初始值 %创建节点导纳矩阵 for i=1:n1 if B1(i,6)==0 %不含变压器的支路 p=B1(i,1); q=B1(i,2); Y(p,q)=Y(p,q)-1/B1(i,3); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4); Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4); else %含有变压器的支路 p=B1(i,1); q=B1(i,2); Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5)); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1/B1(i,3);
Matlab 基础练习题 常量、变量、表达式 1、 MATLAB 中，下面哪些变量名是合法的？（ ） （A ）_num （B ）num_ （C ）num- （D ）-num 2、 在MA TLAB 中，要给出一个复数z 的模，应该使用（ ）函数。 （A ）mod(z) （B ）abs(z) （C ）double(z) （D ）angle(z) 3、 下面属于MATLAB 的预定义特殊变量的是？（ ） （A ）eps （B ）none （C ）zero （D ）exp 4、 判断：在MA TLAB 的内存工作区中，存放一个英文字符 'a' 需要占用1个字节，存放 一个中文字符‘啊’需要占用2个字节。（ 错，都是2个字节 ） 5、 判断：MA TLAB 中，i 和j （ 对 ） 6、 判断：MA TLAB 中，pi 代表圆周率，它等于3.14。（ 错，后面还有很多位小数 ） 7、 在MA TLAB 中，若想计算的5 1)3.0sin(21+= πy 值，那么应该在MA TLAB 的指令窗中 输入的MA TLAB 指令是__y1=2*sin(0.3*pi)/(1+sqrt(5))_。 8、 在MA TLAB 中，a = 1，b = i ，则a 占_8__个字节，b 占_16_个字节，c 占________字 节。 9、 在MA TLAB 中，inf 的含义是__无穷大__，nan 的含义是__非数（结果不定）___。 数组 1、 在MA TLAB 中，X 是一个一维数值数组，现在要把数组X 中的所有元素按原来次序 的逆序排列输出，应该使用下面的（ ）指令。 （A ）X[end:1] （B ）X[end:-1:1] （C ）X (end:-1:1) （D ）X(end:1) 2、 在MA TLAB 中，A 是一个字二维数组，要获取A 的行数和列数，应该使用的MATLAB 的命令是（ ）。 （A ）class(A) （B ）sizeof(A) （C ）size(A) （D ）isa(A) 3、 在MATLAB 中，用指令x=1:9生成数组x 。现在要把x 数组的第二和第七个元素都 赋值为0，应该在指令窗中输入（ ） （A ）x([2 7])=(0 0) （B ）x([2,7])=[0,0] （C ）x[(2,7)]=[0 0] （D ）x[(2 7)]=(0 0) 4、 在MA TLAB 中，依次执行以下指令：clear;A=ones(3,4); A(:)=[-6:5];这时， 若在指令窗中输入指令b=A(:,2)'，那么，MATLAB 输出的结果应该是（ ） （A ）b = -3 -2 -1 （B ）b = -2 -1 0 1 （C ）b = -5 -1 3 （D ）b = -5 -2 1 4 5、 在MA TLAB 中，A = 1:9，现在执行如下指令L1 = ~(A>5)，则MATLAB 的执行结果应 该是L1 =___ 1 1 1 1 1 0 0 0 0___。
%*************************电力系统直角坐标系下的牛顿拉夫逊法潮流计算********** clear clc load E:\data\IEEE014_Node.txt Node=IEEE014_Node; weishu=size(Node); nnum=weishu(1,1); %节点总数 load E:\data\IEEE014_Branch.txt branch=IEEE014_Branch; bwei=size(branch); bnum=bwei(1,1); %支路总数 Y=(zeros(nnum)); Sj=100; %********************************节点导纳矩阵******************************* for m=1:bnum; s=branch(m,1); %首节点 e=branch(m,2); %末节点 R=branch(m,3); %支路电阻 X=branch(m,4); %支路电抗 B=branch(m,5); %支路对地电纳 k=branch(m,6); if k==0 %无变压器支路情形 Y(s,e)=-1/(R+j*X); %互导纳 Y(e,s)=Y(s,e); end if k~=0 %有变压器支路情形 Y(s,e)=-(1/((R+j*X)*k)); Y(e,s)=Y(s,e); Y(s,s)=-(1-k)/((R+j*X)*k^2); Y(e,e)=-(k-1)/((R+j*X)*k); %对地导纳 end Y(s,s)=Y(s,s)-j*B/2; Y(e,e)=Y(e,e)-j*B/2; %自导纳的计算情形 end for t=1:nnum; Y(t,t)=-sum(Y(t,:))+Node(t,12)+j*Node(t,13); %求支路自导纳 end G=real(Y); %电导 B=imag(Y); %电纳 %******************节点分类************************************* * pq=0; pv=0; blancenode=0; pqnode=zeros(1,nnum); pvnode=zeros(1,nnum); for m=1:nnum; if Node(m,2)==3 blancenode=m; %平衡节点编号 else if Node(m,2)==0 pq=pq+1; pqnode(1,pq)=m; %PQ 节点编号 else if Node(m,2)==2 pv=pv+1; pvnode(1,pv)=m; %PV 节点编号 end end end end %*****************************设置电压初值********************************** Uoriginal=zeros(1,nnum); %对各节点电压矩阵初始化 for n=1:nnum Uoriginal(1,n)=Node(n,9); %对各点电压赋初值 if Node(n,9)==0;
东北电力大学课程设计改革试用任务书: 电力系统潮流计算课程设计任务书 设计名称：电力系统潮流计算课程设计 设计性质：理论计算，计算机仿真与验证 计划学时：两周 一、设计目的 1.培养学生独立分析问题、解决问题的能力； 2.培养学生的工程意识，灵活运用所学知识分析工程问题的能力 3.编制程序或利用电力系统分析计算软件进行电力系统潮流分析。 二、原始资料 1、系统图：IEEE14节点。 2、原始资料：见IEEE14节点标准数据库 三、课程设计基本内容： 1.采用PSAT仿真工具中的潮流计算软件计算系统潮流； 1)熟悉PSAT仿真工具的功能； 2)掌握IEEE标准数据格式内容； 3)将IEEE标准数据转化为PSAT计算数据； 2.分别采用NR法和PQ分解法计算潮流，观察NR法计算潮流中雅可比矩阵的变化情况， 分析两种方法计算潮流的优缺点； 3.分析系统潮流情况，包括电压幅值、相角，线路过载情况以及全网有功损耗情况。
4.选择以下内容之一进行分析： 1)找出系统中有功损耗最大的一条线路，给出减小该线路损耗的措施，比较各种措施 的特点，并仿真验证； 2)找出系统中电压最低的节点，给出调压措施，比较各种措施的特点，并仿真验证； 3)找出系统中流过有功功率最大的一条线路，给出减小该线路有功功率的措施，比较 各种措施的特点，并仿真验证； 5.任选以下内容之一作为深入研究：（不做要求） 1)找出系统中有功功率损耗最大的一条线路，改变发电机有功出力，分析对该线路有 功功率损耗灵敏度最大的发电机有功功率，并进行有效调整，减小该线路的损耗； 2)找出系统中有功功率损耗最大的一条线路，进行无功功率补偿，分析对该线路有功 功率损耗灵敏度最大的负荷无功功率，并进行有效调整，减小该线路的损耗； 3)找出系统中电压最低的节点，分析对该节点电压幅值灵敏度最大的发电机端电压， 并有效调整发电机端电压，提高该节点电压水平； 四、课程设计成品基本要求： 1.绘制系统潮流图，潮流图应包括： 1)系统网络参数 2)节点电压幅值及相角 3)线路和变压器的首末端有功功率和无功功率 2.撰写设计报告，报告内容应包括以下几点： 1)本次设计的目的和设计的任务； 2)电力系统潮流计算的计算机方法原理，分析NR法和PQ分解法计算潮流的特点； 3)对潮流计算结果进行分析，评价该潮流断面的运行方式安全性和经济性； 4)找出系统中运行的薄弱环节，如电压较低点或负载较大线路，给出调整措施； 5)分析各种调整措施的特点并比较它们之间的差异； 6)结论部分以及设计心得； 五、考核形式 1.纪律考核：学生组织出勤情况和工作态度等； 2.书面考核：设计成品的完成质量、撰写水平等； 3.答辩考核：参照设计成品，对计算机方法进行电力系统潮流计算的相关问题等进行答辩； 4.采用五级评分制：优、良、中、及格、不及格五个等级。
第6章 MATLAB 数值计算 例6.1 求矩阵A 的每行及每列的最大和最小元素，并求整个矩阵的最大和最小元素。 1356 78256323578255631 01-???? -? ?=???? -??A A=[13,-56,78;25,63,-235;78,25,563;1,0,-1]; max(A,[],2) %求每行最大元素 min(A,[],2) %求每行最小元素 max(A) %求每列最大元素 min(A) %求每列最小元素 max(max(A)) %求整个矩阵的最大元素。也可使用命令：max(A(:)) min(min(A)) %求整个矩阵的最小元素。也可使用命令：min(A(:)) 例6.2 求矩阵A 的每行元素的乘积和全部元素的乘积。 A=[1,2,3,4;5,6,7,8;9,10,11,12]; S=prod(A,2) prod(S) %求A 的全部元素的乘积。也可以使用命令prod(A(:)) 例6.3 求向量X =(1！，2！，3！，…，10！)。 X=cumprod(1:10) 例6.4 对二维矩阵x ，从不同维方向求出其标准方差。 x=[4,5,6;1,4,8] %产生一个二维矩阵x y1=std(x,0,1) y2=std(x,1,1) y3=std(x,0,2) y4=std(x,1,2) 例6.5 生成满足正态分布的10000×5随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。 X=randn(10000,5); M=mean(X) D=std(X) R=corrcoef(X)
例6.6 对下列矩阵做各种排序。 185412613713-?? ??=?? ??-?? A A=[1,-8,5;4,12,6;13,7,-13]; sort(A) %对A 的每列按升序排序 -sort(-A,2) %对A 的每行按降序排序 [X,I]=sort(A) %对A 按列排序，并将每个元素所在行号送矩阵I 例6.7 给出概率积分 2 (d x x f x x -? e 的数据表如表6.1所示，用不同的插值方法计算f (0.472)。 x=0.46:0.01:0.49; %给出x ，f(x) f=[0.4846555,0.4937542,0.5027498,0.5116683]; format long interp1(x,f,0.472) %用默认方法，即线性插值方法计算f(x) interp1(x,f,0.472,'nearest') %用最近点插值方法计算f(x) interp1(x,f,0.472,'spline') %用3次样条插值方法计算f(x) interp1(x,f,0.472,'cubic') %用3次多项式插值方法计算f(x) format short 例6.8 某检测参数f 随时间t 的采样结果如表6.2，用数据插值法计算t =2，7，12，17，22，17，32，37，42，47，52，57时的f 值。 T=0:5:65; X=2:5:57;
习题 1 1. 执行下列指令，观察其运算结果, 理解其意义： (1) [1 2;3 4]+10-2i (2) [1 2; 3 4].*[0.1 0.2; 0.3 0.4] (3) [1 2; 3 4].\[20 10;9 2] (4) [1 2; 3 4].^2 (5) exp([1 2; 3 4]) (6)log([1 10 100]) (7)prod([1 2;3 4]) (8)[a,b]=min([10 20;30 40]) (9)abs([1 2;3 4]-pi) (10) [1 2;3 4]>=[4,3;2 1] (11)find([10 20;30 40]>=[40,30;20 10]) (12) [a,b]=find([10 20;30 40]>=[40,30;20 10]) (提示：a 为行号，b 为列号) (13) all([1 2;3 4]>1) (14) any([1 2;3 4]>1) (15) linspace(3,4,5) (16) A=[1 2;3 4];A(:,2) 2. 执行下列指令，观察其运算结果、变量类型和字节数，理解其意义： (1) clear; a=1,b=num2str(a),c=a>0, a= =b, a= =c, b= =c (2) clear; fun='abs(x)',x=-2,eval(fun),double(fun) 3. 本金K 以每年n 次，每次p %的增值率(n 与p 的乘积为每年增值额的百分比)增加，当增加到rK 时所花费的时间为 ) 01.01ln(ln p n r T += (单位：年) 用MA TLAB 表达式写出该公式并用下列数据计算：r =2, p =0.5, n =12. 4．已知函数f (x )=x 4 -2x 在(-2, 2)内有两个根。取步长h =0.05, 通过计算函数值求得函数的最小值点和两个根的近似解。（提示：求近似根等价于求函数绝对值的最小值点） ? 5. (1) 用z=magic(10)得到10阶魔方矩阵； (2) 求z 的各列元素之和； (3) 求z 的对角线元素之和(提示：先用diag(z)提取z 的对角线)； (4) 将z 的第二列除以3;
毕业设计（论文）MATLAB下的潮流计算实现-稀疏技术
毕业设计（论文）原创性声明和使用授权说明 原创性声明 本人郑重承诺：所呈交的毕业设计（论文），是我个人在指导教师的指导下进行的研究工作及取得的成果。尽我所知，除文中特别加以标注和致谢的地方外，不包含其他人或组织已经发表或公布过的研究成果，也不包含我为获得及其它教育机构的学位或学历而使用过的材料。对本研究提供过帮助和做出过贡献的个人或集体，均已在文中作了明确的说明并表示了谢意。 作者签名：日期： 指导教师签名：日期： 使用授权说明 本人完全了解大学关于收集、保存、使用毕业设计（论文）的规定，即：按照学校要求提交毕业设计（论文）的印刷本和电子版本；学校有权保存毕业设计（论文）的印刷本和电子版，并提供目录检索与阅览服务；学校可以采用影印、缩印、数字化或其它复制手段保存论文；在不以赢利为目的前提下，学校可以公布论文的部分或全部内容。 作者签名：日期：
摘要 电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各母线的电压，各元件中流过的功率，系统的功率损耗等等。在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性、可靠性和经济性。因此潮流计算是研究电力系统的一种很重要和很基础的计算。由于电力系统结构及参数的一些特点，并且随着电力系统不断扩大，潮流问题的方程式阶数越来越高，对这样的方程式并不是任何数学方法都能保证给出正确答案的。这种情况成为促使电力系统计算人员不断寻求新的更可靠方法的重要因素。 本文旨在于研究潮流计算的牛顿—拉夫逊法的基本原理，在Matlab环境中实现牛顿—拉夫逊法潮流计算的数学模型，程序流程以及编制相应程序，并在程序中融合了节点优化编号和稀疏技术，以提高计算效率。最后用IEEE-3O节点标准测试系统验证所编程序。 关键词：潮流计算Newtom-Raphson法节点优化稀疏技术Matlab ABSTRACT Power flow calculation is fundanmental of analysis. Network reconfiguration,fault management,state estimator etc also need the data of electrial system power flow.There is important significance to develop power flow calculation in allusion to traits of distribution network. This paper introduces the principle of Newtom-Raphson algorithm, which is developed for calculation of power flow calculation ,where zero sequence network is open.With this algorithm,the three-phase load is resolved into positive/negative sequence power and coupling power,thus,decoupling three phase power flow into sequencet component power flow.The power flow can be obtained by just finding the positive sequence power flow and then finding the negative sequent component from the coupling https://www.doczj.com/doc/ad8246407.html,pared with the existing methods,the jacobian matrix with the proposed algorithm is of much lower order,thus substantially reducing the computation burden.The proposed algorithm,together with a reference algorithm,has been simulated on an actual IEEE-30 system using statistic load date.And then it will
列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1);
end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x')
一、填空题 1. MATLAB命令窗口中可用clear 命令清除工作区中的变量；用clc 命令清除命令窗口中的内容。 2. MATLAB中的运算符包括数学运算符、矩阵运算符、数组运算符。 3 MATLAB的含义为矩阵运算。 4. MATLAB中的函数类型有嵌套函数、主函数、子函数、私有函数和匿名函数。 5. MATLAB中用系统预定义的特殊变量nargin 表示函数输入参数个数。 6. size命令的功能确定数组大小。 7. MATLAB中用linspace 符号定义单元数组。 8. MATLAB中用hold on 命令打开图形保持功能。 9 MATLAB中如果一个变量名为123，该变量名是否合法否（填是或否）。 10. MATLAB中用inv 函数完成矩阵的求逆运算，用det 函数求矩阵的行列式。 11.已知C={2+5i,’Hi’;1:4,[1 2]} 输入s=C(3) 即s=__'_Hi__'___. 12.先输入P=struct(‘name’,’Jack‘,’num’,43) ，再输入s= 即s=___Jack___. 13.当x=[ 1 2 3; 4 5 6] , y=[ 5 8; 5 8 ] 时，__不能__(能﹨不能)使用plot 命令绘图。 14．sym 或syms _能_______(能﹨不能)创建符号矩阵。 15.已知x=3:2:24,则x有__11_个元素。 结构的开始是“if”命令，结束是_end____命令。 17.已知x为一个向量，计算ln（x）的MATLAB命令是计算_log(x)___。 18.命令：>>subplot(3,3,x)中有__9_个子图。 19.已知s=‘显示“hey”’，则s的元素个数是___9__。 y=f() global W W=2 y=f1(5) function y=f() global W y=w+W 创建以上函数文件，在命令窗口中运行“y=f”命令则显示__出错___。 21. MATLAB中用find 函数查找数组中的非零元素，用sort 函数对数组进行排序。 22. MATLAB中的M文件包括M脚本文件和M函数文件。 23. MATLAB中用系统预定义的特殊变量eps 表示浮点数的精度。 24. MATLAB中8位有符号整数用int8 定义。 25. MATLAB中小括号、转置运算符和冒号运算符，具有最高优先级的为转置运算符。 26. MATLAB中用系统预定义的特殊变量nargout 表示函数输出参数个数。 27. MATLAB中用% 符号表示其后为程序注释。 28. 已知矩阵A,则矩阵A的共轭转置可表示为A’。 29. MATLAB中用grid on 命令显示坐标轴下的网格线。 30. MATLAB中如果一个函数名为test，该函数名是否合法是（填是或否）。