| 查看: 628 | 回复: 1 | ||
[求助]
二维TM波，PML边界，matlab程序。但是仿真结果有问题，求高手解决
|
|
clc; clear; ft=1.5; fs=20; NX=100;%X方向网格数 NY=150;%Y方向网格数 NA=20;%吸收层厚度 bxrhx=1:100; bxirhx=2:150; byrhy=2:100; byirhy=1:150; bxrez=1:100; byrez=1:150; hx=zeros(100,151); hy=zeros(101,150); ezx=zeros(100,150); ezy=zeros(100,150); sigmamy=hx; sigmamx=hy; sigmax=ezx; sigmay=ezy; sigmamax=0.2; sigma=sigmamax*(NA:-1:1)/NA; sigmam=sigmamax*(2*NA-1:-2:1)/2/NA; sigmamx([1:NA,end:-1:end-NA+1],  =repmat([sigmam,sigmam]',1,150);sigmax([1:NA,end:-1:end-NA+1], =repmat([sigma,sigma]',1,150);sigmamy(:,[1:NA,end:-1:end-NA+1])=repmat([sigmam,sigmam],100,1); sigmay(:,[1:NA,end:-1:end-NA+1])=repmat([sigma,sigma],100,1); f1hx=exp(-sigmamy(bxrhx,bxirhx)); f1hy=exp(-sigmamx(byrhy,byirhy)); f2hx=zeros(size(hx)); f2hy=zeros(size(hy)); f2hx(sigmamy == 0)=1/ft; f2hy(sigmamx==0)=1/ft; f2hx(sigmamy~=0)=(1-exp(-sigmamy(sigmamy~=0)))./sigmamy(sigmamy~=0)/ft; f2hy(sigmamx~=0)=(1-exp(-sigmamx(sigmamx~=0)))./sigmamx(sigmamx~=0)/ft; f2hx=f2hx(bxrhx,bxirhx); f2hy=f2hy(byrhy,byirhy); f1ezx=exp(-sigmax); f1ezy=exp(-sigmay); f2ezx=zeros(size(ezx)); f2ezy=zeros(size(ezy)); f2ezx(sigmax==0)=1/ft; f2ezy(sigmay==0)=1/ft; f2ezx(sigmax~=0)=(1-exp(-sigmax(sigmax~=0)))./sigmax(sigmax~=0)/ft; f2ezy(sigmay~=0)=(1-exp(-sigmay(sigmay~=0)))./sigmay(sigmay~=0)/ft; centerx=floor((1+NX)/2); centery=floor((1+NY)/2); for nn=1:150 ezx(centerx,centery)=ezx(centerx,centery)+sin(2*pi*nn/ft/fs); ezy(centerx,centery)=ezy(centerx,centery)+sin(2*pi*nn/ft/fs); hx(bxrhx,bxirhx)=f1hx.*hx(bxrhx,bxirhx)-f2hx.*(ezx(bxrhx,bxirhx)-ezx(bxrhx,bxirhx-1)+ezy(bxrhx,bxirhx)-ezy(bxrhx,bxirhx-1)); hy(byrhy,byirhy)=f1hy.*hx(byrhy,byirhy)+f2hy.*(ezx(byrhy,byirhy)-ezx(byrhy-1,byirhy)+ezy(byrhy,byirhy)-ezy(byrhy-1,byirhy)); ezx=f1ezx.*ezx+f2ezx.*(hy(bxrez+1, -hy(bxrez,);ezy=f1ezy.*ezy-f2ezy.*(hx(:,byrez+1)-hx(:,byrez)); end 程序可以运行。但是自己画出的等位线图不对，contour（ezx+ezy）,contour(hx),contour(hy),都不对，求解决。 |
回复此楼» 猜你喜欢
今天系统多次维护，明天很可能放榜！
已经有5人回复
-
欢迎发来filecode的Mz6后的代码验证其规律
已经有74人回复
-
岩土工程学报什么时候才能终审完呐
已经有4人回复
-
时间戳又变了8-15
已经有25人回复
-
2027广东省杰青
已经有6人回复
-
快农历七夕节了，轻松一下，男人悄悄话，女施主请不要进来。
已经有3人回复
-
哪位老哥知道今年的国自然具体哪一天放榜？
已经有9人回复
-
咱们一起用铁证分析2026国家社科基金中标与否
已经有26人回复
-
我的国基提前知道中了，可是同事的操作让我实在接受不了，怎么会有这样的人
已经有16人回复
-
有时候，自然基金真的不能太认真 （我的申报经验）
已经有12人回复
2楼2015-11-12 22:12:34
25