注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

小明砸乱弹琴

一个工程师的日常

 
 
 

日志

 
 
关于我

喜欢哲学、几何学、心理学、物理学、音乐、简单的文字及一切有趣的事物。热爱幻想,追求思想和学术自由,热衷简单平凡的生活和发现美好的事物。

网易考拉推荐
 
 

基于Matlab的有限元方法求解泊松方程  

2013-05-21 21:01:12|  分类: 编程语言 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

这只是一个作业,提供一个例子,供大家参考:

 
基于Matlab的有限元方法求解泊松方程 - 赵明 - Alexander
 

 

方法步骤:

一、             求出方程的变分形式

二、             区域剖分

三、             构造有限元子空间

四、             建立有限元方程

五、             边界条件处理

六、             求解代数方程组,输出刚度矩阵

七、             输出计算结果


输出矩阵:

基于Matlab的有限元方法求解泊松方程 - 赵明 - Alexander



输出图像


基于Matlab的有限元方法求解泊松方程 - 赵明 - Alexander

程序部分:

clear;%清除所有spacework里面的数据

clear;

 

Xn=1;Yn=1;%[1X1]的空间范围

N=4;

h=1/N;

 

np=[1      7     2     8     3     9     4     10    6     12    7     13    8     14    9     15 11      17    12    18    13    19    14    20    16    22       17    23    18    24    19    25;

    2     6     3     7     4     8     5     9     7     11    8     12    9     13    10    14 12      16    13    17    14    18    15    19    17    21       18    22    19    23    20    24;

    6     2     7     3     8     4     9     5     11    7     12    8     13    9     14    10 16      12    17    13    18    14    19    15    21    17       22    18    23    19    24    20];%局部编码

NE=[1,2,3,4,5,6,11,16,21];%9个第一类节点

NF=[7,8,9,10,12,13,14,15,17,18,19,20,22,23,24,25];%16个第二类节点

 

l=[6,6,6,3,6,6,6,3,6,6,6,3,3,3,3,1];%每一个节点影响元素的个数

 

m=[7,7,7,5,7,7,7,5,7,7,7,5,5,5,5,3];%每一个节点影响节点的个数

Sbe7=[2 3 4 9 10 11];

Sbe8=[4 5 6 11 12 13];

Sbe9=[6 7 8 13 14 15];

Sbe10=[8 15 16];

Sbe12=[10 11 12 17 18 19];

Sbe13=[12 13 14 19 20 21];

Sbe14=[14 15 16 21 22 23];

Sbe15=[16 23 24];

Sbe17=[18 19 20 25 26 27];

Sbe18=[20 21 22 27 28 29];

Sbe19=[22 23 24 29 30 31];

Sbe20=[24 31 32];

Sbe22=[26 27 28];

Sbe23=[28 29 30];

Sbe24=[30 31 32];

Sbe25=[32];

SBE=[Sbe7,Sbe8,Sbe9,Sbe10,Sbe12,Sbe13,Sbe14,Sbe15,Sbe17,Sbe18,Sbe19,Sbe20,Sbe22,Sbe23,Sbe24 Sbe25];%影响元素集,整体编码

 

Sbn7=[2 3 6 7 8 11 12];

Sbn8=[3 4 7 8 9 12 13];

Sbn9=[4 5 8 9 10 13 14];

Sbn10=[4 9 10 14 15];

Sbn12=[7 8 11 12 13 16 17];

Sbn13=[8 9 12 13 14 17 18];

Sbn14=[9 10 13 14 15 18 19];

Sbn15=[10 14 15 19 20];

Sbn17=[12 13 16 17 18 21 22];

Sbn18=[13 14 17 18 19 22 23];

Sbn19=[14 15 18 19 20 23 24];

Sbn20=[15 19 20 24 25];

Sbn22=[17 18 21 22 23];

Sbn23=[18 19 22 23 24];

Sbn24=[19 20 23 24 25];

Sbn25=[20 24 25];

SBN=[Sbn7,Sbn8,Sbn9,Sbn10,Sbn12,Sbn13,Sbn14,Sbn15,Sbn17,Sbn18,Sbn19,Sbn20,Sbn22,Sbn23,Sbn24,Sbn25];%影响节点集,整体编码

x=zeros(1,25);

y=zeros(1,25);

for i=1:25                             %对节点进行坐标编码

    if mod(i,(N+1))==0

        x(i)=1;

    else

