shengbo.rar_PML边界_声波有限差分_声波模拟_频散_高阶差分

close all;clear;clc dt=0.0005; v=2600; time=0.005; f=30; wavelet=0.05;%雷克子波时间长度 dx=5; dy=5; dz=5; pml=20;%边界厚度 R=0.001;%反射系数 a=31+pml;b=31+pml;c=31+pml;%网格中心位置 line=61+2*pml;row=61+2*pml;vertical=61+2*pml;%网格数 p0=zeros(line,row,vertical);%网格点上的赋值 p1=zeros(line,row,vertical); px0=zeros(line,row,vertical);%p0表示前一时刻 px1=zeros(line,row,vertical);%p1表示现在时刻 py0=zeros(line,row,vertical); py1=zeros(line,row,vertical); pz0=zeros(line,row,vertical); pz1=zeros(line,row,vertical); Ax0=zeros(202,202,202);%为边界条件定义的新量 Ax1=zeros(202,202,202); Ay0=zeros(202,202,202); Ay1=zeros(202,202,202); Az0=zeros(202,202,202); Az1=zeros(202,202,202); ddx=zeros(line,row,vertical);%加上边界条件的全部网格 ddy=zeros(line,row,vertical); ddz=zeros(line,row,vertical); %pml边界条件 %八个点 for i=1:line for j=1:row for k=1:vertical if i>=1&&i<=pml&&j>=1&&j<=pml&&k>=1&&k<=pml%x，y，z上边界 x=pml-i;y=pml-j;z=pml-k; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=line-pml&&i<=line&&j>=1&&j<=pml&&k>=1&&k<=pml%x下边界，y，z上边界 x=i+pml-line;y=pml-j;z=pml-k; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=1&&i<=pml&&j>=row-pml&&j<=row&&k>=1&&k<=pml%x上边界，y下边界，z上边界 x=pml-i;y=j+pml-row;z=pml-k; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=1&&i<=pml&&j>=1&&j<=pml&&k>=vertical-pml&&k<=vertical%x，y上边界，z下边界 x=pml-i;y=pml-j;z=k+pml-vertical; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=line-pml&&i<=line&&j>=row-pml&&j<=row&&k>=1&&k<=pml%x，y下边界，z上边界 x=i+pml-line;y=pml+j-row;z=pml-k; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=line-pml&&i<=line&&j>=1&&j<=pml&&k>=vertical-pml&&k<=vertical%x下边界，y上边界，z下边界 x=pml+i-line;y=pml-j;z=pml+k-vertical; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=1&&i<=pml&&j>=row-pml&&j<=row&&k>=vertical-pml&&k<=vertical%x上边界，y，z下边界 x=pml-i;y=pml+j-row;z=pml+k-vertical; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=line-pml&&i<=line&&j>=row-pml&&j<=row&&k>=vertical-pml&&k<=vertical%x，y，z下边界 x=pml+i-line;y=pml+j-row;z=pml+k-vertical; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); % 十二条棱线 elseif i>=pml&&i<=line-pml&&j>=1&&j<=pml&&k>=1&&k<=pml%x方向，y，z=0的棱线 y=pml-j;z=pml-k; ddx(i,j,k)=0; ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=1&&i<=pml&&j>=pml&&j<=row-pml&&k>=1&&k<=pml%y方向，x，z=0的棱线 x=pml-i;z=pml-k; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=0; ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=1&&i<=pml&&j>=1&&j<=pml&&k>=pml&&k<=vertical-pml%z方向，x，y=0的棱线 x=pml-i;y=pml-j; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=0; elseif i>=line-pml&&i<=line&&j>=pml&&j<=row-pml&&k>=1&&k<=pml%y方向，x=max，z=0的棱线 x=pml+i-line;z=pml-k; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=0; ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=pml&&i<=line-pml&&j>=1&&j<=pml&&k>=vertical-pml&&k<=vertical%x方向，y=0，z=max的棱线 y=pml-j;z=pml+k-vertical; ddx(i,j,k)=0; ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=pml&&i<=line-pml&&j>=row-pml&&j<=row&&k>=1&&k<=pml%x方向，y=max，z=0的棱线 y=pml+j-row;z=pml-k; ddx(i,j,k)=0; ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=1&&i<=pml&&j>=pml&&j<=row-pml&&k>=vertical-pml&&k<=vertical%y方向，x=0，z=max的棱线 x=pml-i;z=pml+k-vertical; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=0; ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=1&&i<=pml&&j>=row-pml&&j<=row&&k>=pml&&k<=vertical-pml%z方向，x=0，z=max的棱线 x=pml-i;y=pml+j-row; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=0; elseif i>=line-pml&&i<=line&&j>=1&&j<=pml&&k>=pml&&k<=vertical-pml%z方向，x=max，y=0的棱线 x=pml+i-line;y=pml-j; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=0; elseif i>=pml&&i<=line-pml&&j>=row-pml&&j<=row&&k>=vertical-pml&&k<=vertical%x方向，y=max，z=max的棱线 y=pml+j-row;z=pml+k-vertical; ddx(i,j,k)=0; ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=line-pml&&i<=line&&j>=pml&&j<=row-pml&&k>=vertical-pml&&k<=vertical%y方向，x=max，z=max的棱线 x=pml+i-line;z=pml+k-vertical; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=0; ddz(i,j,k)=-log(R)*3*v*(z*5)^2/(2*(pml*5)^3); elseif i>=line-pml&&i<=line&&j>=row-pml&&j<=row&&k>=pml&&k<=vertical-pml%z方向，x=max，y=max的棱线 x=pml+i-line;y=pml+j-row; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=-log(R)*3*v*(y*5)^2/(2*(pml*5)^3); ddz(i,j,k)=0; %六个面 elseif i>=1&&i<=pml&&j>=pml&&j<=row-pml&&k>=pml&&k<=vertical-pml%x=0，y，z方向的面 x=pml-i; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=0; ddz(i,j,k)=0; elseif i>=line-pml&&i<=line&&j>=pml&&j<=row-pml&&k>=pml&&k<=vertical-pml%x=max，y，z方向的面 x=pml+i-line; ddx(i,j,k)=-log(R)*3*v*(x*5)^2/(2*(pml*5)^3); ddy(i,j,k)=0; ddz(i,j,k)=0; elseif i>=pml&&i<=line-pml&&j>=row-pml&&j<=row&&k>=pml&&k<=vertical-pml%y=max，x，z方向的面 y=pml+j-row; ddx(i,j,k)=0; ddy(i,j,k)=-log(