fin06

 Program 8.5 Finite difference solver for 2D Poisson equation 
   with Dirichlet boundary conditions on a rectangle 
 Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps 
 Output: matrix w holding solution values 
 Example usage: w=poisson(0,1,1,2,4,4)

0001 % Program 8.5 Finite difference solver for 2D Poisson equation
0002 %   with Dirichlet boundary conditions on a rectangle
0003 % Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps
0004 % Output: matrix w holding solution values
0005 % Example usage: w=poisson(0,1,1,2,4,4)
0006 function [w,T]=fin06(xl,xr,yb,yt,z,M,N,P,K,H) 
0007 m=M+1; n=N+1; mn=m*n; 
0008 h=(xr-xl)/M; h2=h^2;
0009 k=(yt-yb)/N; k2=k^2; 
0010 x=xl+(0:M)*h; % set mesh values
0011 y=yb+(0:N)*k; 
0012 A=zeros(mn,mn);
0013 b=zeros(mn,1);
0014 L = 2; % length of input [ cm ]c
0015 %P = 5; % total power [ W ]
0016 % H = 0.005; % convective heat transfer coefficient [ W / (cm^2 * degC) ]
0017 %K = 1.68; % thermal conductivity [ W / (cm^2 * degC) ]
0018 delta = z; % depth in z-direction [ cm ]
0019 for i=2:m-1, % interior points
0020   for j=2:n-1, 
0021     A(i+(j-1)*m,i-1+(j-1)*m)=1/h2;      % u_i-1,j
0022     A(i+(j-1)*m,i+1+(j-1)*m)=1/h2;      % u_i+1,j
0023     A(i+(j-1)*m,i+(j-1)*m)=-2/h2-2/k2-((2*H)/(K*delta));  % u_i,j
0024     A(i+(j-1)*m,i+(j-2)*m)=1/k2;        % u_i,j-1
0025     A(i+(j-1)*m,i+j*m)=1/k2;            % u_i,j+1
0026     b(i+(j-1)*m)=0;
0027   end 
0028 end 
0029 for i=1:m % bottom and top boundary points
0030   j=1;
0031   A(i+(j-1)*m,i) = -3/(2*k)+(H/K);  % u_i,1
0032   A(i+(j-1)*m,i+j*m) = 2/k;       % u_i,2
0033   A(i+(j-1)*m,i+(j+1)*m) = -1/(2*k);    % u_i,3
0034   b(i+(j-1)*m)=0;
0035   j=n;
0036   A(i+(j-1)*m,i+(j-1)*m) = -3/(2*k)+(H/K);  % u_i,n
0037   A(i+(j-1)*m,i+(j-2)*m) = 2/k;       % u_i,n-1
0038   A(i+(j-1)*m,i+(j-3)*m) = -1/(2*k);    % u_i,n-2
0039   b(i+(j-1)*m)=0;
0040 end 
0041 for j=2:n-1 % left and right boundary points
0042   i=1;
0043   A(i+(j-1)*m,i+(j-1)*m) = -3/(2*h)+(H/K);  % u_1,j
0044   A(i+(j-1)*m,i+1+(j-1)*m) = 2/h;       % u_2,j
0045   A(i+(j-1)*m,i+2+(j-1)*m) = -1/(2*h);    % u_3,j
0046   b(i+(j-1)*m)= -P/(L*K*delta); % power input on left side
0047   i=m;
0048   A(i+(j-1)*m,i+(j-1)*m) = -3/(2*h)+(H/K);  % u_m,j
0049   A(i+(j-1)*m,i-1+(j-1)*m) = 2/h;       % u_m-1,j
0050   A(i+(j-1)*m,i-2+(j-1)*m) = -1/(2*h);    % u_m-2,j
0051   b(i+(j-1)*m)=0;
0052 end 
0053 v=A\b; % solve for solution in v labeling
0054 w=reshape(v(1:mn),m,n); %translate from v to w
0055 w=w+20;
0056 T = max(max(w));
0057 %figure(1);
0058 %mesh(x,y,w');
0059 ylabel('Heat Input (W)');