Home > K25104 > PDEs > poisson_9point_hock.m

poisson_9point_hock

PURPOSE ^

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

SYNOPSIS ^

function [ soln, timer ] = poisson_9point_hock( 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 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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