Home > K25104 > PDEs > poisson_5point_hock_FFT.m

poisson_5point_hock_FFT

PURPOSE ^

Solves the poisson equation (grad^2)u(x,y)=f(x,y) on the

SYNOPSIS ^

function [ soln, timer ] = poisson_5point_hock_FFT( f,g,a,b,N )

DESCRIPTION ^

 Solves the poisson equation (grad^2)u(x,y)=f(x,y) on the
 2D domain [a,b]x[a,b], with Dirichlet boundary condition u(x,y)=g(x,y) on
 the boundary of the domain.  This is achieved by overlaying an  
 equispaced grid containing N^2 interior points and applying the 5 point   
 stencil. This function employs the Hockney method with FFT to solve the 
 linear system arising from the application of the stencil, and
 incorporates a timer so its speed can be compared with other solvers.
   f input argument, function handle holding the source function
   g input argument, function handle prescribing boundary condition
   a,b input arguments prescribing the domain
   N input argument, integer value prescribing the number of interior grid
   points in each direction. Note that to employ the Fast Fourier
   Transform, (N+1) must be an integer power of 2.
   soln output argument, array holding the solution values at each of the
   grid points
   timer output argument, time elapsed for comparison with other methods

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ soln, timer ] = poisson_5point_hock_FFT( f,g,a,b,N )
0002 % Solves the poisson equation (grad^2)u(x,y)=f(x,y) on the
0003 % 2D domain [a,b]x[a,b], with Dirichlet boundary condition u(x,y)=g(x,y) on
0004 % the boundary of the domain.  This is achieved by overlaying an
0005 % equispaced grid containing N^2 interior points and applying the 5 point
0006 % stencil. This function employs the Hockney method with FFT to solve the
0007 % linear system arising from the application of the stencil, and
0008 % incorporates a timer so its speed can be compared with other solvers.
0009 %   f input argument, function handle holding the source function
0010 %   g input argument, function handle prescribing boundary condition
0011 %   a,b input arguments prescribing the domain
0012 %   N input argument, integer value prescribing the number of interior grid
0013 %   points in each direction. Note that to employ the Fast Fourier
0014 %   Transform, (N+1) must be an integer power of 2.
0015 %   soln output argument, array holding the solution values at each of the
0016 %   grid points
0017 %   timer output argument, time elapsed for comparison with other methods
0018 
0019 % first check user inputs
0020 if isa(f,'function_handle')==0;
0021     error('f must be a function handle');
0022 elseif isa(g,'function_handle')==0;
0023     error('g must be a function handle');
0024 elseif a>=b;
0025     error('must have a<b');
0026 elseif mod(log2(N+1),1)~=0;
0027     error('to employ the FFT, (N+1) must be an integer power of 2');
0028 end
0029 
0030 h=(b-a)/(N+1);  % store the value of h for hygiene
0031 
0032 B=zeros(N);% initialize array to hold the source function and boundary data
0033 % fill B with the data from the source function on the grid
0034 for j=1:N;
0035     for i=1:N;
0036         B(i,j)=feval(f,a+i*h,a+j*h);
0037     end
0038 end
0039 % let y be a natural ordering of B
0040 y=B(:,1);
0041 for k=2:N;
0042     y=vertcat(y,B(:,k));
0043 end
0044 
0045 y=y*(h^2);  % scale y
0046 
0047 % incorporate boundary conditions at the left and right sides of the grid
0048 for i=1:N;
0049     y(i)=y(i)-feval(g,a,a+i*h);
0050     y((N-1)*N+i)=y((N-1)*N+i)-feval(g,b,a+i*h);
0051 end
0052 
0053 % incorporate boundary conditions on top and bottom sides of the grid
0054 for j=1:N;
0055     y((j-1)*N+1)=y((j-1)*N+1)-feval(g,a+j*h,a);
0056     y(j*N)=y(j*N)-feval(g,a+j*h,b);
0057 end
0058 
0059 tic;    % start stopwatch
0060 
0061 u=hockney_FFT(-4,1,1,0,y);  % solve block TST system by Hockney method
0062                              % with FFT
0063 timer=toc;   % stop stopwatch
0064 
0065 soln(:,1)=u(1:N);    % initialize first column of output
0066 
0067 % reverse the natural ordering so soln is an array representing the N^2
0068 % interior grid points.
0069 for k=1:N-1;
0070     soln=horzcat(soln,u(k*N+1:(k+1)*N));
0071 end
0072 
0073 end
0074 
0075

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005