关 键 词：

卫国 MATLAB 程序设计 应用

资源描述：

刘卫国-MATLAB程序设计与应用（第2版）,卫国,MATLAB,程序设计,应用

内容简介：

例10.1 设系统的微分方程为：试建立系统仿真模型。例10.2 利用Simulink构建函数曲线y=5t2+16。例10.3 利用Simulink仿真求I=。例10.4 有初始状态为0的二阶微分方程x+0.2x+0.4x=0.2u (t), 其中u(t)是单位阶跃函数，试建立系统模型并仿真。例10.5 PID控制器是在自动控制中经常使用的模块，在工程应用中其标准的数学模型为其中采用了一阶环节来近似纯微分动作，为保证有良好的微分近似的效果，一般选N10。试建立PID控制器的模型并建立子系统。例10.6 利用使能子系统构成一个正弦半波整流器。例10.7 利用触发子系统将一锯齿波转换成方波。例10.8 采用S函数实现模块y=nx，即模块的功能是把一个输入信号n倍以后再输出。(1) 利用MATLAB语言编写S函数。程序如下：%*%S函数timesn.m，其输出是输入的n倍%*function sys,x0,str,ts=timesn(t,x,u,flag,n)switch flag,case 0 %初始化sys,x0,str,ts=mdlInitializeSizes;case 3 %计算输出量sys=mdlOutputs(t,x,u,n);case 1,2,4,9 %未使用的flag 值sys=;otherwise %出错处理error(Unhandle flag=,num2str(flag);end%*%mdlInitializeSizes：当flag为0 时进行整个系统的初始化%*function sys,x0,str,ts=mdlInitializeSizes(T)%调用函数simsizes以创建结构体sizessizes=simsizes;%用初始化信息填充结构体sizessizes.NumContStates=0; %无连续状态sizes.NumDiscStates=0; %无离散状态sizes.NumOutputs=1; %有一个输出量sizes.NumInputs=1; %有一个输入信号sizes.DirFeedthrough=1; %输出量中含有输入量sizes.NumSampleTimes=1; %单个采样周期%根据上面的设置设定系统初始化参数sys=simsizes(sizes);%给其他返回参数赋值。x0=; %设置初始状态为零状态str=; %将str变量设置为空字符串ts=-1,0; %假定继承输入信号的采样周期%初始化子程序结束%*%mdlOutputs：当flag值为3时，计算输出量%*function sys=mdlOutputs(t,x,u,n)sys=n*u;%输出量计算子程序结束。(2) 模块的封装与测试。例10.9 采用S函数来构造非线性分段函数。(1) 利用MATLAB语言编写S函数。程序如下：function sys,x0,str,ts=sfunction(t,x,u,flag)switch flag,case 0, sys,x0,str,ts=mdlInitializeSizes;case 3, sys=mdlOutputs(t,x,u);case 1,2,4,9 sys=;otherwise error(Unhandled flag=,num2str(flag);endfunctionsys,x0,str,ts=mdlInitializeSizessizes= simsizes;sizes.NumContStates= 0;sizes.NumDiscStates= 0;sizes.NumOutputs= 1;sizes.NumInputs= 1;sizes.DirFeedthrough= 1;sizes.NumSampleTimes= 1;sys=simsizes(sizes);x0=;str=;ts=0 0;function sys=mdlOutputs(t,x,u)if u=1&u=3&u=4&u=5&u6 sys=2-(u-5)2;else sys=1;end(2) 模块的测试。3例11.1 假定文件textdemo.txt中有以下格式的数据： names, x, y, z = textread (textdemo.txt, %s %d %d %d, 4, headerlines, 1);例11.2 计算当x=0.0,0.1,0.2,1.0时，f(x)=ex的值，并将结果写入文件demo1.txt。x=0:0,1:1; y=x;exp(x);fid=fopen(demo1.txt,w);fprintf(fid,%6.2f %12.8fn,y);fclose(fid);例11.3 将十进制数转换为十六进制数。A=6 10 14 44;fprintf(%9Xn, A+ (A0)*232)例11.4 读出例11.2生成的文件demo1.txt中的数据。 fid=fopen(demo1.txt,r); while 1LINE=fgetl(fid); if LINE0 break, end;disp (LINE);end fclose(fid);例11.5 将例11.1的文件读入到grades中。fid = fopen(textdemo.txt, r);grades = textscan(fid, %s %d %d %d, 3, headerlines, 1);grades:例11.6 假设文件alphabet.txt的内容是按顺序排列的26个大写英文字母，读取前5个字母的ASCII码和这5个字符。fid = fopen(alphabet.txt, r);c = fread(fid, 5);frewind(fid);d = fread(fid, 5, *char);fclose(fid);例11.7 建立一数据文件magic5.dat，用于存放5阶魔方阵。fid=fopen(magic5.dat,w);cnt=fwrite(fid,magic(5),int32);fclose(fid);例11.8 下列程序执行后，变量four、position和three的值是多少？a=1:5;fid=fopen(fdat.bin,w); %以写方式打开文件fdat.binfwrite(fid,a,int16); %将a的元素以双字节整型写入文件fdat.binstatus=fclose(fid); fid=fopen(fdat.bin,r); %以读数据方式打开文件fdat.binstatus=fseek(fid,6,bof); %将文件指针从开始位置向尾部移动6个字节four=fread(fid,1,int16); %读取第4个数据，并移动指针到下一个数据position=ftell(fid); %ftell的返回值为8 status=fseek(fid,-4,cof); %将文件指针从当前位置往前移动4个字节three=fread(fid,1,int16); %读取第3个数据status=fclose(fid);例11.9 创建MAT文件。#include #include /* For strcmp() */#include /* For EXIT_FAILURE, EXIT_SUCCESS */#include mat.h#define BUFSIZE 256int main() MATFile *pmat; /* 定义MAT文件指针*/ mxArray *pa1, *pa2, *pa3; double data9 = 1.0, 4.0, 7.0, 2.0, 5.0, 8.0, 3.0, 6.0, 9.0 ; const char *file = mattest.mat; char strBUFSIZE; int status; /* 打开一个MAT文件，如果不存在则创建一个MAT文件，如果打开失败，则返回 */ printf(Creating file %s.nn, file); pmat = matOpen(file, w); if (pmat = NULL) printf(Error creating file %sn, file); printf(Do you have write permission in this directory?)n);return(EXIT_FAILURE); /* 创建三个mxArray结构体对象，其中pa1、pa2分别为33、22的双精度实型矩阵，*/* pa3为字符串类型的阵列，如果创建失败则返回 */ pa1 = mxCreateDoubleMatrix(3,3,mxREAL); if (pa1 = NULL) printf(%s : Out of memory on line %dn, _FILE_, _LINE_); printf(Unable to create mxArray.n); return(EXIT_FAILURE); pa2 = mxCreateDoubleMatrix(3,3,mxREAL); if (pa2 = NULL) printf(%s : Out of memory on line %dn, _FILE_, _LINE_); printf(Unable to create mxArray.n); return(EXIT_FAILURE); memcpy(void *)(mxGetPr(pa2), (void *)data, sizeof(data); pa3 = mxCreateString(MATLAB: the language of technical computing); if (pa3 = NULL) printf(%s : Out of memory on line %dn, _FILE_, _LINE_); printf(Unable to create string mxArray.n); return(EXIT_FAILURE); /* 向MAT文件中写数据，失败则返回 */ status = matPutVariable(pmat, LocalDouble, pa1); if (status != 0) printf(%s : Error using matPutVariable on line %dn, &_FILE_, _LINE_); return(EXIT_FAILURE); status = matPutVariableAsGlobal(pmat, GlobalDouble, pa2); if (status != 0) printf(Error using matPutVariableAsGlobaln); return(EXIT_FAILURE); status = matPutVariable(pmat, LocalString, pa3); if (status != 0) printf(%s : Error using matPutVariable on line %dn, &_FILE_, _LINE_); return(EXIT_FAILURE); /* 清除矩阵 */ mxDestroyArray(pa1); mxDestroyArray(pa2); mxDestroyArray(pa3); /* 关闭MAT文件 */ if (matClose(pmat) != 0) printf(Error closing file %sn,file); return(EXIT_FAILURE); printf(Donen); return(EXIT_SUCCESS);3FORTRAN语言MAT文件应用(1)常用MAT文件操作函数常用FORTRAN语言MAT函数有：Integer*4 function matOpen(filename,mode)Integer*4 function matClose(mfp)Integer*4 function matDeleteVariable(mfp, name)Integer*4 function matGetDir(mfp,num)Integer*4 function matGetVariable(mfp, name)Integer*4 function matPutVariable(mfp, name, pm)Integer*4 function matPutVariableAsGlobal(mfp, name, pm)其中mfp为指向MAT文件的指针，name是读取或写入MAT文件的矩阵名称；num为MAT文件中矩阵的数目，pm是mxArray结构体的指针，mode为打开文件的方式。(2)应用示例例11.10 读取例11.9产生的MAT文件。 integer matOpen, matGetDir, matGetVariable integer mp, dir, adir(100), pa integer mxGetM, mxGetN, matClose integer ndir, i, stat character*32 names(100) mp = matOpen(mattest.mat, r) if (mp .eq. 0) then write(6,*) Cant open mattest.mat. stop end if! 读取目录 dir = matgetdir(mp, ndir) if (dir .eq. 0) then write(6,*) Cant read directory. stop endif call mxCopyPtrToPtrArray(dir, adir, ndir) do i=1,ndir call mxCopyPtrToCharacter(adir(i), names(i), 32) end do write(6,*) Directory of Mat-file: do i=1, ndir write(6,*) names(i) end do !读取矩阵 write(6,*) Getting array contents: do n=1,ndir pa = matGetVariable(mp, names(n) write(6,*) Retrieved , names(n) write(6,*) With size , mxGetM(pa), -by-, mxGetN(pa) call mxDestroyArray(pa) end do stat = matClose(mp) if (stat .ne. 0) then write(6,*) Error closing mattest.mat. end if end例11.11 创建一个矩阵，然后送到MATLAB计算引擎的工作区中，绘制出结果图。#include #include #include #include engine.h#define BUFSIZE 256int main()Engine *ep; /定义MATLAB引擎变量mxArray *T = NULL, *result = NULL;double time10 = 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 ;/*启动MATLAB计算引擎。如果在本地启动，那么函数所带的参数字符串为空；*/*如果在网络中启动，则需要提供服务器名，即engOpen(“服务器名”)*/if (!(ep = engOpen(0) fprintf(stderr, nCant start MATLAB enginen);return EXIT_FAILURE;/*向新启动的MATLAB工作区放置数据 */T = mxCreateDoubleMatrix(1, 10, mxREAL);memcpy(void *)mxGetPr(T), (void *)time, sizeof(time);engPutVariable(ep, T, T); /*执行MATLAB命令 */engEvalString(ep, D = .5.*(-9.8).*T.2;);engEvalString(ep, plot(T,D););engEvalString(ep, title(Position vs. Time for a falling object););engEvalString(ep, xlabel(Time (seconds););engEvalString(ep, ylabel(Position (meters););/*从MATLAB工作区获取计算结果 */result = engGetVariable(ep,D);printf(The size of result is %d*%dtn,mxGetM(result),mxGetN(result);printf(Hit return to continuenn);fgetc(stdin);printf(Done!n);/*释放内存空间，关闭计算引擎 */mxDestroyArray(T);engEvalString(ep, close;);engClose(ep);return EXIT_SUCCESS; 例11.12 在 FORTRAN源程序中调用 MATLAB 引擎绘制多峰函数矩阵。integer engEvalString, engClose, engOpen integer ep, status ep = engOpen(matlab) if ( ep=0) then write(6,*) Cant start MATLAB engine stop endifstatus = engEvalString(ep,mesh(peaks);) pauseif (engclose(ep)/=0) then write(6,*) Cant close MATLAB engine endif end例11.13 一个MATLAB自带的C语言MEX文件的编辑执行。#include mex.hvoid timestwo(double y, double x) y0 = 2.0*x0;void mexFunction( int nlhs, mxArray *plhs, int nrhs, const mxArray *prhs ) double *x,*y; int mrows,ncols; /* 检查正确的参数数目 */ if(nrhs!=1) mexErrMsgTxt(One input required.); else if(nlhs1) mexErrMsgTxt(Too many output arguments); /* 输入变量必须是非复数类型的标量 */ mrows = mxGetM(prhs0); ncols = mxGetN(prhs0); if( !mxIsDouble(prhs0) | mxIsComplex(prhs0) | !(mrows=1 & ncols=1) ) mexErrMsgTxt(Input must be a noncomplex scalar double.); /* 为返回参数创建矩阵 */ plhs0 = mxCreateDoubleMatrix(mrows,ncols, mxREAL); /* 分配输入输出参数的指针 */ x = mxGetPr(prhs0); y = mxGetPr(plhs0); /* 调用timestwo子例程 */ timestwo(y,x);例11.14 将例11.13用FORTRAN 90实现。subroutine mexFunction(nlhs, plhs, nrhs, prhs)integer plhs(*), prhs(*)integer mxGetPr, mxCreateDoubleMatrixinteger nlhs, nrhsinteger mxGetM, mxGetN, mxIsNumericinteger m, n , x, y! 检查输入、输出的参数个数if(nrhs .ne. 1.or. nlhs .ne. 1) thencall mexErrMsgTxt(One input required.)endif! 获得输入矩阵的大小m = mxGetM(prhs(1)n = mxGetN(prhs(1)! 检查输入数据的类型是否为数值if(mxIsNumeric(prhs(1) .eq. 0) thencall mexErrMsgTxt(Input must be a number.)endif! 创建输出矩阵plhs(1) = mxCreateDoubleMatrix(m,n,0)! 调用计算子例程x = mxGetPr(prhs(1)y = mxGetPr(plhs(1)call timestwo(%VAL(y), %VAL(x)! 将数据送往MATLAB工作空间returnendsubroutine timestwo(y, x)real*8 x, yy = 2.0 * xreturnend6例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);例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);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 自行车轮饰物的运动轨迹。为了使平淡的自行车增添一份美感，同时，也为了增加自行车的安全系数，一些骑车的人及自行车厂家在自行车的辐条上安装一款亮丽夺目的饰物，当有这种饰物的自行车在马路上驶过时，这饰物就象游龙一样，对街边的行人闪过一道波浪形的轨迹。这一波一闪的游龙，其轨迹是什么曲线？试画出它的图形。当自行车在一个抛物线型的拱桥上通过时，或是在一拱一拱的正弦曲线上通过时，这轨迹是什么曲线?试画出其图形。程序如下：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);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(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: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);grid on例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(售价与预期销售量);subplot(2,1,2);plot(x2,y2,o);title(广告费与销售增长因子); 首先利用MATLAB编程计算参数a,b,c,d,e。format 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);利用fminsearch函数寻优：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; %初始位置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);hold on %绘制X坐标线xlabel(x);ylabel(y);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，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=CB例12.8 简支梁荷载情况如图12.17所示。求其弯矩、转角和挠度。已知L=8m，q=500N/m，F0=1000N，M0=800Nm，E=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=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); 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=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+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*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,:)=;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的倒数成比例关系。试按照使液体流速最大确定该管道