load "dfft" load "medit" func complex[int,int] matrixDFFT(real epsi,int n) { // Matirz de Vandermonde // - Si epsi=-1, obtenemos la transformada // - Si epsi=1, obtenemos la inversa // n=NN/2; // Se construye desde fuera hacia dentro, por filas real N=2*n+1; complex[int,int] A(N,N); A=0; complex WN=exp(epsi*2.*pi*1i/N); for(int ii=-n;ii<=0;ii++) { for(int jj=-n;jj<=0;jj++) { A(ii+n,jj+n)=WN^(ii*jj); A(N-ii-n-1,N-jj-n-1)=WN^(ii*jj); A(ii+n,N-jj-n-1)=WN^(-ii*jj); A(N-ii-n-1,jj+n)=WN^(-ii*jj); } } if(epsi<0)A/=N; return A; } func complex[int,int] vector2matrix(complex[int] & u,int n) { complex[int,int] matrixu(n,n); for(int ii=0;ii ZZ,zz; // Creamos el espectro artificialmente real module; if (false) { for(int k=-NN/2;k<=NN/2;k++) { for(int l=-NN/2;l<=NN/2;l++) { module=sqrt(k^2+l^2); if(sqrt(module)<1-1.e-3) { ZZ[][(k+NN/2)+(NN+1)*(l+NN/2)]=0.; } else { ZZ[][(k+NN/2)+(NN+1)*(l+NN/2)]=module^(-(alpha+1)/2); } } }} ZZ[][NN/2*(2+NN)]=1.0; VhP1 zr; zr=real(ZZ); medit("ZZ",Th,zr); // Calculamos la transformada inversa zz[] = dfft(ZZ[],NN+1,1); zr = real(zz); plot(zr,wait=1,fill=1); medit("zr",Th,zr); //plot(zr,zr,wait=1); // solution in the corners of the domain with instability zz[]=DFFT2d(ZZ[],1,NN/2); zr = real(zz); plot(zr,wait=1,fill=1); medit("zrd",Th,zr); //plot(zr,zr,wait=1); // solution in the center of the domain with no instability