波动方程小结.doc

小结 电磁波的正演模型： 初值： 采用中心离散的方法解决这个二维方程， 和分别为空间步长和时间步长，为介质电导率，为介电常数，为真空中的磁导率简记为简记为简记为,则正演模型的离散形式为 = 初值条件离散为 （2-4）其中i, j, m, n, p都是正整数，边界吸收条件离散为 这就是maxwell方程的正演的时域有限差分迭代算法。数值稳定条件是波动方程代码：dt = 59*10(-13); %时间步长dx = 0.25*10(-2); %空间步长dy = 0.25*10(-2);nx = 200; %网格剖分ny = 300; % wavefield initial value for two timeu1 = zeros(nx,ny); %u1初始化u2 = u1;u3 = u1;c = u1; c(：,:) = 2*108;%波速 sigma =0.1*10(-9);% subwave, ricker waveletf = 5000*106; %频率s,tw = ricker(f,dt);%波源方程sn = length(s);%s的长度T = 350;%迭代次数for k=1:Tfor i=2:nx-1 for j=2:ny-1 if(ksn) ss = (s(k+1); else ss = 0; end u3(i,j)=(c(i,j)*dt/dx)2*(u2(i+1,j)-2*u2(i,j)+u2(i-1,j)+. +(c(i,j)*dt/dy)2*(u2(i,j+1)-2*u2(i,j)+u2(i,j-1)+. ss*(i=50)&(j=ny/2)+. +2*u2(i,j)-u1(i,j)+. c(i,j)2/2*dt*sigma*u1(i,j)/(1+c(i,j)2/2*dt*sigma)； end end u3(1,:) =0;%边界条件为0 u3(nx,:) = 0; u3(:,1) =0; u3(:,ny) = 0; up(k,:) = 0; u1 = u2;%迭代赋值 u2 = u3; imagesc(u3); %画出u3 pause(0.01); %暂停0.01s end 矩形子波接收的信号编码：dt = 59*10(-13); %time griddx = 0.25*10(-2); %space griddy = 0.25*10(-2);% c = 1000; % wave velocitynx = 200; %grid numberny = 300; % wavefield initial value for two timeu1 = zeros(nx,ny);u2 = u1;u3 = u1;c = u1;c(1:50,:) = 2*108;c(50:70,1:100) = 2*108;c(50:70,100:200) = 3*108;c(50:70,200:300) = 2*108;c(70:200,:) = 2*108;up=zeros(T,ny); % subwave, ricker waveletf = 7000*106; %frequencys,tw = ricker(f,dt);sn = length(s);sigma =0.1*10(-9);T = 350;for k=1:T k for i=2:nx-1 for j=2:ny-1 if(ksn) ss = (s(k+1); else ss = 0; end u3(i,j)=(c(i,j)*dt/dx)2*(u2(i+1,j)-2*u2(i,j)+u2(i-1,j)+. +(c(i,j)*dt/dy)2*(u2(i,j+1)-2*u2(i,j)+u2(i,j-1)+. ss*(i=2)&(j=j)+. +2*u2(i,j)-u1(i,j)+. c(i,j)2/2*dt*sigma*u1(i,j) )/(1+c(i,j)2/2*dt*sigma); end end u3(1,:) =u3(2,:); u3(nx,:) = u3(nx-1,:); u3(:,1) =u3(:,2); u3(:,ny) =u3(:,ny-1); up(k,:) = 0; u1 = u2; u2 = u3; up(k,:)=u3(2,:); imagesc(up); pause(0.01); end 圆形子波代码 dt = 59*10(-13); %time griddx = 0.25*10(-2); %space griddy = 0.25*10(-2);% c = 1000; % wave velocitynx = 200; %grid numberny = 300;up=zeros(T,ny); % wavefield initial value for two timeu1 = zeros(nx,ny);u2 = u1;u3 = u1;c = u1;c(:,:) = 3*108; % subwave, ricker waveletf = 4000*106; %frequencys,tw = ricker(f,dt);sn = length(s);sigma =0.1*10(-8); T = 300;for k=1:T k for i=2:nx-1 for j=2:ny-1 if(ksn) ss = s(k+1); else ss = 0; end u3(i,j)=(c(i,j)*dt/dx)2*(u2(i+1,j)-2*u2(i,j)+u2(i-1,j)+. +(c(i,j)*dt/dy)2*(u2(i,j+1)-2*u2(i,j)+u2(i,j-1)+. ss*(i=2)&(j=ny/2)+. +2*u2(i,j)-u1(i,j)+. c(i,j)2/2*dt*sigma*u1(i,j) )/(1+c(i,j)2/2*dt*sigma); end end