How to implement a mixed boundary conditions into 2D steady state heat conduction equation?

Hi

I have 2D steady heat conduction equation on the unit square subject to the following mixed Dirichlet/Neumann boundary conditions.

(x,0) =5 T(0,y)=0

T(x,1)=sin(x) T(1,y)=1-y

Use uniform grid in x and y of 21 points in each direction. Apply an initial condition of T(x,y)=0 . Iterate until the maximum change is less than 0.1. And the tolerance value 1.0e-4

here is my code but I'm having difficulties implement the Neumann boundary condition.

for I=1:2
    tic
nx = 21;   %step size x -direction
ny = 21;   %step size y -direction 
Lx = 1; 
Ly = 1;
x = linspace(0,Lx,nx);
y = linspace(0,Ly,nx);
dx = x(2)-x(1);
dy = y(2)-y(1);
k = 0;
% Intial condition 
T = zeros(nx,ny);
% BC 
T(1:nx,nx) = (sin(pi*x(:,1)));
T(1:nx,1) = 5;
T(ny,1:ny) = 0;
T(1,1:ny) = (1-y(1,:));
Told = T; 
tol = 1.0e-4;
error = 1;
k = k+1;
% iterative solver 
 iterative_solver = input('iterative_solver(1. jacobi 2. gauss) = ');
 % Jacobi 
 if iterative_solver ==1
     jacobi_iter = 1;
     while(error>tol)
         for i = 2:nx-1
             for j = 2:ny-1
                  T(i,j) = 0.25*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1) + Told(i,j+1));
             end 
         end 
         error = max(max(abs(Told-T)));
         jacobi_iter = jacobi_iter+1;
         Told=T;
     end 
     contourf(x,y,T,'ShowText','on')
     colorbar;
      title_text=sprintf('total jacobi iteration = %d @steady state',jacobi_iter);
 title(title_text)
 xlabel('x-axis');
 ylabel('y-axis');
 end 
 % Gauss 
 if iterative_solver ==2
 gs_iter =1;
 while (error>tol)
     for i = 2:nx-1
         for j = 2:ny-1
             T(i,j) = (0.25)*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1));
         end 
     end 
     error = max(max(abs(Told-T)));
     gs_iter=gs_iter+1;
     Told = T;
 end 
 contourf(x,y,T,'ShowText','on')
 colorbar;
 xlabel('x-axis')
 ylabel('y-axis')
 title_text=sprintf('total gauss seidel iteration @steady state = %d',gs_iter);
 title(title_text)
 end
 if I==1
     ct_j =toc;
     fprintf('%d jacobi iterations for computational time of %0.3g secondsn',k,ct_j)
 elseif I==2 
     ct_gs = toc;
     fprintf('%d gauss seidel iterations for computational time of %0.3g secondsn',k,ct_gs)
 end 
 
end

Answer 1

[Alan Stevens]

Like this? I've used: (T(1:nx,2) - T(1:nx,1))/dy = 5 and rearranged this to get T(1:nx,1). In addition I've set the y values to decrease from ny-1, instead of increasing from 2 and shifted the calculation of T(1:nx,1) to the end of the loop. Check carefully!!

for I=1:2
    tic
nx = 21;   %step size x -direction
ny = 21;   %step size y -direction 
Lx = 1; 
Ly = 1;
x = linspace(0,Lx,nx);
y = linspace(0,Ly,nx);
dx = x(2)-x(1);
dy = y(2)-y(1);
k = 0;
% Intial condition 
T = zeros(nx,ny);
% BC 
T(1:nx,nx) = (sin(pi*x(:,1)));
%T(1:nx,1) = 5;
T(ny,1:ny) = 0;
T(1,1:ny) = (1-y(1,:));
Told = T; 
tol = 1.0e-4;
error = 1;
k = k+1;
% iterative solver 
 iterative_solver = input('iterative_solver(1. jacobi 2. gauss) = ');
 % Jacobi 
 if iterative_solver ==1
     jacobi_iter = 1;
     while(error>tol)
         for i = 2:nx-1
             for j = ny-1:-1:2 %%%%%%%%%%%%%%%%
                  T(i,j) = 0.25*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1) + Told(i,j+1));
             end 
         end 
         T(1:nx,1) = T(1:nx,2) - 5*dy; %%%%%%%%%%%%%
         error = max(max(abs(Told-T)));
         jacobi_iter = jacobi_iter+1;
         Told=T;
     end 
     contourf(x,y,T,'ShowText','on')
     colorbar;
      title_text=sprintf('total jacobi iteration = %d @steady state',jacobi_iter);
 title(title_text)
 xlabel('x-axis');
 ylabel('y-axis');
 end 
 % Gauss 
 if iterative_solver ==2
 gs_iter =1;
 while (error>tol)
     for i = 2:nx-1
         for j = ny-1:-1:2 %%%%%%%%%%%%%%%%%%%
             T(i,j) = (0.25)*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1));
         end 
     end 
     T(1:nx,1) = T(1:nx,2) - 5*dy; %%%%%%%%%%%%%%%%%%%
     error = max(max(abs(Told-T)));
     gs_iter=gs_iter+1;
     Told = T;
 end 
 contourf(x,y,T,'ShowText','on')
 colorbar;
 xlabel('x-axis')
 ylabel('y-axis')
 title_text=sprintf('total gauss seidel iteration @steady state = %d',gs_iter);
 title(title_text)
 end
 if I==1
     ct_j =toc;
     fprintf('%d jacobi iterations for computational time of %0.3g secondsn',k,ct_j)
 elseif I==2 
     ct_gs = toc;
     fprintf('%d gauss seidel iterations for computational time of %0.3g secondsn',k,ct_gs)
 end 
 
end