Home > K25104 > PDEs > poisson_5point.m

poisson_5point

PURPOSE ^

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

SYNOPSIS ^

function [ soln, timer ] = poisson_5point( 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 uses MATLAB library 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. Note that MATLAB puts a 
 limit on the size of an array as a percentage of RAM, this method forms 
 an (N^2)x(N^2) array, so if N is chosen too large it may return an error.
 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( 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 uses MATLAB library to solve the linear system
0007 % arising from the application of the stencil, and incorporates a timer so
0008 % its speed can be compared with other solvers. Note that MATLAB puts a
0009 % limit on the size of an array as a percentage of RAM, this method forms
0010 % an (N^2)x(N^2) array, so if N is chosen too large it may return an error.
0011 % Input arguments:
0012 %   f, function handle holding the source function
0013 %   g, function handle prescribing boundary condition
0014 %   a,b prescribing the domain
0015 %   N, integer value prescribing the number of interior grid points in
0016 %     each direction.
0017 % Output arguments:
0018 %   soln, array holding the solution values at each of the interior
0019 %     grid points
0020 %   timer, time elapsed for comparison with other methods
0021 
0022 % first check user inputs
0023 if isa(f,'function_handle')==0;
0024     error('f must be a function handle');
0025 elseif isa(g,'function_handle')==0;
0026     error('g must be a function handle');
0027 elseif a>=b;
0028     error('must have a<b');
0029 elseif mod(N,1)~=0;
0030     error('N must be an integer value');
0031 end
0032 
0033 h=(b-a)/(N+1);  % store the value of h for hygiene
0034 
0035 % initialize array to hold the source function and boundary data
0036 B=zeros(N);
0037 
0038 % fill B with the data from the source function on the grid
0039 for j=1:N;
0040     for i=1:N;
0041         B(i,j)=feval(f,a+i*h,a+j*h);
0042     end
0043 end
0044 % let y be a natural ordering of B
0045 y=B(:,1);
0046 for k=2:N;
0047     y=vertcat(y,B(:,k));
0048 end
0049 
0050 y=y*(h^2);  % scale y
0051 
0052 % incorporate boundary conditions at the left and right sides of the grid
0053 for i=1:N;
0054     y(i)=y(i)-feval(g,a,a+i*h);
0055     y((N-1)*N+i)=y((N-1)*N+i)-feval(g,b,a+i*h);
0056 end
0057 
0058 % incorporate boundary conditions on top and bottom sides of the grid
0059 for j=1:N;
0060     y((j-1)*N+1)=y((j-1)*N+1)-feval(g,a+j*h,a);
0061     y(j*N)=y(j*N)-feval(g,a+j*h,b);
0062 end
0063 
0064 tic;    % start stopwatch
0065 % we need to construct the block tridiagonal matrix which arises by
0066 % applying the 5 point formula with natural ordering
0067 
0068 % the fastest way to do this is to create a TST matrix first
0069 A=-4*diag(ones(N*N,1))+diag(ones(N*N-1,1),1)+diag(ones(N*N-1,1),-1)+...
0070     diag(ones(N*(N-1),1),N)+diag(ones(N*(N-1),1),-N);
0071 
0072 % then remove the unwanted entries
0073 for k=1:N-1;
0074     A(N*k,N*k+1)=0;
0075     A(N*k+1,N*k)=0;
0076 end
0077 
0078 u=A\y;      % calculate u
0079 
0080 timer=toc;      % stop stopwatch
0081 
0082 soln(:,1)=u(1:N);    % initialize first column of output
0083 
0084 % reverse the natural ordering so soln is an array representing the N^2
0085 % interior grid points.
0086 for k=1:N-1;
0087     soln=horzcat(soln,u(k*N+1:(k+1)*N));
0088 end
0089 
0090 end
0091

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