function [X,Y,u,residNorm] = nma_P2DDIRSOR(N,WB,SB,EB,NB,C,rhs) %Solve 2D poisson PDE on unit square. Dirichlet B.C % file name: nma_P2DDIRSOR.m % % Solve the 2D poisson PDE on unit square with Dirichlet B.C % % INPUT: % N = number of unknowns in one dimension (this is the same % as internal grid points in one direction % WB = west edge boundary condition. Must be a number. % SB = south edge boundary condition. Must be a number. % EB = east edge boundary condition. Must be a number. % NB = north edge boundary condition. Must be a number. % C = a constant to multiply the h^2 with in order % to generate tolerance to use for terminating iterative % solver. h is the spacing between grid points. % typical values of C is 1 or 0.1 or 0.01. The smaller this % number, the longer it will take to converge % f = a handle to the RHS function. The function must be defined to % accept 2 inputs X,Y and return back the function result % OUTPUT: % X,Y,Z = matrices which gives the x,y coordinates of the grid, with % the solution being the z matrix, giving the value of the % solution at each coordinate. % residNorm = a vector of the residual at each iteration. The length % of the vector gives the number of iterations needed % to converge. % % EXAMPLE 1 % f = @(X,Y)-exp(-(X-0.25).^2-(Y-0.6).^2) % [x,y,z,r] = nma_2DPoissonDirichlet(15,0,0,0,0,1,f); % surf(x,y,z) % plot(1:length(r),r) % % EXAMPLE 2 % close all; clear all; % N=15; C=0.001; % f=@(X,Y) -2*((1-6*X.^2).*Y.^2.*(1-Y.^2)+(1-6*Y.^2).*X.^2.*(1-X.^2)); % [x,y,z,r] = nma_2DPoissonDirichlet(N,0,0,0,0,C,f); % surf(x,y,z) % h=1/(N+1); % [X,Y] = meshgrid(0:h:1, 0:h:1); % z2 = (X.^2-X.^4).*(Y.^4-Y.^2); % figure; surf(X,Y,z2) % % by Nasser M. Abbasi % DOPLOTS = false; %set to false if do not want to see plots h = 1/(N+1); w = 2*(1-pi*h); % evaluate f(x,y) on grid [X,Y] = meshgrid(0:h:1, 0:h:1); f = rhs(X,Y); nPoints = size(f,1); fv = reshape(f,nPoints^2,1); % use grid vector norm normf = sqrt(h) * norm(fv,2); % initialize space (grid) for residual calculation and for solution resid = zeros(nPoints,nPoints); u = zeros(nPoints,nPoints); unew = u; done = false; %flag set to true in loop below when it converges tolerance = C*h^2; % set tolerance k = 1; % initialize iteration counter r = zeros(100,1); while not(done) for i = 2 : nPoints-1 for j = 2 : nPoints-1 resid(i,j)= f(i,j) - 1/h^2 * ( u(i-1,j) + u(i+1,j) + ... u(i,j-1) + u(i,j+1) - 4*u(i,j) ); unew(i,j) = w/4* ( unew(i-1,j) + u(i+1,j) + unew(i,j-1)+... u(i,j+1) - h^2*f(i,j)) + (1-w)*u(i,j); end end u = unew; residv = reshape(resid,nPoints^2,1); % use grid vector norm r(k) = sqrt(h) * norm(residv,2); if mod(k,100)==0 tmp=r; r=zeros(length(tmp)+100,1); r(1:length(tmp))=tmp; end if ( r(k) / normf)