Helmholtz方程的有限元解法..doc

发布时间：2022-06-08 发布人：admin 分类：说明书资料大小：0.88M 资料格式：doc 举报版权申诉

第1页 / 共7页

baidu_31133973-9087525-Helmholtz方程的有限元解法.doc.pdf-第2页.png

第2页 / 共7页

baidu_31133973-9087525-Helmholtz方程的有限元解法.doc.pdf-第3页.png

第3页 / 共7页

baidu_31133973-9087525-Helmholtz方程的有限元解法.doc.pdf-第4页.png

第4页 / 共7页

baidu_31133973-9087525-Helmholtz方程的有限元解法.doc.pdf-第5页.png

第5页 / 共7页

baidu_31133973-9087525-Helmholtz方程的有限元解法.doc.pdf-第6页.png

第6页 / 共7页

baidu_31133973-9087525-Helmholtz方程的有限元解法.doc.pdf-第7页.png

第7页 / 共7页

文本预览

问题重述：   2 u k u  ，于 1   (0,1)*(0,1) u  ，于 1   0 { x  0,0   y 1} {0    x 1, y  1} y x y 1} 于 x 1, y  {0 { x   0, 1, y 0,0  1}  1,0   2   x 1} {0        0} {    1 u  n  其中 k=1,5,15,20 采用有限元方法，这里对[0,1]*[0,1]采用 n*n 的划分，n 为偶数。程序中的变量 can 为问题中的 k，为了方便比较，采用了画图的方式。在下面在 n=10 时，分别考察 can=k=1,5,15,20 时的情况。 can=k=1 时 can=k=5 时

can=k=10 时 can=k=15 can=k=20 时时

当划分变细时的效果以 k=20 为例进行研究 n=20 n=64 时时

% n 是行与列划分的格子数，对整个[0,1]*[0,1]有 n^2 个划分,can 为方程中的参数 k %这里我们用 1,2,3 按逆时针来表示一个三角形的各个顶点 % s 是一个 n^2*10 的关联矩阵，s(i,1)表示第 i 个三角形，s(i,2),s(i,3),s(i,4)分别表示第 i 个三角形的 1,2,3 所对应的顶点 % s(i,5),s(i,6)；s(i,7)，s(i,8)；s(i,9)，s(i,10)分别表示顶点 1,2,3 所代表的坐标 %生成关联矩阵 s %A 是总刚矩阵 %声明符号变量 x y n=9; can=5; s=zeros(2*n^2,10); h=1/n; st=1/(2*n^2); A=zeros((n+1)^2,(n+1)^2); syms x y; for k=1:1:2*n^2 s(k,1)=k; q=fix(k/(2*n)); r=mod(k,(2*n)); if (r~=0) r=r; else r=2*n;q=q-1; end if (r<=n) s(k,2)=q*(n+1)+r;

s(k,3)=q*(n+1)+r+1; s(k,4)=(q+1)*(n+1)+r+1; s(k,5)=(r-1)*h; s(k,6)=q*h; s(k,7)=r*h; s(k,8)=q*h; s(k,9)=r*h; s(k,10)=(q+1)*h; s(k,2)=q*(n+1)+r-n; s(k,3)=(q+1)*(n+1)+r-n+1; s(k,4)=(q+1)*(n+1)+r-n; s(k,5)=(r-n-1)*h; s(k,6)=q*h; s(k,7)=(r-n)*h; s(k,8)=(q+1)*h; s(k,9)=(r-n-1)*h; s(k,10)=(q+1)*h; else end end %下面生成基函数 L(i)表示第 i 个点顶点的基函数 %生成单刚矩阵 d %生成单刚矩阵并将其加入总纲矩阵 d=zeros(3,3); B=zeros((n+1)^2,1); b=zeros(3,1); %生成 A 的总刚 for k=1:1:2*n^2 L(1)=(1/(2*st))*((s(k,7)*s(k,10)-s(k,9)*s(k,8))+(s(k,8)-s(k,10))*x+(s (k,9)-s(k,7))*y); L(2)=(1/(2*st))*((s(k,9)*s(k,6)-s(k,5)*s(k,10))+(s(k,10)-s(k,6))*x+(s (k,5)-s(k,9))*y); L(3)=(1/(2*st))*((s(k,5)*s(k,8)-s(k,7)*s(k,6))+(s(k,6)-s(k,8))*x+(s(k, 7)-s(k,5))*y); for i=1:1:3 for j=i:3 d(i,j)=int(int(((((diff(L(i),x))*(diff(L(j),x)))+((diff(L(i),y))*(dif f(L(j),y))))-((can^2)*L(i)*L(j))),x,0,1),y,0,1); d(j,i)=d(i,j); end

end for i=1:1:3 for j=1:1:3 end end for i=1:1:3 end b(i)=int(int((L(i)),x,0,1),y,0,1); B(s(k,(i+1)),1)=B(s(k,(i+1)),1)+b(i); end %下面根据边界条件来求解有限元方程组 Mx=B,齐次边界条件约掉了很多项 M=zeros((n+1)^2,n^2); j=1; for i=1:1:(n^2+n) A(s(k,(i+1)),s(k,(j+1)))=A(s(k,(i+1)),s(k,(j+1)))+d(i,j); if ((mod(i,(n+1)))~=1) M(:,j)=A(:,i); j=j+1; else end continue end %preanswer 是未知点的值是（n+1）^2*(n^2)的 preanswer=M\B; %得到所有节点的值 answer=zeros((n+1)^2,1); j=1; for i=1:1:(n^2+n) if ((mod(i,(n+1)))~=1) answer(i)=preanswer(j); j=j+1; else answer(i)=0; end end %生成求解后的图 Z=zeros((n+1),(n+1)); for i=1:1:(n+1)^2 s=fix(i/(n+1))+1; r=mod(i,(n+1)); if(r==0) r=n+1; s=s-1; else

end Z(r,s)=answer(i); end [X,Y]=meshgrid(1:-h:0,0:h:1); surf(X,Y,Z)

分享到：

赞收藏

资料库

Helmholtz方程的有限元解法..doc

相关推荐

课程资源

热门标签

最新资料