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
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