Home > K25104 > PDEs > poisson_5point_hock.m

poisson_5point_hock

PURPOSE ^

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

SYNOPSIS ^

function [ soln, timer ] = poisson_5point_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 5 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_5point_hock( f,g,a,b,N )
0002 % Solves the poisson equation (grad^2)u(x,y)=f(x,y) on the 2D domain
0003 % [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 without FFT to solve
0007 % the 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)-feval(g,a,a+i*h);
0051     y((N-1)*N+i)=y((N-1)*N+i)-feval(g,b,a+i*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)-feval(g,a+j*h,a);
0057     y(j*N)=y(j*N)-feval(g,a+j*h,b);
0058 end
0059 
0060 tic;    % start stopwatch
0061 
0062 u=hockney(-4,1,1,0,y);  % solve block TST system by Hockney method
0063 
0064 timer=toc;   % stop stopwatch
0065 
0066 soln(:,1)=u(1:N);    % initialize first column of output
0067 
0068 % reverse the natural ordering so soln is an array representing the N^2
0069 % interior grid points.
0070 for k=1:N-1;
0071     soln=horzcat(soln,u(k*N+1:(k+1)*N));
0072 end
0073 
0074 end
0075 
0076

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