第12章例题源程序

1、例12.1 三相不平衡交流电路分析。计算图12.1所示三相不平衡交流电路各支路电流()并绘制相量图。其中ra= rb= rc=5, rab=6, rbc=10, rca=15,E=220V 。用MATLAB先求，进而求各支路电流:a=16,-5,-5;-5,20,-5;-5,-5,25; %系数矩阵b=220;-110-110*sqrt(3)*i;-110+110*sqrt(3)*i;I=inv(a)*b; %解方程Ia=I(1)-I(3) %求各支路电流Ib=I(2)-I(1)Ic=I(3)-I(2)h=compass(Ia,Ib,Ic); %绘制相量图set(h,LineWidth,2);

2、例12.2 调谐振荡电路分析。分析图12.3所示的调谐振荡电路（iL=f(v)=，0，0），要求绘制振荡波形和相轨迹。 基于以上状态方程建立函数文件vdpol.m：function ydot=vdpol(t,y)ydot(1)=0.1*(1-y(2)2)*y(1)-y(2); %的值可以任意变化，此处取0.1ydot(2)=y(1);ydot=ydot; 求解微分方程，并绘制振荡波形(t,y)和相轨迹(y,dy/dt)：t0=0;tf=60; %确定积分区间y0=0;0.25; %确定初始条件t,y=ode45(vdpol,t0,tf,y0); %求解微分方程subplot(1,2,1);pl

