Labfans是一个针对大学生、工程师和科研工作者的技术社区。 | 论坛首页 | 联系我们(Contact Us) |
![]() |
![]() |
#1 |
初级会员
注册日期: 2008-05-04
年龄: 39
帖子: 7
声望力: 0 ![]() |
![]()
function [zeta,s,U]=diffraction(m,n)
dzeta=1/m; ds=6/n; ta=dzeta/(ds*ds); zeta=zeros(m+1,1); s=zeros(n-1,1); U=zeros(m+1,n-1); b=zeros(n-1,1); ll=zeros(n-1,1); ul=zeros(n-1,1); zeta=[0:dzeta:1]; s=[-(n-2)/2*ds:ds ![]() %r=3; %thet0=4.5*10^-3; f=inline('sqrt(0.1/sqrt(pi)/(4*10^-3)*exp(-(2*10^-3).^2/(4*10^-3).^2))*exp(-s.^2/(0.69*0.69))','s'); for j=1:n-1 U(1,j)=f(s(j)); end; ne=2.3; lambda0=0.5*10^-6; ke=2*pi/lambda0*ne; thet=-2*10^-3; x0=18.265*10^-6; alpha=thet*ke*x0; g0=alpha/2*ds*ta; g1=1/4*sqrt(-1)*ta; for k=1:n-2 A(k,k+1)=g0-g1; A(k+1,k)=-g1; A(k,k)=-g0+1+2*g1; end A(1,1)=A(1,1)-g1; A(n-1,n-1)=1+g1; % 分解矩阵A ul=A(1,1); for l=2:n-1; ll(l)=A(l,l-1)/ul(l-1); ul(l)=A(l,l)-A(l-1,l)*ll(l); end for k=1:m; b(1)=(1+g0-g1)*U(k,1)+(-g0+g1)*U(k,2); for j=2:n-2; b(j)=g1*U(k,j-1)+(1+g0-2*g1)*U(k,j)+(-g0+g1)*U(k,j+1); end b(n-1)=g1*U(k,n-2)+(1-g1)*U(k,n-1); %用追赶法求解线性方程组AU=b y(1)=b(1); for l=2:n-1 y(l)=b(l)-ll(l)*y(l-1); end; U(k+1,n-1)=y(n-1)/ul(n-1); for l=n-2 ![]() U(k+1,l)=(y(l)+(-g0+g1)*U(k+1,l+1))/ul(l); end; end mesh(s,zeta,abs(U).^2) xlabel('s') ylabel('\xi') zlabel('|U|.^2') |
![]() |
![]() |