function A=nma_laplaceRectDirchletBendCorner(nx,ny,bot,right,top,left,alpha1,beta1,alpha2,beta2,lambda,tol) % solves laplace PDE for rectangular region Dirclet BC %function A=nma_laplaceRectDirchletBendCorner(nx,ny,bot,right,top,left,alpha1,beta1,alpha2,beta2,lambda,tol) % % Function that solves the laplace PDE for rectangular region for % Dirclet boundary conditions with the lower left corner bend per % parameters alpha and beta. % %INPUT: % nx: number of grid points in the x-direction % ny: number of grid points in the y-direction % bot: bottom edge boundary value % right: right edge boundary value % top: top edge boundary value % left: left edge boundary value % alpha1: bend parameter. see problem notes for more info % alpha2: bend parameter. see problem notes for more info % beta1: bend parameter. see problem notes for more info % beta2: bend parameter. see problem notes for more info % lambda: the value lambda to use for the relaxation method of Liebmann % tol: the absolute tolerance to check for convergance. % %OUTPUT: % the solution matrix. % % %Example run: % nx=3;ny=3;bot=0;right=50;top=100;left=75; % alpha1=0.732; beta1=0.732; alpha2=1; beta2=1; lambda=1.5; tol=1; % A=nma_laplaceRectDirchletBendCorner(nx,ny,bot,right,top,left,alpha1,beta1,alpha2,beta2,lambda,tol) % %Author: Nasser Abbasi %May 25, 2003 TRUE = 1; FALSE = 0; A = zeros(nx+2,ny+2); A_old = A; % apply dirchlet conditions A(1,:) = top; A(:,1) = left; A(end,:) = bot; A(:,end) = right; %patch the 4 corner, just for cosmotics purposes A(1,1) = 0; A(end,1) = 0; A(1,end) = 0; A(end,end) = 0; A_old = A; converged = FALSE; nIter = 0; while(~converged) nIter = nIter+1; %fprintf('Iteration number %d\n',nIter); for(i=2+nx-1:-1:2) for(j=2:1:2+ny-1) if(i==2+nx-1 & j==2) % edge node, handle special case. Term1 = A(i+1,j)/( (1+beta1/beta2)*(1+ (beta1*beta2)/(alpha1*alpha2) ) ); Term2 = A(i-1,j)/( (1+beta2/beta1)*(1+ (beta1*beta2)/(alpha1*alpha2) ) ); Term3 = A(i,j+1)/( (1+alpha2/alpha1)*(1+ (alpha1*alpha2)/(beta1*beta2) ) ); Term4 = A(i,j-1)/( (1+alpha1/alpha2)*(1+ (alpha1*alpha2)/(beta1*beta2) ) ); A(i,j) = Term1+Term2+Term3+Term4; else A(i,j) = A(i,j+1) + A(i,j-1) + A(i-1,j) + A(i+1,j); A(i,j) = A(i,j)/4; end A(i,j) = lambda*A(i,j) + (1-lambda)*A_old(i,j); end end [max_row , max_row_index] = max(A(2:end-1,2:end-1)); [max_element , max_col_index] = max(max_row); pos = [max_row_index(1)+1,max_col_index+1]; epsilonA = abs(A(pos(1),pos(2)) - A_old(pos(1),pos(2)))/A(pos(1),pos(2))*100; if( epsilonA < tol ) converged = TRUE; end A_old = A; %A %epsilonA end %A=A(2:end-1,2:end-1);