3、ot(t,y(:,2); %绘制振荡波形subplot(1,2,2);plot(y(:,2),y(:,1); %绘制相轨迹例12.3 一典型线性反馈控制系统结构如图12.6所示。 解法2：利用MATLAB控制系统工具箱中已经定义的一些LTI仿真函数，编写如下程序：G=tf(4,1,2,3,4); Gc=tf(1,-3,1,3); H=tf(1,0.01,1); G_o = Gc*G ; %构造开环系统的传递函数。G_c = feedback(G_o, H); %构造闭环系统的传递函数。step(G_o); %求开环系统的阶跃响应并绘制相应的曲线。figure;step(G_c); %求闭环系统

4、的阶跃响应并绘制相应的曲线。例12.4 自行车轮饰物的运动轨迹。为了使平淡的自行车增添一份美感，同时，也为了增加自行车的安全系数，一些骑车的人及自行车厂家在自行车的辐条上安装一款亮丽夺目的饰物，当有这种饰物的自行车在马路上驶过时，这饰物就象游龙一样，对街边的行人闪过一道波浪形的轨迹。这一波一闪的游龙，其轨迹是什么曲线？试画出它的图形。当自行车在一个抛物线型的拱桥上通过时，或是在一拱一拱的正弦曲线上通过时，这轨迹是什么曲线?试画出其图形。程序如下：clear;%第(1)类情况的实现x0=0:0.01:2;R=0.1;r=0.075;x1=x0-r*sin(x0/R); %计算f(x)=0时p点运

5、动轨迹y1=R-r*cos(x0/R);subplot(3,1,1);plot(x1,y1,x0,0); %绘制运动轨迹曲线和f(x)曲线xlabel(x1);ylabel(y1);grid on%第(2)类情况的实现x0=-1:0.01:1;R=0.1;r=0.1;y0=0.2-0.2*x0.2; %计算路面曲线fai=atan(-0.4*x0); %求int=inline(sqrt(1+(-0.4*x).2);%定义的积分函数for k=1:length(x0) theta1(k)=quad(int,0,x0(k)/R; %调用quad函数求endx2=x0+R*0.4*x0./sqrt(

6、1+(-0.4*x0).2)-r*sin(theta1-fai);%运动轨迹方程y2=y0+R./sqrt(1+(-0.4*x0).2)-r*cos(theta1-fai);subplot(3,1,2);plot(x2,y2,x0,y0) %绘制运动轨迹曲线和f(x)曲线xlabel(x2);ylabel(y2);grid on%第(3)类情况的实现x0=0:0.01:10;R=0.1;r=0.075;y0=0.3*sin(x0); %计算路面曲线fai=atan(0.3*cos(x0); %求int=inline(sqrt(1+(0.3*cos(x).2);%定义的积分函数for k=1:l

7、ength(x0) theta2(k)=quad(int,0,x0(k)/R; %调用quad函数求end%p点运动轨迹方程x3=x0-0.3*R*cos(x0)./sqrt(1+(0.3*cos(x0).2)-r*sin(theta2-fai); y3=0.3*sin(x0)+R./sqrt(1+(0.3*cos(x0).2)-r*cos(theta2-fai);subplot(3,1,3);plot(x3,y3,x0,y0)xlabel(x3),ylabel(y3);grid on例12.5 广告费用与效应。某装饰材料公司欲以每桶2元的价钱购进一批彩漆。一般来说，随着彩漆售价的提高，预期销

8、售量将减少，并对此进行了估算，见表12.1。现在的问题是装饰材料公司采取怎样的营销战略使得预期的利润最大?x1=2:0.5:6;y1=41000,38000,34000,32000,29000,28000,25000,22000,20000;x2=0:10000:70000;y2=1.0,1.4,1.7,1.85,1.95,2.00,1.95,1.8;subplot(2,1,1);plot(x1,y1,o);title(售价与预期销售量);subplot(2,1,2);plot(x2,y2,o);title(广告费与销售增长因子); 首先利用MATLAB编程计算参数a,b,c,d,e。form

9、at long;x1=2:0.5:6;y1=41000,38000,34000,32000,29000,28000,25000,22000,20000;x2=0:10000:70000;y2=1.0,1.4,1.7,1.85,1.95,2.00,1.95,1.8;a1=polyfit(x1,y1,1);a2=polyfit(x2,y2,2);disp(a1),disp(a2)建立函数文件p.m：function f=p(x)f=x(2)-(1.01875+4.09226e-5*x(2)-4.2560e-10*x(2)2)* .(50422.2-5133.33*x(1)*(x(1)-2);利用f

10、minsearch函数寻优：format short ex,f=fminsearch(p,0,0) %寻优从坐标原点开始例12.6 考虑空气阻力时抛射体质心的飞行轨迹问题。假设空气阻力的方向与速度向量相反，大小与速度的平方成正比。抛射体的受力情况如图12.13所示。计算质点飞行的轨迹和距离。 (1)函数文件zf.m：function rdot=zf(t,r)c=0.02;g=9.8;m=1;vm=sqrt(r(3)2+r(4)2);rdot=r(3);r(4);-c*vm*r(3)/m; -c*vm*r(4)/m-g;(2)主程序ex12_6.m：clear;y0=0;x0=0; %初始位置v

11、0=input(请输入初始输入速度(m/s)：);rho= input(请输入初始方向(度)：);tf=input(请输入飞行时间(s)：);vx0=v0*cos(rho*pi/180); %计算x,y方向的初始速度vy0=v0*sin(rho*pi/180);t,r=ode45(zf,0,tf,0;0;vx0;vy0); %解微分方程H=max(r(:,2) %求轨道的最高点T=t(find(r(:,2)=H) %到最高点的所需时间L=min(r(find(r(:,2)0),1) %计算射程plot(0,100,0,0);hold on %绘制X坐标线xlabel(x);ylabel(y);

12、plot(r(:,1),r(:,2); %绘制运动轨迹例12.7 静不定问题。假设有一结构如图12.15(a)所示，5根直杆连接于节点A，B，C，D，E和O。O点受一个垂直向下的力F作用，各杆具有相同的弹性模量E和横截面积A。假设连接点允许旋转，所以杆上没有力矩作用。设5根杆的长度分别为L1,L2,L3,L4,L5。这里，L1=L5且L2=L4。5根杆上受到的力分别为F1, F2, F3, F4,F5，如图12.15(b)所示。由结构的对称性可知，F1=F5，F2=F4。所以总共有3个未知力F1,F2,F3。设杆L1和L3所构成的夹角AOC=1，L2和L3所构成的夹角BOC=2。同样EOC=1

13、，DOC=2。求F1,F2,F3。 alpha1=atan(1/2);alpha2=atan(1/4);L3=1000;F=25000;E=;A=100;L1=L3/cos(alpha1);L2=L3/cos(alpha2);C=2*cos(alpha1),2*cos(alpha2),1,0;L1/(E*A),0,0,-cos(alpha1);.0,L2/(E*A),0,-cos(alpha2);0,0,L3/(E*A),-1;B=F;0;0;0;X=CB例12.8 简支梁荷载情况如图12.17所示。求其弯矩、转角和挠度。已知L=8m，q=500N/m，F0=1000N，M0=800Nm，E=

14、200109N/m，I=210-6m4。 由此编写MATLAB程序如下：L=8;q=500;F0=1000;M0=800;E=200e9;I=2e-6;N1=(3*q*L2/8+F0*L/2-M0)/L; N2=(q*L2/8+F0*L/2+M0)/L; x=linspace(0,L,101); dx=L/100;M1=N1*x(1:51)-q*x(1:51).2/2; %分3段用数组列出M的表达式M2=N2*(L-x(52:76)-M0;M3=N2*(L-x(77:101);M=M1,M2,M3;A0=cumtrapz(M)*dx/(E*I); %由M积分求转角(未记积分常数)Y0=cumt

15、rapz(A0)*dx; %由转角积分求挠度(未记积分常数)C=0,1;L,1-Y0(1);-Y0(101); %由边界条件求积分常数A=A0+C(1);Y=Y0+C(1)*x+C(2); %求转角和挠度的完整值subplot(3,1,1); plot(x, M); grid on;subplot(3,1,2); plot(x, A); grid on;subplot(3,1,3); plot(x, Y); grid on; 例12.9 计算图12.19所示直梁的自振频率(基频)。已知E=21011Pa，=7860kg/m3，A=1.1410-3m2，I=1.0310-7m4。编写程序如下：E

16、=2e11;I=1.03e-7;A=1.14e-3;rho=7860;L=1;X=12/L/L/L,-6/L/L,-12/L/L/L,-6/L/L;Ke=E*I*X;-6/L/L,4/L,6/L/L,2/L;-X;-6/L/L,2/L,6/L/L,4/L;K=Ke;KK=Ke;for i=2:4 %生成整体刚度矩阵 K=K,zeros(2*i,1),zeros(2*i,1); K=K;zeros(1,2*i+2);zeros(1,2*i+2); KK=zeros(2*i,1),zeros(2*i,1),KK; KK=zeros(1,2*i+2);zeros(1,2*i+2);KK; K=K+K

17、K;endMe=156,-22*L,54,13*L;-22*L,4*L*L,-13*L,-3*L*L;54,.-13*L,156,22*L;13*L,-3*L*L,22*L,4*L*L;Me=rho*A*L/420*Me;M=Me;MM=Me;for i=2:4 %生成整体质量矩阵 M=M,zeros(2*i,1),zeros(2*i,1); M=M;zeros(1,2*i+2);zeros(1,2*i+2); MM=zeros(2*i,1),zeros(2*i,1),MM; MM=zeros(1,2*i+2);zeros(1,2*i+2);MM; M=M+MM;end%处理约束条件K(1,:

18、)=;K(:,1)=; K(8,:)=;K(:,8)=; K(8,:)=;K(:,8)=; M(1,:)=;M(:,1)=; M(8,:)=;M(:,8)=;M(8,:)=;M(:,8)=; AA=inv(M)*K;W=eig(AA); %求特征值W=sqrt(W);disp(W);例12.10 无约束优化问题。已知梯形截面管道(图12 .21所示)的参数是：底边长度是c，高度是h，斜边与底边的夹角是，横截面积是A=64516mm2。管道内液体的流速与管道截面的周长s的倒数成比例关系。试按照使液体流速最大确定该管道的参数。这是一个二维无约束非线性优化问题。(1)建立目标函数文件fminx .m。function f=fminx(x)a=64516;f=a/x(1)-x(1)/tan(x(2)*pi/180)+2*x(1)/sin(x(2)*pi/180);(2) 调用fminsearch函数求解。x0=25;45;x,fval=fminsearch(fminx,x0)例12.11 证券投资组合问题。设金融市场上有2种风险证券A和B，它们的期望收益率分别为rA=12%，rB=18%，方差分别为，。同时市场上还有一种无风险证券，其收益率为rf=6%，设计一种投资组合方案，使得风险最