为了理解EIT正问题利用有限元求解的方法，自己通过在网上查找程序，结合EIDORS软件包所建立的EIT模型，编写了求解节点电势的matlab程序，并且进行了验证。
1.利用EIDORE建立模型，建立模型后，能够得到模型划分的单元、节点坐标等信息。程序如下：
imdl= mk_common_model('b2c',16);%建立二维模型，该函数生成的是inv_model结构。
fmdl= imdl.fwd_model;%赋值正问题模型
fmdl.stimulation= mk_stim_patterns(16,1,[0 1],[0 1], {},1);%设置激励模式
img.fwd_model= fmdl;
img.elem_data=ones(size(fmdl.elems,1),1);
img.elem_data([22:26,36:40])=3; %对介质的电导率信息进行赋值
node_v= calc_all_node_voltages(img);%求解的节点电位值，方便验证
figure;
show_fem(fmdl,[1 1 2]);%显示模型
2.将第一步中所得到的的相关信息进行提取，写入txt文件，为有限元计算提供数据
程序如下：
fid=fopen('D:\matlabanzhuang\bin\FEM\data.txt','wt');%写入文件路径
%%%%存取单元的个数、节点的个数%%%%%
[m,n]=size(fmdl.elems); %m为单元数，q为节点个数
q=size(fmdl.nodes,1);
fprintf(fid,'%d %d\n',q,m);
%%%%%存取每个单元的介质信息--电导率
for i=1:1:m
fprintf(fid,'%d\n',img.elem_data(i));
end
%%%%%存储单元的节点（此处建立的模型为二维，三角形剖分）
for i=1:1:m
for j=1:1:n
if j==n %如果一行的个数达到n个则换行，否则空格
fprintf(fid,'%d\n',fmdl.elems(i,j));
else
fprintf(fid,'%d ',fmdl.elems(i,j));
end
end
end
%%%%存储节点的坐标
[m,n]=size(fmdl.nodes);
for i=1:1:m
for j=1:1:n
if j==n %如果一行的个数达到n个则换行，否则空格
fprintf(fid,'%d\n',fmdl.nodes(i,j));
else
fprintf(fid,'%d ',fmdl.nodes(i,j));
end
end
end
fclose(fid);
3.数据得到后，即可进行有限元求解
程序如下：
%有限元方法求解各节点电势（正问题求解）
%通过相关的软件，得到模型划分的相应单元数，单元的节点数，节点的坐标以及每个单元所对应的节点
%打开数据文件 FP1数据文件指针
%读入初始数据
clear
FP1=fopen('data.txt','rt');
NPOIN=fscanf(FP1,'%d',1);%总节点数
NELEM=fscanf(FP1,'%d',1);%划分的单元数NJIEZHI=fscanf(FP1,'%f',[1,NELEM]);%单元的电导率信息
LNODS=fscanf(FP1,'%f',[3,NELEM])';%单元定义数组，代表每个单元所包含的节点
COORD=fscanf(FP1,'%f',[2,NPOIN])';%结点坐标数组 每个节点的坐标
ASTIF=ASSEMBLE(NPOIN,NELEM,NJIEZHI,COORD,LNODS);%生成总刚度矩阵
%%%%添加边界条件
m=fix(NPOIN/4);%m输入电流节点，自己设定的
n=fix(NPOIN/3);%n输出电流节点，自己设定的，如果极板（输入电流区域）不为点电极，即相应的极板对应的节点均为I/A
%%因为是进行电势的求解，需要设定一个电势参考点，即设定某一点的电势为0.
c=1;%c表示电势为0的节点，此处设定节点一为0电势点
ASTIF(c,1:c-1)=0;
ASTIF(c,c+1:NPOIN)=0;
ASTIF(c,c)=1;%c参考电势为0
b=zeros(NPOIN,1);%b为激励模式，即边界条件
b(114,1)=-1;%设置激励模式，输入节点为1，输出节点电流激励为-1
b(116,1)=1;
%b(m,1)=-1;
%b(n,1)=1;
[Q1,Q2,Q3,Q4]=gaus(ASTIF,b);%Q4即为所求电势解。
scatter(COORD(:,1),COORD(:,2),5,Q4);%散点图
[X,Y,Z]=griddata(COORD(:,1),COORD(:,2),Q4,linspace(-1,1,100)',linspace(-1,1,100),'v4');%将数据进行拟合并践行插值，x、y坐标，电势大小作为z值
figure(1)
contourf(X,Y,Z,20); %等值线图，等势线
figure(2)
surf(X,Y,Z);%三维曲面
figure(3)
scatter(COORD(:,1),COORD(:,2),5,Q4);%散点图
fclose(FP1);%关闭文件
求得各节点的电势，并画出相应的图。
其图像为：
得到的节点电势的值与用EIDORS软件得到的电势值误差不大，说明计算有效。此处激励采用的是1.2电极激励时，各节点的电势大小。
第三步程序中相关函数，可进行下载。