matlab程序的有限元法解泊松方程

1、.基于matlab編程的有限元法一、待求問題：泛定方程：邊界條件：以（0,-1），（0,1），（1,0）為頂點的三角形區(qū)域邊界上2、 編程思路及方法1、 給節(jié)點和三角形單元編號，并設定節(jié)點坐標畫出以（0,-1），（0,1），（1,0）為頂點的三角形區(qū)域figure1由于積分區(qū)域規(guī)則，故采用特殊剖分單元，將區(qū)域沿水平豎直方向分等份，此時所有單元都是等腰直角三角形，剖分單元個數由自己輸入，但豎直方向份數（用jmax表示）必須是水平方向份數（imax）的兩倍，所以用戶只需輸入水平方向的份數imax。采用上述剖分方法，節(jié)點位置也比較規(guī)則。然后利用循環(huán)從區(qū)域內部（非邊界）的節(jié)點開始編號，格式為nn(i,

2、j)=n1，i，j分別表示節(jié)點所在列數與行數，并將節(jié)點坐標存入相應矩陣x(n1)，y(n1)。由于區(qū)域上下兩部分形狀不同因此，分兩個循環(huán)分別編號賦值，然后再對邊界節(jié)點編號賦值。然后再每個單元的節(jié)點進行局部編號，由于求解區(qū)域和剖分單元的特殊性，分別對內部節(jié)點對應左上角正方形的兩個三角形單元，上左，左上，下斜邊界節(jié)點要對應三個單元，上左，左上，左下，右頂點的左下、左上，右上邊界的左上，分別編號以保證覆蓋整個區(qū)域。2、 求解泊松方程首先一次獲得每個單元節(jié)點的整體編號，然后根據其坐標求出每個三角形單元的面積。利用有限元方法的原理，分別求出系數矩陣和右端項，并且由于邊界條件特殊，邊界上，因此做積分時只需

3、對場域單元積分而不必對邊界單元積分。求的兩個矩陣后很容易得到節(jié)點電位向量，即泊松方程的解。3、 畫解函數的平面圖和曲面圖由節(jié)點單位向量得到，j行i列節(jié)點的電位，然后調用繪圖函數imagesc(nnv)與surf(x1,y1,nnv)分別得到解函數的平面圖figure2和曲面圖figure3。4、 將結果輸出為文本文件輸出節(jié)點編號，坐標，電位值3、 計算結果1、積分區(qū)域：2、f=1，x方向75份，y方向150份時，解函數平面圖和曲面圖對比：當f=1時，界函數平面圖3、輸出文本文件由于節(jié)點多較大，列在本文最末四、結果簡析由于三角形區(qū)域分布的是正電荷，因此必定電位最高點在區(qū)域中部，且沿x軸對稱，三角

4、形邊界電位最低等于零。當f=1時，發(fā)現(xiàn)電位最高點向x軸負方向移動了，這是由于此時電荷在三角形區(qū)域上均勻分布，而f=x時，x越大面電荷密度越大，附近相應電位越高，所得圖像與實際情況相符。五、matlab源程序1、finite_element_tri.m文件function finite_element_tri(imax)% 用有限元法求解三角形形區(qū)域上的possion方程,其中a為1，f=xjmax=2*imax;% 其中imax jmax分別表示x軸和y軸方向的網格數，其中jmax等于imax的兩倍% 定義一些全局變量global ndm nel na % ndm 總節(jié)點數% nel 基元數%

5、 na 表示活動節(jié)點數v=0; j=0;x0=1/imax;y0=x0;%v=0為邊界條件domain_tri % 調用函數畫求解區(qū)域x,y,nn,ne=setelm_tri(imax,jmax); % 給節(jié)點和三角形元素編號，并設定節(jié)點坐標% 以下求解有限元方程的求系數矩陣t=zeros(ndm,ndm);for n=1:nel n1=ne(1,n); n2=ne(2,n); n3=ne(3,n);%整體編號 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

