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 9 point stencil. This function employs the Hockney method without 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. Input arguments: f, function handle holding the source function g, function handle prescribing boundary condition a,b prescribing the domain N, integer value prescribing the number of interior grid points in each direction. Output arguments: soln, array holding the solution values at each of the interior grid points timer, time elapsed for comparison with other methods
0001 function [ soln, timer ] = poisson_9point_hock( 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 9 point 0006 % stencil. This function employs the Hockney method without 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 % Input arguments: 0010 % f, function handle holding the source function 0011 % g, function handle prescribing boundary condition 0012 % a,b prescribing the domain 0013 % N, integer value prescribing the number of interior grid points in 0014 % each direction. 0015 % Output arguments: 0016 % soln, array holding the solution values at each of the interior 0017 % grid points 0018 % timer, time elapsed for comparison with other methods 0019 0020 % first check user inputs 0021 if isa(f,'function_handle')==0; 0022 error('f must be a function handle'); 0023 elseif isa(g,'function_handle')==0; 0024 error('g must be a function handle'); 0025 elseif a>=b; 0026 error('must have a<b'); 0027 elseif mod(N,1)~=0; 0028 error('N must be an integer value'); 0029 end 0030 0031 h=(b-a)/(N+1); % store the value of h for hygiene 0032 0033 B=zeros(N);% initialize array to hold the source function and boundary data 0034 % fill B with the data from the source function on the grid 0035 for j=1:N; 0036 for i=1:N; 0037 B(i,j)=feval(f,a+i*h,a+j*h); 0038 end 0039 end 0040 % let y be a natural ordering of B 0041 y=B(:,1); 0042 for k=2:N; 0043 y=vertcat(y,B(:,k)); 0044 end 0045 0046 y=y*(h^2); % scale y 0047 0048 % incorporate boundary conditions at the left and right sides of the grid 0049 for i=1:N; 0050 y(i)=y(i)-(1/6)*feval(g,a,a+(i-1)*h)-(2/3)*feval(g,a,a+i*h)-(1/6)*feval(g,a,a+(i+1)*h); 0051 y((N-1)*N+i)=y((N-1)*N+i)-(1/6)*feval(g,b,a+(i-1)*h)-(2/3)*feval(g,b,a+i*h)-(1/6)*feval(g,b,a+(i+1)*h); 0052 end 0053 0054 % incorporate boundary conditions on top and bottom sides of the grid 0055 for j=1:N; 0056 y((j-1)*N+1)=y((j-1)*N+1)-(1/6)*feval(g,a+(j-1)*h,a)-(2/3)*feval(g,a+j*h,a)-(1/6)*feval(g,a+(j+1)*h,a); 0057 y(j*N)=y(j*N)-(1/6)*feval(g,a+(j-1)*h,b)-(2/3)*feval(g,a+j*h,b)-(1/6)*feval(g,a+(j+1)*h,b); 0058 end 0059 0060 % this point we've overcompensated on the corner elements, so fix this 0061 y(1)=y(1)+(1/6)*feval(g,a,a); 0062 y(N)=y(N)+(1/6)*feval(g,a,b); 0063 y(N*(N-1)+1)=y(N*(N-1)+1)+(1/6)*feval(g,b,a); 0064 y(N*N)=y(N*N)+(1/6)*feval(g,b,b); 0065 0066 tic; % start stopwatch 0067 0068 u=hockney(-10/3,2/3,2/3,1/6,y); % solve block TST system by Hockney method 0069 0070 soln(:,1)=u(1:N); % initialize first column of output 0071 0072 timer=toc; % stop stopwatch 0073 0074 % reverse the natural ordering so soln is an array representing the N^2 0075 % interior grid points. 0076 for k=1:N-1; 0077 soln=horzcat(soln,u(k*N+1:(k+1)*N)); 0078 end 0079 0080 end 0081 0082