文档介绍：程序执行方法：(just for mat lab new users) 7. 1以上版本
打开matlab,点开newM-file,将上述源程序复制粘贴到M-文件中，修改蓝色部分的格式, 保存。按F5即可执行~
%%%% A 99 LIl*Ue'*KE*Ue;%计算目标函数柔度的值
de (ely, elx)二-penal*x(ely, elx)" (penalT)*Ue‘*KE*Ue; %灵敏度分析的结果，参 考文献公式
end
end % FILTERING OF SENSITIVITIES
[de] = check (nelx, nely, rmin, x, de) ; %灵敏度过滤，为了边界光顺
% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
[x] = 0C(nelx, nely, x, volfrac, de) ; %优化准则法更新设计变量
% PRINT RESULTS
change = max (max (abs (x~xold))) ; %计算目标函数的该变量
disp([' It. : ' sprintf C %4i , loop) ' Obj. : ' sprintf ('%10. 4f', c)...
'Vol. : ' sprintf ('%6. 3f', sum (sum (x))/(nelx*nely)) ...
'ch. : sprintf C %6. 3f，, change )]) %屏幕显示迭代信息
% PLOT DENSITIES
colormap(gray); imagesc(~x); axis equal; axis tight; axis off;pause (le~6); end
%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xnew]=0C (nelx, nely, x, volfrac, de) %oc 算法子程序
11 = 0; 12 = 100000; move=0. 2;%11、12用于体积约束的拉格朗日乘子
while (12~11 > le-4)
Imid = 0. 5*(12+11);
xnew 二 max (0. 001, max(x-move, min (1., min (x+move, x. *sqrt (~dc. /Imid))))) ;%核£'部分 参见林文
if sum (sum (xnew)) - volfrac*nelx*nely > 0; % 采用 了二乘法更新拉格朗日乘子
二 Imid;
else
二 Imid;
end
end
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [den] =check(nelx, nely, rmin, x, de) %敏度过滤技术子程序
dcn=zeros(nely, nelx);
for i = 1:nelx
for j = 1:nely
sum=0. 0;
for k