6、(n1,n2)=t(n1,n2)+(y(n2)-y(n3)*(y(n3)-y(n1)+(x(n3)-x(n2)*(x(n1)-x(n3)/(4*s); t(n2,n1)=t(n1,n2); t(n1,n1)=t(n1,n1)+(y(n2)-y(n3)2+(x(n3)-x(n2)2)/(4*s);%v=0則邊界積分為零，非零時積分編程類似，再加邊界積分。 end k=n1;n1=n2;n2=n3;n3=k; % 輪換坐標將值賦入3階主子矩陣中 endendm=t(1:na,1:na);% 求有限元方程的右端項f=x;%場源函數g=zeros(na,1);for n=1:nel n1=ne(1,n

7、); n2=ne(2,n); n3=ne(3,n); 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*f(n1)+f(n2)+f(n3)*s/12;%f在單元上為線性差值時場域單元的積分公式 end n4=n1; n1=n2; n2=n3; n3=n4; % 輪換坐標標 endendf=mg; % 求解方程得結果nnv=zeros(imax+1,jmax+1);fi=zeros(ndm,1);fi(1:na)=f(1:na);fi(na+1:ndm)=v;

8、for j=0:jmax for i=0:imax n=nn(i+1,j+1); if n=0 n=na+1; end nnv(i+1,j+1)=fi(n); endendfigure(2)imagesc(nnv);%畫解函數的平面圖x1=zeros(1,imax+1);y1=zeros(1,jmax+1);for i=1:imax+1 x1(i)=(i-1)*x0;endfor i=1:jmax+1 y1(i)=(i-1)*y0;endfigure(3)surf(x1,y1,nnv);% 畫解函數的曲面圖% 以下是結果的輸出fid=fopen(finite_element_tri.txt,w

9、);fprintf(fid,n *有限元法求解三角形區(qū)域上possion方程的結果* n n);l=1:ndm;fprintf(fid,nn 節(jié)點編號 坐標分量x 坐標分量y u（x,y）的值nn);for i=1:ndm fprintf(fid,%8d%14.5f%14.5f%14.5fn,l(i),x(i),y(i),fi(i);endfclose(fid);end2、domain_tri.m文件function domain_tri% 畫求解區(qū)域xy=0 1;0 -1;1 0;a=zeros(3,3);a(1,1)=2; a(1,2)=-1;a(1,3)=-1;a(2,2)=2; a(2

10、,1)=-1;a(2,3)=-1;a(3,3)=2; a(3,2)=-1;a(3,1)=-1;a=sparse(a);figure(1);gplot(a,xy);3、setelm_tri.m文件function x,y,nn,ne=setelm_tri(imax,jmax) % 給節(jié)點和三角形單元編號，并設定節(jié)點坐標 % 定義一些全局變量global ndm nel na% i1 i2 j1 j2 imax jmax分別描述網線縱向和橫向數目的變量% x y表示節(jié)點坐標% nn描述節(jié)點編號% ne 描述各單元局部節(jié)點編號與總體編號對應的矩陣% ndm 總節(jié)點數% nel 單元數% na 表示不

11、含邊界的節(jié)點數nlm=imax*jmax;dx=1/imax;dy=1/jmax;x=nlm:1;y=nlm:1;nn=zeros(imax+1,jmax+1);n1=0; for j=3:jmax/2 for i=2:j-1 n1=n1+1; nn(i,j)=n1; %x=i列，y=j行處節(jié)點 x(n1)=(i-1)*dx; y(n1)=-1+(j-1)*dy; endendk=jmax/2+1;for j=jmax/2+1:jmax-1 %三角形區(qū)域上下兩部分節(jié)點坐標分別求 k=k-1; for i=2:k n1=n1+1; nn(i,j)=n1; x(n1)=(i-1)*dx; y(n1

12、)=1+(j-jmax-1)*dy; endendna=n1;%不含邊界節(jié)點數for j=jmax+1:-1:jmax/2+1 %降序 n1=n1+1; nn(1,j)=n1; x(n1)=0; y(n1)=1+(j-jmax-1)*dy;endfor j=jmax/2:-1:1 n1=n1+1; nn(1,j)=n1; x(n1)=0; y(n1)=-1+(j-1)*dy;end %for i=2:imax+1 n1=n1+1; nn(i,i)=n1; x(n1)=(i-1)*dx; y(n1)=-1+(i-1)*dy;endk=0;for i=imax:-1:2 k=k+2; n1=n1+

13、1; nn(i,i+k)=n1; x(n1)=(i-1)*dx; y(n1)=1+(i+k-jmax-1)*dy;end% 以上四個循環(huán)為對邊界節(jié)點進行編號ndm=n1; ne=zeros(3,2*ndm); n1=0;for j=3:jmax/2 for i=2:j-1 n1=n1+1; ne(1,n1)=nn(i,j); ne(2,n1)=nn(i-1,j+1); ne(3,n1)=nn(i-1,j); n1=n1+1; ne(1,n1)=nn(i,j); ne(2,n1)=nn(i,j+1); ne(3,n1)=nn(i-1,j+1); endendk=jmax/2+1;for j=jm

14、ax/2+1:jmax-1 k=k-1; for i=2:k n1=n1+1; ne(1,n1)=nn(i,j); ne(2,n1)=nn(i-1,j+1); ne(3,n1)=nn(i-1,j); n1=n1+1; ne(1,n1)=nn(i,j); ne(2,n1)=nn(i,j+1); ne(3,n1)=nn(i-1,j+1); endend %內部節(jié)點對應左上角正方形的兩個三角形單元，上左，左上for i=2:imax n1=n1+1; ne(1,n1)=nn(i,i); ne(2,n1)=nn(i-1,i); ne(3,n1)=nn(i-1,i-1); n1=n1+1; ne(1,n

15、1)=nn(i,i); ne(2,n1)=nn(i-1,i+1); ne(3,n1)=nn(i-1,i); n1=n1+1; ne(1,n1)=nn(i,i); ne(2,n1)=nn(i,i+1); ne(3,n1)=nn(i-1,i+1);end %下斜邊界節(jié)點要對應三個單元，上左，左上，左下n1=n1+1;ne(1,n1)=nn(imax+1,imax+1);ne(2,n1)=nn(imax,imax+1);ne(3,n1)=nn(imax,imax);%右頂點的左下n1=n1+1;ne(1,n1)=nn(imax+1,imax+1);ne(2,n1)=nn(imax,imax+2);n

16、e(3,n1)=nn(imax,imax+1);%右頂點的左上k=0;for i=imax:-1:2 k=k+2; n1=n1+1; ne(1,n1)=nn(i,i+k); ne(2,n1)=nn(i-1,i+k+1); ne(3,n1)=nn(i-1,i+k);end %右上邊界的左上nel=n1;%此時n1值為總的單元個數六、輸出文本文件 *有限元法求解三角形區(qū)域上possion方程的結果* 節(jié)點編號 坐標分量x 坐標分量y u（x,y）的值 1 0.01333 -0.98667 0.00000 2 0.01333 -0.98000 0.00001 3 0.02667 -0.98000 0

17、.00001 4 0.01333 -0.97333 0.00002 5 0.02667 -0.97333 0.00003 6 0.04000 -0.97333 0.00002 7 0.01333 -0.96667 0.00002 8 0.02667 -0.96667 0.00004 9 0.04000 -0.96667 0.00005 10 0.05333 -0.96667 0.00003 11 0.01333 -0.96000 0.00003 12 0.02667 -0.96000 0.00006 13 0.04000 -0.96000 0.00007 14 0.05333 -0.96000

18、 0.00007 15 0.06667 -0.96000 0.00005 16 0.01333 -0.95333 0.00004 17 0.02667 -0.95333 0.00008 18 0.04000 -0.95333 0.00010 19 0.05333 -0.95333 0.00011 20 0.06667 -0.95333 0.00010 21 0.08000 -0.95333 0.00007 22 0.01333 -0.94667 0.00005 23 0.02667 -0.94667 0.00010 24 0.04000 -0.94667 0.00014 25 0.05333

19、-0.94667 0.00016 26 0.06667 -0.94667 0.00016 27 0.08000 -0.94667 0.00013 28 0.09333 -0.94667 0.00008 29 0.01333 -0.94000 0.00006 30 0.02667 -0.94000 0.00012 31 0.04000 -0.94000 0.00017 32 0.05333 -0.94000 0.00020 33 0.06667 -0.94000 0.00022 34 0.08000 -0.94000 0.00021 35 0.09333 -0.94000 0.00017 36

20、0.10667 -0.94000 0.00010 37 0.01333 -0.93333 0.00007 38 0.02667 -0.93333 0.00015 39 0.04000 -0.93333 0.00021 40 0.05333 -0.93333 0.00025 41 0.06667 -0.93333 0.00028 42 0.08000 -0.93333 0.00028 43 0.09333 -0.93333 0.00026 44 0.10667 -0.93333 0.00021 45 0.12000 -0.93333 0.00012 46 0.01333 -0.92667 0.0

21、0009 47 0.02667 -0.92667 0.00017 48 0.04000 -0.92667 0.00024 49 0.05333 -0.92667 0.00030 50 0.06667 -0.92667 0.00034 51 0.08000 -0.92667 0.00036 52 0.09333 -0.92667 0.00036 53 0.10667 -0.92667 0.00032 54 0.12000 -0.92667 0.00025 55 0.13333 -0.92667 0.00015 56 0.01333 -0.92000 0.00010 57 0.02667 -0.9

22、2000 0.00020 58 0.04000 -0.92000 0.00028 59 0.05333 -0.92000 0.00036 60 0.06667 -0.92000 0.00041 61 0.08000 -0.92000 0.00045 62 0.09333 -0.92000 0.00046 63 0.10667 -0.92000 0.00044 64 0.12000 -0.92000 0.00038 65 0.13333 -0.92000 0.00030 66 0.14667 -0.92000 0.00017 67 0.01333 -0.91333 0.00011 68 0.02

23、667 -0.91333 0.00022 69 0.04000 -0.91333 0.00032 70 0.05333 -0.91333 0.00041 71 0.06667 -0.91333 0.00048 72 0.08000 -0.91333 0.00053 73 0.09333 -0.91333 0.00056 74 0.10667 -0.91333 0.00056 75 0.12000 -0.91333 0.00052 76 0.13333 -0.91333 0.00045 77 0.14667 -0.91333 0.00034 78 0.16000 -0.91333 0.00019

24、 79 0.01333 -0.90667 0.00013 80 0.02667 -0.90667 0.00025 81 0.04000 -0.90667 0.00037 82 0.05333 -0.90667 0.00047 83 0.06667 -0.90667 0.00055 84 0.08000 -0.90667 0.00062 85 0.09333 -0.90667 0.00066 86 0.10667 -0.90667 0.00068 87 0.12000 -0.90667 0.00066 88 0.13333 -0.90667 0.00061 89 0.14667 -0.90667

25、 0.00052 90 0.16000 -0.90667 0.00039 91 0.17333 -0.90667 0.00022 92 0.01333 -0.90000 0.00014 93 0.02667 -0.90000 0.00028 94 0.04000 -0.90000 0.00041 95 0.05333 -0.90000 0.00053 96 0.06667 -0.90000 0.00063 97 0.08000 -0.90000 0.00071 98 0.09333 -0.90000 0.00077 99 0.10667 -0.90000 0.00080 100 0.12000

26、 -0.90000 0.00080 101 0.13333 -0.90000 0.00077 102 0.14667 -0.90000 0.00070 103 0.16000 -0.90000 0.00059 104 0.17333 -0.90000 0.00044 105 0.18667 -0.90000 0.00024 106 0.01333 -0.89333 0.00016 107 0.02667 -0.89333 0.00031 108 0.04000 -0.89333 0.00045 109 0.05333 -0.89333 0.00059 110 0.06667 -0.89333

27、0.00071 111 0.08000 -0.89333 0.00080 112 0.09333 -0.89333 0.00088 113 0.10667 -0.89333 0.00093 114 0.12000 -0.89333 0.00095 115 0.13333 -0.89333 0.00093 116 0.14667 -0.89333 0.00089 117 0.16000 -0.89333 0.00080 118 0.17333 -0.89333 0.00067 119 0.18667 -0.89333 0.00049 120 0.20000 -0.89333 0.00027 12

28、1 0.01333 -0.88667 0.00017 122 0.02667 -0.88667 0.00034 123 0.04000 -0.88667 0.00050 124 0.05333 -0.88667 0.00065 125 0.06667 -0.88667 0.00078 126 0.08000 -0.88667 0.00090 127 0.09333 -0.88667 0.00099 128 0.10667 -0.88667 0.00106 129 0.12000 -0.88667 0.00110 130 0.13333 -0.88667 0.00110 131 0.14667

29、-0.88667 0.00107 132 0.16000 -0.88667 0.00101 133 0.17333 -0.88667 0.00090 134 0.18667 -0.88667 0.00074 135 0.20000 -0.88667 0.00055 136 0.21333 -0.88667 0.00030 137 0.01333 -0.88000 0.00019 138 0.02667 -0.88000 0.00037 139 0.04000 -0.88000 0.00055 140 0.05333 -0.88000 0.00071 141 0.06667 -0.88000 0

30、.00086 142 0.08000 -0.88000 0.00100 143 0.09333 -0.88000 0.00111 144 0.10667 -0.88000 0.00119 145 0.12000 -0.88000 0.00125 146 0.13333 -0.88000 0.00127 147 0.14667 -0.88000 0.00126 148 0.16000 -0.88000 0.00122 149 0.17333 -0.88000 0.00113 150 0.18667 -0.88000 0.00100 151 0.20000 -0.88000 0.00082 152

31、 0.21333 -0.88000 0.00060 153 0.22667 -0.88000 0.00033 154 0.01333 -0.87333 0.00020 155 0.02667 -0.87333 0.00040 156 0.04000 -0.87333 0.00060 157 0.05333 -0.87333 0.00078 158 0.06667 -0.87333 0.00095 159 0.08000 -0.87333 0.00109 160 0.09333 -0.87333 0.00122 161 0.10667 -0.87333 0.00133 162 0.12000 -

32、0.87333 0.00140 163 0.13333 -0.87333 0.00145 164 0.14667 -0.87333 0.00146 165 0.16000 -0.87333 0.00143 166 0.17333 -0.87333 0.00137 167 0.18667 -0.87333 0.00126 168 0.20000 -0.87333 0.00111 169 0.21333 -0.87333 0.00091 170 0.22667 -0.87333 0.00066 171 0.24000 -0.87333 0.00036 172 0.01333 -0.86667 0.

33、00022 173 0.02667 -0.86667 0.00044 174 0.04000 -0.86667 0.00065 175 0.05333 -0.86667 0.00084 176 0.06667 -0.86667 0.00103 177 0.08000 -0.86667 0.00120 178 0.09333 -0.86667 0.00134 179 0.10667 -0.86667 0.00146 180 0.12000 -0.86667 0.00156 181 0.13333 -0.86667 0.00162 182 0.14667 -0.86667 0.00165 183

34、0.16000 -0.86667 0.00165 184 0.17333 -0.86667 0.00161 185 0.18667 -0.86667 0.00152 186 0.20000 -0.86667 0.00139 187 0.21333 -0.86667 0.00122 188 0.22667 -0.86667 0.00099 189 0.24000 -0.86667 0.00071 190 0.25333 -0.86667 0.00038 191 0.01333 -0.86000 0.00024 192 0.02667 -0.86000 0.00047 193 0.04000 -0

35、.86000 0.00070 194 0.05333 -0.86000 0.00091 195 0.06667 -0.86000 0.00111 196 0.08000 -0.86000 0.00130 197 0.09333 -0.86000 0.00146 198 0.10667 -0.86000 0.00160 199 0.12000 -0.86000 0.00172 200 0.13333 -0.86000 0.00180 201 0.14667 -0.86000 0.00185 202 0.16000 -0.86000 0.00187 203 0.17333 -0.86000 0.0

36、0185 204 0.18667 -0.86000 0.00179 205 0.20000 -0.86000 0.00168 206 0.21333 -0.86000 0.00153 207 0.22667 -0.86000 0.00133 208 0.24000 -0.86000 0.00108 209 0.25333 -0.86000 0.00077 210 0.26667 -0.86000 0.00041 211 0.01333 -0.85333 0.00025 212 0.02667 -0.85333 0.00050 213 0.04000 -0.85333 0.00075 214 0

37、.05333 -0.85333 0.00098 215 0.06667 -0.85333 0.00120 216 0.08000 -0.85333 0.00140 217 0.09333 -0.85333 0.00158 218 0.10667 -0.85333 0.00174 219 0.12000 -0.85333 0.00188 220 0.13333 -0.85333 0.00198 221 0.14667 -0.85333 0.00205 222 0.16000 -0.85333 0.00209 223 0.17333 -0.85333 0.00209 224 0.18667 -0.85333 0.00205 225 0.20000 -0.85333 0.00197 226 0.21333 -0.85333 0.00184 227 0