fem.m

给出了线性方程组的CG法、Gauss_Seidel、SOR、Jacobi、fem求解

function [A,b]=fem(u,v,a0,b0)
% a 是x轴上的边界长度 ,b是y轴上的边界长度
% m是x轴上的等分数，n是y轴上的等分数，Nnode为节点总数
% node全局编码数组，为一维数组，存放节点全局编号
% Elematrix单元号与节点局部编码和全局编码的二维数组
% clear all
% clc

% u=32;v=32;a0=1;b0=1;
hx=a0/u;hy=b0/v;
N_node=(u+1)*(v+1); N_elem=2*u*v;
delta_e=hx*hy*0.5;

%生成节点总体编码与各单元的节点局部编码的矩阵
%本程序中的节点编码及其排序是根据《偏微分方程数值解法》（陆金甫，关治）P240中例子 

for e=1:N_elem
    alphax(e)=1;
    alphay(e)=1;
    beta(e)=0;
    t=2*u;
    if mod(e,t)==0
        t=e/t-1;
        x(e)=ceil(e/2)+t;
        M_ju_zong(1,e)=x(e)+1;
        M_ju_zong(2,e)=x(e)+u+2;
        M_ju_zong(3,e)=x(e);
    else
        x(e)=ceil(e/2)+floor(e/t);
        if mod(e,2)==1
            M_ju_zong(1,e)=x(e)+u+1;
            M_ju_zong(2,e)=x(e);
            M_ju_zong(3,e)=x(e)+u+2;
        else
            M_ju_zong(1,e)=x(e)+1;
            M_ju_zong(2,e)=x(e)+u+2;
            M_ju_zong(3,e)=x(e);
        end
    end
end


%初始化总体刚度矩阵
for i=1:N_node
    b(i)=0;
    for j=1:N_node
        A(i,j)=0;
    end
end

%计算ai,aj,am,bi,bj,bm

for e=1:N_elem
      i(e)=M_ju_zong(1,e);
      j(e)=M_ju_zong(2,e);
      m(e)=M_ju_zong(3,e);
%确定i（e）的坐标
      r=ceil(i(e)/(u+1));%第i个节点的行号(取右边整数)
      s=mod(i(e),u+1);   %第i个节点的列号
      if s==0
         s=u+1;
      end
      node(i(e))=(u+1)*(r-1)+s;
      x(i(e))=(s-1)*hx;
      y(i(e))=(r-1)*hy;
%       fe(i(e))=0;
      fe(i(e))=-2*pi^2*sin(pi*x(i(e)))*sin(pi*y(i(e)));
%确定j（e）的坐标
      r=ceil(j(e)/(u+1));%第i个节点的行号(取右边整数)
      s=mod(j(e),u+1);   %第i个节点的列号
      if s==0
         s=u+1;
      end
      node(j(e))=(u+1)*(r-1)+s;
      x(j(e))=(s-1)*hx;
      y(j(e))=(r-1)*hy;
%       fe(i(e))=0;
      fe(j(e))=-2*pi^2*sin(pi*x(j(e)))*sin(pi*y(j(e)));
%确定m（e）的坐标
      r=ceil(m(e)/(u+1));%第i个节点的行号(取右边整数)
      s=mod(m(e),u+1);   %第i个节点的列号
      if s==0
         s=u+1;
      end
      node(m(e))=(u+1)*(r-1)+s;
      x(m(e))=(s-1)*hx;
      y(m(e))=(r-1)*hy;
%       fe(i(e))=0;
      fe(m(e))=-2*pi^2*sin(pi*x(m(e)))*sin(pi*y(m(e)));
      
     be(1)=y(j(e))-y(m(e));
     be(2)=y(m(e))-y(i(e));
     be(3)=y(i(e))-y(j(e));
     ce(1)=x(m(e))-x(j(e));
     ce(2)=x(i(e))-x(m(e));
     ce(3)=x(j(e))-x(i(e));
end
deltae=0.5*(be(1)*ce(2)-be(2)*ce(1));
%计算delta_ij

for i=1:3
    for j=1:3
        if (i~=j)
            del_ij=1.0;
        else
            del_ij=0.0;
        end
    end
end
    
for e=1:N_elem
%计算单元刚度矩阵
     for i=1:3
         b(i)=deltae*fe(m(e))/3;
         for j=1:3
             kk(i,j)=(alphax(e)*be(i)*be(j)+alphay(e)*ce(i)*ce(j))/(4*deltae)+beta(e)*(1.0+del_ij)*deltae/12;
         end
     end
%合成单元刚度矩阵到刚度矩阵
     for i=1:3
        for j=1:3
            A(M_ju_zong(i,e),M_ju_zong(j,e))=A(M_ju_zong(i,e),M_ju_zong(j,e))+kk(i,j);
            b(M_ju_zong(i,e))=b(M_ju_zong(i,e))+b(i);
        end
     end
end

%边界条件的处理

 A(1:u+1,:)=[];%去掉矩阵前u+1行
 A(:,1:u+1)=[];%去掉矩阵前u+1列
 A=A(1:N_node-2*(u+1),:);%去掉矩阵后u+1行
 A=A(:,1:N_node-2*(u+1));%去掉矩阵后u+1列
 b(1:u+1)=[];%去掉右端向量前u+1行
 b(N_node-2*(u+1)+1:N_node-(u+1))=[];%去掉右端向量后u+1行

 %去掉左边界结点所对应的行列
 
 g=N_node-2*(u+1);p=0;
 for i=1:(g-p)
     if mod(i+p,u+1)==1&i+p<=g
         A(i,:)=[];
         A(:,i)=[];
         b(i)=[];
         p=p+1;
     end
 end
%去掉右边界结点所对应的行列

g1=g-v+1;p=0;
 for i=1:(g1-p)
     if mod(i+p,u)==0&i+p<=g1
         A(i,:)=[];
         A(:,i)=[];
         b(i)=[];
         p=p+1;
     end
 end
 b=b';
% size(k);size(b);
% spy(k);
% disp(b');
% k;