Home > K25104 > PDEs > poisson_9point.m

poisson_9point

PURPOSE ^

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

SYNOPSIS ^

function [ soln, timer ] = poisson_9point( 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 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_9point( 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 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 B=zeros(N);% initialize array to hold the source function and boundary data
0036 % fill B with the data from the source function on the grid
0037 for j=1:N;
0038     for i=1:N;
0039         B(i,j)=feval(f,a+i*h,a+j*h);
0040     end
0041 end
0042 % let y be a natural ordering of B
0043 y=B(:,1);
0044 for k=2:N;
0045     y=vertcat(y,B(:,k));
0046 end
0047 
0048 y=y*(h^2);  % scale y
0049 
0050 % incorporate boundary conditions at the left and right sides of the grid
0051 for i=1:N;
0052     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);
0053     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);
0054 end
0055 
0056 % incorporate boundary conditions on top and bottom sides of the grid
0057 for j=1:N;
0058     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);
0059     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);
0060 end
0061 
0062 % this point we've overcompensated on the corner elements, so fix this
0063 y(1)=y(1)+(1/6)*feval(g,a,a);
0064 y(N)=y(N)+(1/6)*feval(g,a,b);
0065 y(N*(N-1)+1)=y(N*(N-1)+1)+(1/6)*feval(g,b,a);
0066 y(N*N)=y(N*N)+(1/6)*feval(g,b,b);
0067 
0068 tic;    % start stopwatch
0069 % we need to construct the block tridiagonal matrix which arises by
0070 % applying the 9 point formula with natural ordering
0071 
0072 % the fastest way to do this is to create a TST matrix first
0073 A=-(10/3)*diag(ones(N*N,1))+(2/3)*diag(ones(N*N-1,1),1)+(2/3)*diag(ones(N*N-1,1),-1)+...
0074     (2/3)*diag(ones(N*(N-1),1),N)+(1/6)*diag(ones(N*(N-1)+1,1),N-1)+(1/6)*diag(ones(N*(N-1)-1,1),N+1)+(2/3)*diag(ones(N*(N-1),1),-N)...
0075     +(1/6)*diag(ones(N*(N-1)+1,1),-N+1)+(1/6)*diag(ones(N*(N-1)-1,1),-N-1);
0076 
0077 % then remove the unwanted entries
0078 for k=1:N-1;
0079     A(N*k,N*k+1)=0;
0080     A(N*k+1,N*k)=0;
0081 end
0082 for k=1:N-2;
0083     A(N*k+1,N*(k+1))=0;
0084     A(N*k,N*(k+1)+1)=0;
0085     A(N*(k+1)+1,N*k)=0;
0086     A(N*(k+1),N*k+1)=0;
0087 end
0088 A(1,N)=0;
0089 A(N,1)=0;
0090 A(N*N,N*(N-1)+1)=0;
0091 A(N*(N-1)+1,N*N)=0;
0092 
0093 A_inv=inv(A);   % calculate the inverse of A
0094 
0095 u=zeros(N*N,1); % initialize vector to hold the solution
0096 u=A_inv*y;      % calculate u
0097 
0098 timer=toc;      % stop stopwatch
0099 
0100 soln(:,1)=u(1:N);    % initialize first column of output
0101 
0102 % reverse the natural ordering so soln is an array representing the N^2
0103 % interior grid points.
0104 for k=1:N-1;
0105     soln=horzcat(soln,u(k*N+1:(k+1)*N));
0106 end
0107 
0108 end
0109

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