电脑桌面
添加小米粒文库到电脑桌面
安装后可以在桌面快捷访问

偏微分方程的有限元法求解

偏微分方程的有限元法求解_第1页
1/11
偏微分方程的有限元法求解_第2页
2/11
偏微分方程的有限元法求解_第3页
3/11
精品文档---下载后可任意编辑精品文档---下载后可任意编辑精品文档---下载后可任意编辑精品文档---下载后可任意编辑精品文档---下载后可任意编辑精品文档---下载后可任意编辑精品文档---下载后可任意编辑精品文档---下载后可任意编辑精品文档---下载后可任意编辑精品文档---下载后可任意编辑%对于 d2u/dx2=f 的 FEM 解算器,其中 f=x*(1-x)%%边界条件 u(0)=0, u(1)=0.%精确解用以比对xx=linspace(0,1,101);%产生 0-1 之间的均分指令,101 为元素个数uex=(1/6).*xx.^3-(1/12).*xx.^4-(1/12).*xx;%对力项设置高斯点的数目NGf=2;if (NGf==2)xiGf=[-1/sqrt(3);1/sqrt(3)];%ξ1、ξ2的值aGf=[1 1];else,NGf=1;xiGf=[0.0];aGf=[2.0];精品文档---下载后可任意编辑end%单元数目Ne=5;%建立网格节点x=linspace(0,1,Ne+1);%零刚性矩阵K=zeros(Ne+1,Ne+1);b=zeros(Ne+1,1);%对所有单元循环计算刚性和残差for ii=1:Ne,kn1=ii;kn2=ii+1;x1=x(kn1);x2=x(kn2);dx=x2-x1;%每一个单元的长度dxidx=2/dx;%dξ/dxdxdxi=1/dxidx;%dx/dξdN1dxi=-1/2;%dζ1/dξdN2dxi=1/2;%dζ2/dξdN1dx=dN1dxi*dxidx;%-1/(xj-xj-1)dN2dx=dN2dxi*dxidx;%1/(xj-xj-1)K(kn1,kn1)=K(kn1,kn1)-2*dN1dx*dN1dx*dxdxi;%Rj 的第二项K(kn1,kn2)=K(kn1,kn2)-2*dN1dx*dN2dx*dxdxi;K(kn2,kn1)=K(kn2,kn1)-2*dN2dx*dN1dx*dxdxi;K(kn2,kn2)=K(kn2,kn2)-2*dN2dx*dN2dx*dxdxi;%用高斯积分估量力项的积分for nn=1:NGf%NGf=2xiG=xiGf(nn);%得到高斯点的 ξN1=0.5*(1-xiG);%求 N1 和 N2(即在 xiG 的权重/插值) 形状函数在 ξ 的值N2=0.5*(1+xiG);%ζ 值fG=xiG*(1-xiG);%对 ξ 点求 fgG1=N1*fG*dxdxi;%在节点处估量权函数在高斯点的被积函数gG2=N2*fG*dxdxi;%估量是个积分值b(kn1)=b(kn1)+aGf(nn)*gG1;% aGf 为 1b(kn2)=b(kn1)+aGf(nn)*gG2;endend%在 x=0 处设置 Dirichlet 条件kn1=1;K(kn1,:)=zeros(size(1,Ne+1));K(kn1,kn1)=1;b(kn1)=0;%在 x=1 处设置 Dirichlet 条件kn1=1;K(kn1,:)=zeros(size(1,Ne+1));K(kn1,kn1)=1;b(kn1)=0;%求解方程v=K\b;%v 为 Kx=b 的解plot(x,v,'*-');%画图并比较hold on;plot(xx,uex);hold off;xlabel('x');ylabel('u');

1、当您付费下载文档后,您只拥有了使用权限,并不意味着购买了版权,文档只能用于自身使用,不得用于其他商业用途(如 [转卖]进行直接盈利或[编辑后售卖]进行间接盈利)。
2、本站所有内容均由合作方或网友上传,本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺!文档内容仅供研究参考,付费前请自行鉴别。
3、如文档内容存在违规,或者侵犯商业秘密、侵犯著作权等,请点击“违规举报”。

碎片内容

偏微分方程的有限元法求解

确认删除?
VIP
微信客服
  • 扫码咨询
会员Q群
  • 会员专属群点击这里加入QQ群
客服邮箱
回到顶部