首页 | 期刊简介 | 主要栏目 | 投稿指南 | 编委会 | 广告合作及广告客户 | 期刊订阅 | 书籍推荐 | 联系我们
一个简单二阶常微分方程的有限元法求解过程(并附上有限元求解程序)

本文转自湖南大学崔向阳团队微信

本次为大家分享的是一个简单二阶常微分方程的有限元法求解过程(并附上有限元求解程序)

源程序代码:


% 一个简单二阶常微分方程的有限元法求解过程

% HNU崔向阳团队 2017/03/14


clear

clc


nel=50;

nnode=nel+1;

nnel=2;

ndof=1;

sdof=nnode*ndof;

for inode=1:nnode

    xm(1,inode)=(inode-1)*1/nel;

    xm(2,inode)=0;

end

for iel=1:nel

    nodes(iel,1)=iel;

    nodes(iel,2)=iel+1;

end


bcdof=[1 nnode];

bcval=[0 0];


ff=sparse(sdof,1);

kk=sparse(sdof,sdof);

for iel=1:nel

    for inode=1:nnel

        nd(inode)=nodes(iel,inode);

        x(inode)=xm(1,nd(inode));

        y(inode)=xm(2,nd(inode));

    end

    [k,f]=cal_k_f(x,y);

    [index]=feeldof(nd,nnel,ndof);

    kk=feasmb_k(kk,k,index);

    ff=feasmb_f(ff,f,index);

end


[kk,ff]=feaply(kk,ff,bcdof,bcval);


[LL UU]=lu(kk);

utemp=LL\ff;

disp=UU\utemp;


x0=xm(1,:);

y0=disp(:,1)';

xref=(0:0.01:1);

for inode=1:101

yref(inode)=sin(xref(inode))/sin(1)-xref(inode);

end


figure(1)

hold on

plot(x0,y0,'bs','LineWidth',2)

plot(xref,yref,'r-','LineWidth',2)


function [k,f]=cal_k_f(x,y)

x1=x(1);

x2=x(2);

y1=y(1);

y2=y(2);

L=sqrt((x2-x1)^2+(y2-y1)^2);

k=[1/L -1/L;-1/L 1/L;]-[L/3 L/6;L/6 L/3];

f=[(x2*x2+x1*x2-2*x1*x1)/6;(2*x2*x2-x1*x2-x1*x1)/6];

return

end


function [kk]=feasmb_k(kk,k,index)

edof=length(index);

for i=1:edof

    ii=index(i);

    for j=1:edof

        jj=index(j);

        kk(ii,jj)=kk(ii,jj)+k(i,j);

    end

end

return

end


function [ff]=feasmb_f(ff,f,index)

edof=length(index);

for i=1:edof

    ii=index(i);

    ff(ii)=ff(ii)+f(i);

end

return

end


function [kk,ff]=feaply(kk,ff,bcdof,bcval)

n=length(bcdof);

sdof=size(kk);

for i=1:n

    c=bcdof(i);

    for j=1:sdof

       kk(c,j)=0;

       kk(j,c)=0;

    end

    kk(c,c)=1;

    ff(c)=bcval(i);

end

return

end


function [index]=feeldof(nd,nnel,ndof)

 edof = nnel*ndof;

 k=0;

   for i=1:nnel

     start = (nd(i)-1)*ndof;

       for j=1:ndof

         k=k+1;

         index(k)=start+j;

       end

   end

return

end

 

联系我们:


联系人:崔向阳

E-mail:cui435@163.com

 

欢迎关注《计算机辅助工程》微信公众平台,微信公众号:CAEChina

发布时间:2017年03月21日

您是第5121653位访问者
地址:上海市浦东新区海港大道1550号上海海事大学图书馆B511 邮编:201306
联系电话:021-38284908 传真:021-38284916 E-mail:smucae@163.com
微信号:CAEChina   QQ:274226135   QQ群(中国CAE技术交流):204762677
本系统由北京勤云科技发展有限公司设计
沪ICP备11028865号-2