yuanhanshu.rar_matlab 正演模拟_地震 正演_地震模拟_正演 边界_边界条件

function [a]=yuanhanshu %该函数样点较多，效果较好，但有频散现象，待解决 Nx=101; Nz=101; Nt=500;hx=0.01;hz=0.01;dt=0.002; c=1; % 模拟参数 for i=1:Nx for j=1:Nz v(i,j)=1; end end for i=1:Nx for j=1:Nz r(j,i)=v(i,j)*dt/hx; r2(j,i)=(v(i,j)*dt/hx)^2; s(j,i)=2-4*(v(i,j)*dt/hx)^2; end end % 简缩“常量”使后面的量更容易代入 u=zeros(Nz,Nx,Nt); % 分布空间，预值充零 for k=1:2 for j=1:Nz for i=1:Nx u(j,i,k)=0; end end end % 初始条件 %波动方程的差分方程及震源函数和边界吸收条件 for k=2:Nt-1 for j=1:Nz for i=2:Nx-1 if j==1 u(j,i,k+1)=(2-2*r(j,i)-r2(j,i))*u(j,i,k)+2*r(j,i)*(1+r(j,i))*u(j+1,i,k)-r2(j,i) *u(j+2,i,k)+(2*r(j,i)-1)*u(j,i,k-1)-2*r(j,i)*u(j+1,i,k-1); % 上边界吸收 elseif j==Nz u(j,i,k+1)=(2-2*r(j,i)-r2(j,i))*u(j,i,k)+2*r(j,i)*(1+r(j,i))*u(j-1,i,k)-r2(j,i) *u(j-2,i,k)+(2*r(j,i)-1)*u(j,i,k-1)-2*r(j,i)*u(j-1,i,k-1); % 下边界吸收 elseif j== 51&i==51&(k-1)*dt <= 6*pi/100 %炮点 S(51,51)，子波持续时间条件 u(j,i,k+1)=dt^2*sin(50*(k-1)*dt)*exp(-188*((k-1)*dt-3*pi/100)^2)+r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1); else u(j,i,k+1)=r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k)) +s(j,i)*u(j,i,k)-u(j,i,k-1); end end end for i=1:Nx for j=2:Nz-1 if i==1 u(j,i,k+1)=(2-2*r(j,i)-r2(j,i))*u(j,i,k)+2*r(j,i)*(1+r(j,i))*u(j,i+1,k)-r2(j,i) *u(j,i+2,k)+(2*r(j,i)-1)*u(j,i,k-1)-2*r(j,i)*u(j,i+1,k-1); % 左边界吸收 elseif i==Nx u(j,i,k+1)=(2-2*r(j,i)-r2(j,i))*u(j,i,k)+2*r(j,i)*(1+r(j,i))*u(j,i-1,k)-r2(j,i) *u(j,i-2,k)+(2*r(j,i)-1)*u(j,i,k-1)-2*r(j,i)*u(j,i-1,k-1); % 右边界吸收 elseif j== 51&i==51&(k-1)*dt <= 6*pi/100 % 同上 u(j,i,k+1)=dt^2*sin(50*(k-1)*dt)*exp(-188*((k-1)*dt-3*pi/100)^2)+r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k))+s(j,i)*u(j,i,k)-u(j,i,k-1); else u(j,i,k+1)=r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k)) +s(j,i)*u(j,i,k)-u(j,i,k-1); end if j==81&i==51 a(c)=u(j,i,k+1); c=c+1; end end end end %角点的处理 for k=1:Nt u(1,1,k)=1/2*(u(1,2,k)+u(2,1,k)); u(1,Nx,k)=1/2*(u(1,Nx-1,k)+u(2,Nx,k)); u(Nz,1,k)=1/2*(u(Nz-1,1,k)+u(Nz,2,k)); u(Nz,Nx,k)=1/2*(u(Nz,Nx-1,k)+u(Nz-1,Nx,k)); end %显示图像 save U u ; % 保存数值结果 % 波场快照图显示 for k=1:20:Nt % 表示所有时刻，也可以是等间隔时间如 k=1:10:Nt figure(k) imagesc(u(:,:,k)) ; colormap(gray);colorbar; xlabel('x'); ylabel('z’'); title('Wave Field Snapshot'); axis square end clf; shg; % M=moviein(Nt); %预设画面矩阵 % for k=1:Nt % imagesc(u(:,:,k)) % colormap(hot) % xlabel('x');ylabel('z'); % title('Wave Field Movie'); % axis([0 Nx 0 Nz]); % M(:,k)=getframe; % 捕获画面 % end %plot(a); xin_xi_ti_shi=('***************** 该程序已经运行结束********************')