有限元法求解泊松方程，用matlab编写的计算程序

是基于MATLAB的有限元磁场分析，泊松方程为别接条件的，

% 用有限元法求解三角形形区域上的Possion方程
function Finite_element_tri（ImaxJmax）
global ndm nel na
% ndm 总节点数
% nel 基元数
% na 活动节点数
Imax=30;Jmax=60;%设定网格数
V=0; J=0;X0=1/Imax;Y0=X0;
domain_tri   
[XYNNNE]=setelm_tri（ImaxJmax）;  % 给节点和三角形元素编号，并设定节点坐标
% 求解有限元方程的求系数矩阵
T=zeros（ndmndm）;
for n=1:nel
    n1=NE（1n）;   n2=NE（2n）;   n3=NE（3n）;
    s=abs（（X（n2）-X（n1））*（Y（n3）-Y（n1））-（X（n3）-X（n1））*（Y（n2）-Y（n1）））/2;
    for k=1:3
        if n1<=na|n2<=na
        T（n1n2）=T（n1n2）+（（Y（n2）-Y（n3））*（Y（n3）-Y（n1））+（X（n3）-X（n2））*（X（n1）-X（n3）））/（4*s）;
        T（n2n1）=T（n1n2）;
        T（n1n1）=T（n1n1）+（（Y（n2）-Y（n3））^2+（X（n3）-X（n2））^2）/（4*s）;
        end
        k=n1;n1=n2;n2=n3;n3=k; 
    end
end
M=T（1:na1:na）;
% 求有限元方程的右端项
G=zeros（na1）;
for n=1:nel
    n1=NE（1n）;   n2=NE（2n）;   n3=NE（3n）;  
    s=abs（（X（n2）-X（n1））*（Y（n3）-Y（n1））-（X（n3）-X（n1））*（Y（n2）-Y（n1）））/2;
    for k=1:3
        if n1<=na
            G（n1）=G（n1）+（2*X（n1）+X（n2）+X（n3））*s/12;
        end
        n4=n1;  n1=n2;  n2=n3;  n3=n4; 
    end
end

%  求解方程得结果
F=M\G; 
NNV=zeros（Imax+1Jmax+1）;
fi=zeros（ndm1）;
fi（1:na）=F（1:na）;
fi（na+1:ndm）=V;
for j=0:Jmax
    for i=0:Imax
        n=NN（i+1j+1）;
        if n<=0
            n=na+1;
        end
        NNV（i+1j+1）=fi（n）;
    end
end

% 画等电势线
X1=zeros（1Imax+1）;
Y1=zeros（1Jmax+1）;
for i=1:Imax+1
    X1（i）=（i-1）*X0;
end
for i=1:Jmax+1
    Y1（i）=（i-1）*Y0;
end

% 画解函数的曲面图
figure（2）
surf（X1Y1NNV‘）;
fid=fopen（‘Finite_element_tri.txt‘‘w‘）;
fprintf（fid‘\n *********有限元法求解三角形区域上Possion方程的结果********** \n \n‘）;
fprintf（fid‘\n 节点编号 \n \n‘）;
Nna=fliplr（NN）;
fprintf（fid‘%4d%4d%4d%4d%4d%4d%4d%4d%4d\n‘Nna）;fprintf（fid‘\n 各节点的电势 \n \n‘）;
NNV=fliplr（NNV）;
fprintf（fid‘%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f\n‘NNV）;
L=[1:ndm]‘;
fprintf（fid‘\n\n    节点编号    坐标分量x   坐标分量y      u（xy）的值\n\n‘）;
for i=1:ndm
    fprintf（fid‘%8d%14.5f%14.5f%14.5f\n‘L（i）X（i）Y（i）fi（i））;
end
fclose（fid）;
end 
 
function [XYNNNE]=setelm_tri（ImaxJmax）
% 给节点和三角形元素编号，并设定节点坐标
global  ndm nel na
% I1 I2 J1 J2 Imax Jmax分别描述网线纵向和横向数目的变量
% X Y表示节点坐标
% NN描述节点编号
% NE 描述各基点局域节点的矩阵
% ndm 总节点数
% nel 基元数
% na 表示活动节点数
nlm=Imax*Jmax;
dx=1/Imax;
dy=1/Jmax;
X=nlm:1;
Y=nlm:1;
NN=zeros（Imax+1Jmax+1）;
n1=0; 
% 活动节点编号
for j=3:Jmax/2
    for i=2:j-1
        n1=n1+1;
        NN（ij）=n1;
        X（n1）=（i-1）*dx;
        Y（n1）=-1+（j-1）*dy;
    end
end
k=Jmax/2+1;
for j=Jmax/2+1:Jmax-1
    k=k-1;
    for i=2:k
        n1=n1+1;
        NN（ij）=n1;
        X（n1）=（i-1）*dx;
        Y（n1）=1+（j-Jmax-1）*dy;
    end
end
na=n1;
for j=Jmax+1:-1:Jmax/2+1
    n1=n1+1;
    NN（1j）=n1;
    X（n1）=0;
    Y（n1）=1+（j-Jmax-1）*dy;
end
for j=Jmax/2:-1:1
    n1=n1+1;
    NN（1j）=n1;
    X（n1）=0;
    Y（n1）=-1+（j-1）*dy;
end
for i=2:Imax+1
    n1=n1+1;
    NN（ii）=n1;
    X（n1）=（i-1）*dx;
    Y（n1）=-1+（i-1）*dy;
end
K=0;
for i=Imax:-1:2
    K=K+2;
    n1=n1+1;
    NN（ii+K）=n1;
    X（n1）=（i-

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       4717  2007-01-13 21:11  Finite_element_tri.m

----------- ---------  ---------- -----  ----

                 4935                    2