function [ A,f ] = nma_GENP2D(n,force) % generate A and f for the \$Au=f\$, to solve for u on 2D based on 5 points laplacian. %--------------------------------------------- % filename: nma_GENP2D.m % function which generate A system matrix, and vector f for the %problem Au=f, to solve for u on 2D. A is Poisson discrete A matrix %based on use of standard 5 points laplacian. % % INPUT: % n: grid size. This is the number of internal grid points on one edge of % the square. Be carful, internal points only. Does not include the % boundary. So if there are total of 5 points on the edge, call this % this with n=3. i.e. n is the number of unknowns along one dimension. % force: function handler to be called to evaluate the RHS. Must % accept 2 arguments and return result. vectored version. See example % below. % % OUTPUT: % A: The system matrix, sparse matrix. use full() to convert to dense % if needed. Size is n^2xn^2 % f: column vector size is n % % note, ordering used is LEX ordering. i.e. f(1) is force on point (1,1) % which is located on the lower left edge of the square, not on the top % left corner. % % EXAMPLE % % myforce = @(X,Y) -exp( -(X-0.25).^2 - (Y-0.6).^2 ); % [A,f]=nma_GENP2(3,myforce); % full(A) % % -4 1 0 1 0 0 0 0 0 % 1 -4 1 0 1 0 0 0 0 % 0 1 -4 0 0 1 0 0 0 % 1 0 0 -4 1 0 1 0 0 % 0 1 0 1 -4 1 0 1 0 % 0 0 1 0 1 -4 0 0 1 % 0 0 0 1 0 0 -4 1 0 % 0 0 0 0 1 0 1 -4 1 % 0 0 0 0 0 1 0 1 -4 % f % % -0.8005 % -0.9301 % -0.6554 % -0.8005 % -0.9301 % -0.6554 % -0.4855 % -0.5641 % -0.3975 % Solve % u=A\f % 0.5412 % 0.7063 % 0.4904 % 0.6579 % 0.8637 % 0.5998 % 0.4261 % 0.5609 % 0.3895 % z=reshape(ans,3,3) % 0.5412 0.6579 0.4261 % 0.7063 0.8637 0.5609 % 0.4904 0.5998 0.3895 % mesh(z) A = lap2d(n,n); h = 1/(n+1); %mash spacing [X,Y] = meshgrid(h:h:1-h,1-h:-h:h); f = force(X,Y); f = flipud(f); f = reshape(f',n^2,1); %need normal coordinates end %---------------------------------------------------------- function L2 = lap2d(nx,ny) Lx=lap1d(nx); Ly=lap1d(ny); Ix=speye(nx); Iy=speye(ny); L2=kron(Iy,Lx)+kron(Ly,Ix); end %------------------------------------- % % % function L=lap1d(n) e=ones(n,1); L=spdiags([e -2*e e],[-1 0 1],n,n); end