% BVP_2D_Poisson_FD_cg.m % This MATLAB routine solves a 2-D boundary value % problem on a rectangular domain using finite differences. % Zero Dirichlet boundary conditions are used. % This program uses the conjugate gradient method to solve % the linear system, and does not even store in memory % the matrix of the problem. % K. Beers. MIT ChE. 10/17/03 function iflag_main = ... BVP_2D_Poisson_FD_cg(fun_name,L,H,num_pts); iflag_main = 0; % set number of points in grid Nx = num_pts; Ny = num_pts; Ntot = Nx*Ny; % First, place a 2-D computational grid. x = linspace(0,L,Nx); y = linspace(0,H,Ny); % Form a 2-D regular mesh such that % XG(i,j) = x(i), YG(i,j) = y(j) [XG,YG] = form_2D_mesh(x,y); % compute value of f(x,y) at each grid point FG = feval(fun_name,XG,YG,L,H); % set RHS vector of linear system b = zeros(Ntot,1); % b only has non-zero values for interior points for i=2:(Nx-1) for j=2:(Ny-1) n = get_label(i,j,Nx,Ny); b(n) = FG(i,j); end end % Now, solve the program using the conjugate gradient method, % where we suppy a routine Poisson_2D_calc_Av() that % multiplies an input vector v by the matrix A. We supply % the name of this routine in lieu of storing the matrix % in memory. CG_iter_max = round(1.1*Ntot); tol = 1e-6; phi = pcg(@Poisson_2D_calc_Av,b,tol,CG_iter_max,[],[],zeros(Ntot,1),x,y); % We now extract the results of the calculation % into the 2-D grid format and make a filled % contour plot. PHIG = zeros(size(XG)); for i=1:Nx for j=1:Ny n = get_label(i,j,Nx,Ny); PHIG(i,j) = phi(n); end end figure; contourf(XG,YG,PHIG,min(Nx,Ny)); colorbar; xlabel('x'); ylabel('y'); title('Poisson BVP (FD): \phi(x,y)'); % save the results to a .mat file save BVP_2D_Poisson_FD_cg.mat; iflag_main = 1; return; %-------------------------------------- % This routine returns the product of the A matrix with an % input vector v. It is used in the CG solver in lieu of % storing the matrix A explicitly. function Av = Poisson_2D_calc_Av(v,x,y); Av = zeros(size(v)); Nx =length(x); dx = x(2) - x(1); Ny = length(y); dy = y(2) - y(1); % We next specify equations for each boundary point. % BC # 1 i = 1; for j=1:Ny n = get_label(i,j,Nx,Ny); Av(n) = Av(n) + v(n); end % BC # 2 i = Nx; for j=1:Ny n = get_label(i,j,Nx,Ny); Av(n) = Av(n) + v(n); end % BC # 3 j = 1; for i=2:(Nx-1) n = get_label(i,j,Nx,Ny); Av(n) = Av(n) + v(n); end % BC # 4 j = Nx; for i=2:(Nx-1) n = get_label(i,j,Nx,Ny); Av(n) = Av(n) + v(n); end % We now set the linear equations for the interior points. factor_x = 1/(dx^2); factor_y = 1/(dy^2); factor_cent = 2*(factor_x + factor_y); for i=2:(Nx-1) for j=2:(Ny-1) n = get_label(i,j,Nx,Ny); Av(n) = Av(n) - factor_x*v(n-Ny) - factor_x*v(n+Ny); Av(n) = Av(n) - factor_y*v(n-1) - factor_y*v(n+1); Av(n) = Av(n) + factor_cent*v(n); end end return; %-------------------------------------- function [XG,YG] = form_2D_mesh(x,y); Nx = length(x); Ny = length(y); XG = zeros(Nx,Ny); YG = zeros(Nx,Ny); for k=1:Nx XG(k,:) = x(k)*ones(1,Ny); end for k=1:Ny YG(:,k) = y(k)*ones(Nx,1); end return; %-------------------------------------- function n = get_label(i,j,Nx,Ny); n = (i-1).*Ny + j; return;