function nma_P2DDIRJCB_S() %script solves 2D Poission PDE on unit square using Jacobian iterative method % file name: nma_P2DDIRJCB_S.m % % Small matlab script which solves the 2D Poission PDE % on unit square using the Jacobian iterative method % with stopping criteria that of the relative residual. % % This is a standalone script which require no additional % toolboxes to run. % % There are 3 parameters to change only below. Edit these and rerun % to see the solution % % n: THis is overall number of grid points on one edge. % c: Tolerance constant. Use 1 or 0.1 or 0.01, used for checking for % convergence. % force: this is the RHS of the PDE, edit to change the function % % by Nasser M. Abbasi 11/17/2010 % %---------- CONFIGURATION PARAMETERS --------------------------- n = 9; % all grid points, including internal points, on one edge c = 0.01; % constant of tolerance myforce = @(X,Y) -exp( -(X-0.25).^2 - (Y-0.6).^2 ); % RHS of PDE %---------- END CONFIGURATION PARAMETERS ------------------------ h = 1/(n-1); % grid spacing tol = c * h^2; % tolerance res = zeros(n,n); % storage for residual u = res; % storage for solution A = [0 1 0;1 0 1;0 1 0]; % 5 point stencil for u calculation B = [0 1 0;1 -4 1;0 1 0]; % 5 point stencil for residue calculation k = 0; % loop counter [X,Y] = meshgrid(0:h:1,0:h:1); % coordinates f = myforce(X,Y); normf = norm( reshape(f(2:end-1,2:end-1),1,(n-2)^2),2); % find norm % Now start the iteration solver, stop when relative residual < tolerance figure; done = false; i = 2:n-1; j = 2:n-1; %the iterations vectorized while ~done k = k+1; u(i,j) = (1/4)*( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - h^2 * f(i,j) ); res(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) ); mesh(X,Y,u); hold on; title(sprintf('solution at iteration %d',k)); %zlim([0,.07]); drawnow; if abs(normf)