x(i)=(mod(i,(N+1))-1)*h;

    end

    if mod(i,(N+1))==0

        y(i)=(fix(i/(N+1))-1)*h;

    else

    y(i)=fix(i/(N+1))*h;

    end

end

a=zeros(25,25);%生成零矩阵

for i=1:16

    for j=1:m(i)

    a(NF(i),SBN(sum(m(1:(i-1)))+j))=intsbe(NF(i),SBN(sum(m(1:(i-1)))+j),np,x,y,h);%初步求取刚度矩阵

    end

end

 

 

a(1:6,:)=[];%截取矩阵,剔除无关节点

a(5,:)=[];

a(9,:)=[];

a(13,:)=[];

 

f0=-int('x^2',0,1)-int('2*x',1,0);%边界积分的数值

f=zeros(16,1);%生成零矩阵

for i=1:16   %16个未知节点F的值

f(i)=f0+((x(NF(i))^2)+(y(NF(i))^2)*x(NF(i)))*h*h*l(i)/6;

end

f

u=a\f;

u          %输出节点的值

a(:,1:6)=[];%截取矩阵,剔除已知节点

a(:,5)=[];

a(:,9)=[];

a(:,13)=[];

a     %输出刚度矩阵

U=[u(1),u(2),u(3),u(4),u(5);u(6),u(7),u(8),u(9),u(10);u(11),u(12),u(13),u(14),u(15);u(16),u(17),u(18),u(19),u(20);u(21),u(22),u(23),u(24),u(25)];%生成节点矩阵

x=0:0.25:1; y=0:0.25:1;

[X,Y]=meshgrid(x,y);

surf(X, Y, U) % 画出立体网状图

 

积分求取函数:

function Val_int=intsbe(XX,YY,np,x,y,h)

Val_int=0;

syms s t r1 r2 xx yy

for p=1:32                                         %坐标变换

if and(find(np(:,p:p)==XX),find(np(:,p:p)==YY))

    Sbe=p;

    ii=np(1,Sbe);jj=np(2,Sbe);kk=np(3,Sbe);

C1=[s,t,1;x(jj),y(jj),1;x(kk),y(kk),1];

C2=[x(ii),y(ii),1;s,t,1;x(kk),y(kk),1];

r1=det(C1)/(h*h);

r2=det(C2)/(h*h);

J=[x(ii)-x(kk),x(jj)-x(kk);y(ii)-y(kk),y(jj)-y(kk)];

    switch find(np(:,p:p)==XX)

        case 1

            Rx=r1;

         switch find(np(:,p:p)==YY)

        case 1

            Ry=r1;

            Val_int=Val_int+det(J)*quad2d(@(s,t)2*s.*s,0,1,0,@(s)1-s);

        case 2

            Ry=r2;

            Val_int=Val_int+det(J)*quad2d(@(s,t)2*s.*t,0,1,0,@(s)1-s);

        case 3

            Ry=1-r1-r2;

            Val_int=Val_int+det(J)*quad2d(@(s,t)2*s.*(1-s-t),0,1,0,@(s)1-s);

         end

        case 2

            Rx=r2;

           switch find(np(:,p:p)==YY)

        case 1

             Ry=r1;

            Val_int=Val_int+det(J)*quad2d(@(s,t)2*t.*s,0,1,0,@(s)1-s);

        case 2

            Ry=r2;

            Val_int=Val_int+det(J)*quad2d(@(s,t)2*t.*t,0,1,0,@(s)1-s);

        case 3

            Ry=1-r1-r2;

            Val_int=Val_int+det(J)*quad2d(@(s,t)2*t.*(1-s-t),0,1,0,@(s)1-s);

           end

        case 3

            Rx=1-r1-r2;

            switch find(np(:,p:p)==YY)

        case 1

             Ry=r1;

            Val_int=Val_int+det(J)*quad2d(@(s,t)2*(1-s-t).*s,0,1,0,@(s)1-s);

        case 2

            Ry=r2;

            Val_int=Val_int+det(J)*quad2d(@(s,t)2*(1-s-t).*t,0,1,0,@(s)1-s);

        case 3

            Ry=1-r1-r2;

           Val_int=Val_int+det(J)*quad2d(@(s,t)2*(1-s-t).*(1-s-t),0,1,0,@(s)1-s);

            end

    end

    

       Val_diff=diff(Rx,s)*diff(Ry,s)+diff(Rx,t)*diff(Ry,t);

       Val_int=Val_int+det(J)*Val_diff/2;

end

end

程序下载地址
  评论这张
 
阅读(977)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017