您好,欢迎访问三七文档
例12.1三相不平衡交流电路分析。计算图12.1所示三相不平衡交流电路各支路电流(cbaIII...,,)并绘制相量图。其中ra=rb=rc=5Ω,rab=6Ω,rbc=10Ω,rca=15Ω,E=220V。用MATLAB先求cabcabIII...,和,进而求各支路电流: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);例12.2调谐振荡电路分析。分析图12.3所示的调谐振荡电路(iL=f(v)=3vv,β0,γ0),要求绘制振荡波形和相轨迹。基于以上状态方程建立函数文件vdpol.m:functionydot=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);plot(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);%求闭环系统的阶跃响应并绘制相应的曲线。例12.4自行车轮饰物的运动轨迹。为了使平淡的自行车增添一份美感,同时,也为了增加自行车的安全系数,一些骑车的人及自行车厂家在自行车的辐条上安装一款亮丽夺目的饰物,当有这种饰物的自行车在马路上驶过时,这饰物就象游龙一样,对街边的行人闪过一道波浪形的轨迹。这一波一闪的游龙,其轨迹是什么曲线?试画出它的图形。当自行车在一个抛物线型的拱桥上通过时,或是在一拱一拱的正弦曲线上通过时,这轨迹是2什么曲线?试画出其图形。程序如下:clear;%第(1)类情况的实现x0=0:0.01:2;R=0.1;r=0.075;x1=x0-r*sin(x0/R);%计算f(x)=0时p点运动轨迹y1=R-r*cos(x0/R);subplot(3,1,1);plot(x1,y1,x0,0);%绘制运动轨迹曲线和f(x)曲线xlabel('x1');ylabel('y1');gridon%第(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)');%定义θ的积分函数fork=1:length(x0)theta1(k)=quad(int,0,x0(k))/R;%调用quad函数求θendx2=x0+R*0.4*x0./sqrt(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');gridon%第(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)');%定义θ的积分函数fork=1:length(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');gridon例12.5广告费用与效应。某装饰材料公司欲以每桶2元的价钱购进一批彩漆。一般来说,随着彩漆售价的提高,预期销售量将减少,并对此进行了估算,见表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('售价与预期销售量');3subplot(2,1,2);plot(x2,y2,'o');title('广告费与销售增长因子');首先利用MATLAB编程计算参数a,b,c,d,e。formatlong;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:functionf=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);利用fminsearch函数寻优:formatshorte[x,f]=fminsearch('p',[0,0])%寻优从坐标原点开始例12.6考虑空气阻力时抛射体质心的飞行轨迹问题。假设空气阻力的方向与速度向量相反,大小与速度的平方成正比。抛射体的受力情况如图12.13所示。计算质点飞行的轨迹和距离。(1)函数文件zf.m:functionrdot=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;%初始位置v0=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]);holdon%绘制X坐标线xlabel('x');ylabel('y');plot(r(:,1),r(:,2));%绘制运动轨迹例12.7静不定问题。假设有一结构如图12.15(a)所示,5根直杆连接于节点A,B,4C,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,∠DOC=α2。求F1,F2,F3。alpha1=atan(1/2);alpha2=atan(1/4);L3=1000;F=25000;E=205000;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=C\B例12.8简支梁荷载情况如图12.17所示。求其弯矩、转角和挠度。已知L=8m,q=500N/m,F0=1000N,M0=800N·m,E=200×109N/m,I=2×10-6m4。由此编写MATLAB程序如下:L=8;q=500;F0=1000;M0=800;E=200e9;I=2e-6;N1=(3*q*L^2/8+F0*L/2-M0)/L;N2=(q*L^2/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=cumtrapz(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);gridon;subplot(3,1,2);plot(x,A);gridon;subplot(3,1,3);plot(x,Y);gridon;例12.9计算图12.19所示直梁的自振频率(基频)。已知E=2×1011Pa,ρ=7860kg/m3,A=1.14×10-3m2,I=1.03×10-7m4。编写程序如下:E=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;fori=2:4%生成整体刚度矩阵K=[K,zeros(2*i,1),zeros(2*i,1)];K=[K;zeros(1,2*i+2);zeros(1,2*i+2)];5KK=[zeros(2*i,1),zeros(2*i,1),KK];KK=[zeros(1,2*i+2);zeros(1,2*i+2);KK];K=K+KK;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
本文标题:第12章例题源程序
链接地址:https://www.777doc.com/doc-2242627 .html