使用谱方法和有限差分法根据位移求解应变

使用谱方法和有限差分法，根据位移求解应变

clear;close all;
Nx = 128;Ny=Nx;dx=1;dy=dx;
delkx = (2.0 *pi ) / (Nx * dx);
delky = (2.0 *pi ) / (Ny * dy);
kx = delkx * [0:Nx/2-1 -Nx/2:-1];
ky = delky * [0:Ny/2-1 -Ny/2:-1];
[kX,kY] = meshgrid(kx,ky);

%kX=kX';kY=kY';%需要转置;%不需要转置

X = 0 : 1 : Nx-1;Y = 0 : 1: Ny-1;
fangcha= 100;fangcha2 = 220;
u1 = zeros(Nx, Ny);
u2 = zeros(Nx,Ny); 
A = 10 ^ 3;
center_u1 = 50;center_u2 = 55;
for row = 1 : 1 : Nx
    for col = 1 : 1 : Ny
        u1(row, col) = (X(row) - center_u1) .* (X(row)-center_u1) + (Y(col) - center_u1) .* (Y(col) - center_u1);
        u2(row, col) = (X(row) - center_u2) .* (X(row)-center_u2) + (Y(col) - center_u2) .* (Y(col) - center_u2);
    end
end

u1 = -u1/(2*fangcha);
u1 = A * exp(u1) / (sqrt(2*pi) * sqrt(fangcha));
u2 = -u2/(2*fangcha);
u2 = A * exp(u2) / (sqrt(2*pi) * sqrt(fangcha2));

surf(X, Y, u1);surf(X, Y, u2);
%% fourier spectral method
e11 = real(ifft2(1i .* kX .* fft2(u1))); 
e22 = real(ifft2(1i .* kY .* fft2(u2)));
e12 = real(ifft2(1/2 *(1i .* kX .* fft2(u2) + 1i .* kY .* fft2(u1))));
%% Difference method 周期性边界条件
for i=1:Nx
for j=1:Ny
    jp=j+1;jm=j-1;
    ip=i+1;im=i-1;
if(im == 0)
    im=Nx;
end
if(ip == (Nx+1))
  ip=1;
end
if(jm == 0) 
  jm = Ny;
end
if(jp == (Ny+1))
  jp=1;
end
hne=u1(ip,j);hnw=u1(im,j);hns=u1(i,jm);hnn=u1(i,jp);hnc=u1(i,j);
hne1=u2(ip,j);hnw1=u2(im,j);hns1=u2(i,jm);hnn1=u2(i,jp);hnc1=u2(i,j);

e11_diff(i,j) = (hnn - hns) / 2;
e22_diff(i,j) = (hne1 - hnw1) / 2;
e12_diff(i,j) = 1/2*((hne-hnw)/2 + (hnn1-hns1)/2);
% lap_con(i,j) =(hnw + hne + hns + hnn -4.0*hnc)/(dx*dy);
end;
end;
error = e11 - e11_diff;
err = max(error(:));
subplot(3,3,1);imagesc(e11);colorbar;title('e11 spectral');
subplot(3,3,2);imagesc(e11_diff);colorbar;title('e11 diff');
subplot(3,3,3);imagesc(e22);colorbar;title('e22 spectral');
subplot(3,3,4);imagesc(e22_diff);colorbar;title('e22 diff');
subplot(3,3,5);imagesc(e12);colorbar;title('e12 spectral');
subplot(3,3,6);imagesc(e12_diff);colorbar;title('e12 diff');
subplot(3,3,7);imagesc(u1);colorbar;title('u1');
subplot(3,3,8);imagesc(u2);colorbar;title('u2');