1、工程结构优化设计编程作业1.用黄金分割法求方程 在区间-1,1上的解。33+1=0在 MATLAB 中没有专门的函数实现黄金分割法求解线性方程，可通过编写 glodf.m 函数实现黄金分割法求解，其代码如下:function x=glodf(f,a,b,eps)if(nargin=3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(f1=0)x=a;endif(f2=0)x=b;endif(f1*f20)disp(两端点函数值乘积大于0!);return;elset1=a+(b
2、-a)*0.382;t2=a+(b-a)*0.618;f_1=subs(sym(f),findsym(sym(f),t1);f_2=subs(sym(f),findsym(sym(f),t2);tol=abs(t1-t2);while(toleps)if(f_1*f_20)a=t2;elseb=t1;endendt1=a+(b-a)*0.382;t2=a+(b-a)*0.618;f_1=subs(sym(f),findsym(sym(f),t1);f_2=subs(sym(f),findsym(sym(f),t2);tol=abs(t2-t1);endx=(t1+t2)/2;end其实现的 M
3、ATLAB 代码如下: clear all; tic x=glodf(x3-3*x+1,-1,1,1e-6)x =0.34732.用抛物线法求解 的极值点，取值范围0.5，1，允许误差为 1e-6。=3+20.2在 MATLAB 中没有专门的函数实现抛物线法求解线性方程，可通过编写 paowuxian.m 函数实现抛物线法求解，其代码如下:function x=paowuxian(f,x0,x1,x2,e)if nargineh1=y-z;h2=x-y;c1=(feval(f,y)-feval(f,z)/h1;c2=(feval(f,x)-feval(f,y)/h2;d=(c1-c2)/(h2
4、+h1);w=c2+h2*d;xi=x-(2*feval(f,x)/(w+(w/abs(w)*sqrt(w2-4*feval(f,x)*d);z=y;y=x;x=xi;end其实现的 MATLAB 代码如下: fun=inline(x3+2*x-0.2); paowuxian(fun,0,0.5,1,1e-6),format shortans =0.09953.用齿形法对下图的结构进行杆件截面优化设计：节点处作用两种工况：（1） P1=2000kN，P2=0kN ；（2）P1=0kN，P2=2000kN。各杆采用同一材料制成，弹性模量 E 为常数，容许应力分别为： +=2107kPa ，-=1
5、.5107kPa。此结构为对称结构，工况也对称，所以 A1=A3，计算时只取两个变量 A1 和 A2，一种工况 P1=2000kN，P2=0kN。目标函数取为结构总重量：在 MATLAB 中没有专门的函数实现齿形法求解线性方程，可通过编写 Zigzag.m 函数实现齿形法求解，其代码如下:function newA1,newA2,newA3,W,exitflag =Zigzag( A1,A2,maxstep)if narginw1newA1=A11;newA2=A22;newA3=A11;W=vpa(w1*p*l,4);exitflag=1;lAW)2(1return;endu1=max(U)
6、;A=A.*u1;A11=A(1);A22=A(2);U=U./u1;U=roundn(U,-4);w=sqrt(2)*2*A11+A22;w1=w;if w1w2newA1=A13;newA2=A24;newA3=A13;W=vpa(w2*p*l,4);exitflag=0;return;endA=A.*U;A1=A(1);A2=A(2);U=(1e-4)*(sqrt(2)*A1+A2)/(sqrt(2)*A12+2*A1*A2);(1e-4)*(sqrt(2)*A1)/(sqrt(2)*A12+2*A1*A2);u1=max(U);A=A.*u1;A13=A(1);A24=A(2);U=U
7、./u1;U=roundn(U,-4);w=sqrt(2)*2*A13+A24;w2=w;A=A.*U;A1=A(1);A2=A(2);U=(1e-4)*(sqrt(2)*A1+A2)/(sqrt(2)*A12+2*A1*A2);(1e-4)*(sqrt(2)*A1)/(sqrt(2)*A12+2*A1*A2);endend其实现的MATLAB代码如下: clear;A1=1e-4;A2=1e-4;Anew1,Anew2,Anew3,W,exitflag=Zigzag(A1,A2)Anew1 =7.7346e-05Anew2 =4.5309e-05Anew3 =7.7346e-05W =0.0
8、002641*l*pexitflag =04. 用修改齿形法对结构进行杆件截面优化设计。在 MATLAB 中没有专门的函数实现修改齿形法求解线性方程，可通过编写 IZigzag.m 函数实现修改齿形法求解，其代码如下:function newA1,newA2,newA3,W,exitflag =IZigzag( A1,A2,v,maxstep)if narginw1newA1=A11;newA2=A22;newA3=A11;W=vpa(w1*p*l,4);exitflag=1;return;endu1=max(U);A=A.*u1;A11=A(1);A22=A(2);U=U./u1;x1=U(
9、1);x2=U(2);y1=1-v*(1-x1);y2=1-v*(1-x2);U=y1;y2;U=roundn(U,-4);w=sqrt(2)*2*A11+A22;w1=w;if w1w2newA1=A13;newA2=A24;newA3=A13;W=vpa(w2*p*l,4);exitflag=0;return;endA=A.*U;A1=A(1);A2=A(2);U=(1e-4)*(sqrt(2)*A1+A2)/(sqrt(2)*A12+2*A1*A2);(1e-4)*(sqrt(2)*A1)/(sqrt(2)*A12+2*A1*A2);u1=max(U);A=A.*u1;A13=A(1);
10、A24=A(2);U=U./u1;x1=U(1);x2=U(2);y1=1-v*(1-x1);y2=1-v*(1-x2);U=y1;y2;U=roundn(U,-4);w=sqrt(2)*2*A13+A24;w2=w;A=A.*U;A1=A(1);A2=A(2);U=(1e-4)*(sqrt(2)*A1+A2)/(sqrt(2)*A12+2*A1*A2);(1e-4)*(sqrt(2)*A1)/(sqrt(2)*A12+2*A1*A2);endend其实现的 MATLAB 代码如下: clear;A1=1e-4;A2=1e-4;v=0.5;Anew1,Anew2,Anew3,W,exitfla
11、g=IZigzag(A1,A2,v)Anew1 =7.8168e-05Anew2 =4.2839e-05Anew3 =7.8168e-05W =0.0002639*l*pexitflag =05. 用满应力法对结构进行杆件截面优化设计在 MATLAB 中没有专门的函数实现满应力法求解线性方程，可通过编写Full_Stressed.m 函数实现满应力法求解，其代码如下:function newA1,newA2,newA3,exitflag =Full_Stressed( A1,A2,wucha,maxstep)if nargin clear all;close all;A1=1e-4;A2=1e
12、-4;newA1,newA2,newA3,exitflag=Full_Stressed (A1,A2)newA1 =9.8988e-05newA2 =1.4373e-06newA3 =9.8988e-05exitflag =-16. 用加松弛因子满应力法对结构进行杆件截面优化设计在 MATLAB 中没有专门的函数实现加松弛因子满应力法求解线性方程，可通过编写SRF_Full_Stressed.m 函数实现满应力法求解，其代码如下:function newA1,newA2,newA3,exitflag =SRF_Full_Stressed( A1,A2,B,wucha,maxstep)if nargin clear all;close all;A1=1e-4;A2=1e-4;B=1.1;wucha=0.05;newA1,newA2,newA3,exitflag=SRF_Full_Stressed (A1,A2,B,wucha)newA1 =9.5527e-05newA2 =6.3262e-06newA3 =9.5527e-05
