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