function u = fish2d(f) % FISH2D Poisson solver in 2D based on matlab fft, % square geometry, homogeneous Dirichlet boundary conditions. % % U = FISH2D(F) is the solution of LAMPF(U) = F. % % Requires uniform mesh (2^k-1)^2 interior nodes. % % C. T. Kelley, January 6, 1994 % % This code comes with no guarantee or warranty of any kind. % % Solves - u_xx - u_yy = f % % function u = fishsd(f) % % global fish_lal fish_ual ua; n=length(f); m=sqrt(n); hh=.5/(m+1); h2=(m+1)*(m+1); e=ones(m,1); % % Don't compute eigenvalues, factor tridiagonal matrices % or allocate space for the 2-D representation of the solution % more than once. Note the use of global variables to store % fish_ual, fish_lal, and ua. % if exist('fish_ual')==1 [f1, f2]=size(fish_ual); m1=sqrt(f1); else m1=m; end if (exist('fish_ual') == 0) | (m1 ~= m) % % allocate room for ua, the 2D representation of the solution % ua=zeros(m,m); % % set up diagonal of system using eigenvalues of 1D operator % x=1:m; x=pi*hh*x'; lam=sin(x); lam=h2*(2*ones(m,1) + 4*lam.*lam); en=h2*ones(n,1); fn=en; gn=en; dx=zeros(n,1); for i=1:m fn(i*m)=0; gn((i-1)*m+1)=0; dx((i-1)*m+1:i*m)=lam(i)*ones(m,1); end al=spdiags([-fn, dx, -gn], [-1,0,1], n,n); % % factor the tridiagonal matrix for solution in the x direction % [fish_lal,fish_ual]=lu(al); end % % Map f into a rectangular array. % Take sine transform in y direction. % ua(:)=f; ua=sintv(ua).'; % % Solve in the x direction % ua(:)=fish_ual\(fish_lal\ua(:)); % % Take inverse transform. % ua=isintv((ua.')); % % Map uz into a linear array. % u=ua(:); % function u=sintv(z) % sine transform, assume that size(z)=[nx,ny]=[2^k-1,ny] % % C. T. Kelley, May 1994 % % This code comes with no guarantee or warranty of any kind. % [nx,ny]=size(z); ww=-4*fft([zeros(ny,1), z']',2*nx+2); u=imag(ww(2:nx+1,:)); % function u=isintv(z) % inverse sine transform, assume that size(z)=[nx,ny]=[2^k-1,ny] % % C. T. Kelley, May 1994 % % This code comes with no guarantee or warranty of any kind. % [nx,ny]=size(z); ww=ifft([zeros(ny,1), z']',2*nx+2); u=imag(ww(2:nx+1,:));