function u_val = poisson(pos,rectangle,f,bdry,solver,h); % u_val = poisson(pos,rectangle,f,bdry,h,solver) % % solves poisson equation with dirichlet boundary conditions % on rectangle [-a,a] x [-b,b] using a five-point stencil % % pos [x,y] a point of the rectangle % rectangle [a,b] % 2a length in x-direction % 2b length in y-direction % f right-hand side of -\Delta u(x) = f (constant) % bdry Dirichlet data [x0,xa,y0,yb] % x0 boundary value on {-a} x [-b,b] (constant) % xa boundary value on { a} x [-b,b] (constant) % y0 boundary value on [-a,a] x {-b} (constant) % yb boundary value on [-a,a] x { b} (constant) % solver 'Cholesky' sparse Cholesky solver % 'FFT' FFT based fast solver % h h discretization parameter % % u_val solution at pos = [x,y] % the grid n = ceil(1/h); a = rectangle(1); b = rectangle(2); nx = 2*ceil(a)*n-1; ny = 2*ceil(b)*n-1; x=1:nx; y=1:ny; dx = (2*a)/(nx+1); dy = (2*b)/(ny+1); % the right hand side r = f*ones(nx,ny); r(1,y) = r(1,y)+bdry(1)/dx^2; r(nx,y) = r(nx,y)+bdry(2)/dx^2; r(x,1) = r(x,1)+bdry(3)/dy^2; r(x,ny) = r(x,ny)+bdry(4)/dy^2; % solve it switch solver case 'Cholesky' % [Dem87,§6.3.3] Ax = spdiags(ones(nx,1)*[-1 2 -1]/dx^2,-1:1,nx,nx); Ay = spdiags(ones(ny,1)*[-1 2 -1]/dy^2,-1:1,ny,ny); A = kron(Ay,speye(nx)) + kron(speye(ny),Ax); u = A\r(:); u = reshape(u,nx,ny); case 'FFT' % [Dem87,§6.7] u = dst2(r); d = 4*(sin(x'/2*pi/(nx+1)).^2*ones(1,ny)/dx^2+ ... ones(nx,1)*sin(y/2*pi/(ny+1)).^2/dy^2); u = u./d; u = dst2(u); otherwise error('MATLAB:badopt','%s: no such solver known',solver); end % the solution u_val = u(round((pos(1)+a)/dx),round((pos(2)+b)/dy)); return % subroutines for 1D and 2D fast sine transform [Dem97,p.324] function y = dst(x) n = size(x,1); m = size(x,2); y = [zeros(1,m);x]; y = imag(fft(y,2*n+2)); y = sqrt(2/(n+1))*y(2:n+1,:); return function y = dst2(x) y = dst(dst(x)')'; return