(1)
(2)
方程如下
Ne=8;Np=9;
alpha_x=ones(Ne,1)*1;
alpha_y=ones(Ne,1)*1;
beta=zeros(Ne,1);
f=ones(Ne,1)*(-1);
nie=[4 2 1;4 5 2;5 3 2;5 6 3;7 5 4;7 8 5;8 6 5;8 9 6];
x1=2;x2=1;x3=0;
y1=0;y2=1;y3=2;
pos=[x1,y1;x2,y1;x3,y1; ... %postion of points
x1,y2;x2,y2;x3,y2; ...
x1,y3;x2,y3;x3,y3];
b1=pos(nie(:,2),2)-pos(nie(:,3),2);
b2=pos(nie(:,3),2)-pos(nie(:,1),2);
b3=pos(nie(:,1),2)-pos(nie(:,2),2);
c1=pos(nie(:,3),1)-pos(nie(:,2),1);
c2=pos(nie(:,1),1)-pos(nie(:,3),1);
c3=pos(nie(:,2),1)-pos(nie(:,1),1);
s=1/2*(b1.*c2-b2.*c1);
%get K and b
K=zeros(Np,Np);
b=zeros(Np,1);
for ie=1:Ne
Ke=1/4/s(ie)*(alpha_x(ie)*[b1(ie);b2(ie);b3(ie)]*[b1(ie),b2(ie),b3(ie)]+ ...
alpha_y(ie)*[c1(ie);c2(ie);c3(ie)]*[c1(ie),c2(ie),c3(ie)])+...
s(ie)/12*beta(ie)*(1+eye(3));
K(nie(ie,:),nie(ie,:))=K(nie(ie,:),nie(ie,:))+Ke;
be=s(ie)*f(ie)/3;
b(nie(ie,:))= b(nie(ie,:))+be;
end
%apply boundary condition on AB and CD
%phi(1:3)=1;phi(7:9)=0;
K(1:3,:)=0;K(1,1)=1;K(2,2)=1;K(3,3)=1;b(1:3)=1;
K(7:9,:)=0;K(7,7)=1;K(8,8)=1;K(9,9)=1;b(7:9)=0;
%resolve
phi=K^(-1)*b;