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
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