xd=input('');

xn=size(xd,1);

AA=1.0e30;                        % AA为大数

for i=1:xn

    K1(2*xd(i,1)-1,2*xd(i,1)-1)=K1(2*xd(i,1)-1,2*xd(i,1)-1)*AA;

    F(2*xd(i,1)-1)=xd(i,2)*K1(2*xd(i,1)-1,2*xd(i,1)-1);

end

fprintf('请输入y方向受到位移约束的结点编号及其约束位移值(形式为[node1 v1;node2 v2;……;noden vn]) \n ');

yd=input('');

yn=size(yd,1);

for i=1:yn

    K1(2*yd(i,1),2*yd(i,1))=K1(2*yd(i,1),2*yd(i,1))*AA;

    F(2*yd(i,1))=yd(i,2)*K1(2*yd(i,1),2*yd(i,1));

end

u=K1\F;                            % 解整体刚度方程KU=P

fprintf('边界修正后的整体刚度矩阵为K1=\n');

disp(K1);

fprintf('边界修正后的外荷载列为F=\n');

disp(F);

fprintf('在外荷载作用下各结点的位移为u=\n');

disp(u);

〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓

function k=k(co,E,NU,D,hou)

% co为单元的结点坐标，为8行2列矩阵

% 函数返回单元的刚度矩阵

syms s t;         % s,t表示局部坐标

N=sym(zeros(8,1));

N(1)=(1-s)*(1-t)*(-s-t-1)/4;

N(2)=(1+s)*(1-t)*(s-t-1)/4;

N(3)=(1+s)*(1+t)*(s+t-1)/4;

N(4)=(1-s)*(1+t)*(-s+t-1)/4;

N(5)=(1-t)*(1+s)*(1-s)/2;

N(6)=(1+s)*(1+t)*(1-t)/2;

N(7)=(1+t)*(1+s)*(1-s)/2;

N(8)=(1-s)*(1+t)*(1-t)/2;      

x=0;

y=0;

for i=1:8

    x=x+N(i)*co(i,1);

    y=y+N(i)*co(i,2);

end