熱傳模組的建立（續） 穩態邊界所有情形

     case 2
                disp('Now we will divide the 2-D plane for FDE')
        %too choose the divide element by setting the side
        r=input('Please input m row of the matrix you want ==>   ');
        c=input('Please input n column of the matrix you want ==>   ');
        zeros(r,c);
        for t=1:r*c
           v(t)=t;
        end
        V=reshape(v,r,c);
        %set up a matrix to number every element of the matrix
        N=r*c;
        A=zeros(N);
        %after observing the FDE matrix, derive down rule code
        for i=1:N
          A(i,i)=-4;
          if i+r0
              A(i,i-r)=1;
          end
          if i+1> ')
        H=input('The convection coeffcient is>> ')
        K=input('The conduction coeffcient is>> ') 
        L=input('The length of element>> ')
               %First of all,we deal the boundary exclude the corner
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
               %Top boundary
               switch Bu
                   case 1
                       Tu=input('Temperature of top boundary>> ')
                        for j=2:c-1  %it may be a bug here
                            B((j-1)*r+1,1)=-Tu;
                        end
                   case 2
                       Tu=input('The surrounding temperature of the top>> ')
                    
                        for j=2:c-1
                            B((j-1)*r+1,1)=-(2*H*Tu*L/K+G*L^2/K);
                            A((j-1)*r+1,(j-1)*r+1)=-(4+2*H*L/K);% retune the matrix
                            A((j-1)*r+1,(j-1)*r+2)=2;
                        end
                   case 3
                       Qu=input('-dT/dy  (toward the center is positive)  =')
                       for j=2:c-1
                            B((j-1)*r+1,1)=-(2*Qu*L/K+G*L^2/K);
                            A((j-1)*r+1,(j-1)*r+1)=-4 ;% retune the matrix
                            A((j-1)*r+1,(j-1)*r+2)=2;  
                       end
               end
              
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
               %left boundary
               switch Bl
                   case 1
                       Tl=input('Temperature of left boundary>> ')
                       for k=2:r-1
                           B(k,1)=-Tl;
                       end
                   case 2
                       Tl=input('The surrounding temperature of the left>> ')
                       for k=2:r-1
                           B(k,1)=-(2*H*Tl*L/K+G*L^2/K);
                           A(k,k)=-(4+2*H*L/K);
                           A(k,k+r)=2;
                       end
                      
                   case 3
                       Ql=input('-dT/dy  (toward the center is positive)  =')
                       for k=2:r-1
                           B(k,1)=-(2*Ql*L+G*L^2);
                           A(k,k)=-4;
                           A(k,k+r)=2;
                       end
               end
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
              
                    %down boundary
                   
               switch Bd
                   case 1
                       Td=input('Temperature of down boundary>> ')
                       for h=2:c-1
                         B(r*h,1)=-Td;
                       end
                   case 2
                       Td=input('The surrounding temperature of the bottom>> ')
                       for h=2:c-1
                         B(r*h,1)=-(2*H*Td*L/K+G*L^2/K);% retune the matrix
                         A(r*h,r*h)=-(4+2*H*L/K);
                         A(r*h,r*h-1)=2;
                       end
                   case 3
                       Qd=input('-dT/dy  (toward the center is positive)  =')
                       for h=2:c-1
                         B(r*h,1)=-(2*Qd*L/K+G*L^2/K);% retune the matrix
                         A(r*h,r*h)=-4;
                         A(r*h,r*h-1)=2;
                       end
               end
              
              
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
                    %right boundary
               switch Br
                   case 1
                       Tr=input('Temperature of right boundary>> ')
                       for p=2:r-1
                         B(r*(c-1)+p,1)=-Tr;
                       end
                   case 2
                       Tr=input('The surrounding temperature of the right>> ')
                       for p=2:r-1
                           B(r*(c-1)+p,1)=-(2*H*Tr*L/K+G*L^2/K);
                           A(r*(c-1)+p,r*(c-1)+p)=-(4+2*H*L/K);
                           A(r*(c-1)+p,r*(c-2)+p)=2;
                       end
                      
                   case 3
                       Qr=input('-dT/dy  (toward the center is positive)  =')
                       for p=2:r-1
                           B(r*(c-1)+p,1)=-(2*H*Qr*L/K+G*L^2/K);
                           A(r*(c-1)+p,r*(c-1)+p)=-4;
                           A(r*(c-1)+p,r*(c-2)+p)=2;
                       end
               end
              
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
      %four corner's temp need another expression (6 situation)
   
      % the first angle B(1,1)
       if Bu==1&Bl==1     %both constant Temp
          B(1,1)=-(G*L^2/(2*K)+(Tu+Tl));
          A(1,1)=-4;
       elseif Bu==2&Bl==2  %both convection
           B(1,1)=-(G*L^2/(2*K)+H*L*(Tu+Tl)/K);
           A(1,1)=-(2+2*H*L/K);
       elseif Bu==3&Bl==3  %both heat flux or insulated
           B(1,1)=-(G*L^2/(2*K)+2*(Qu+Ql)*L/K);
           A(1,1)=-2;
       elseif Bu==1&Bl==2 
           B(1,1)=-(G*L^2/(2*K)+H*L*Tl/K+Tu);
           A(1,1)=-(3+H*L);
       elseif Bu==2&Bl==1
           B(1,1)=-(G*L^2/(2*K)+H*L*Tu/K+Tl);
           A(1,1)=-(3+H*L);
       elseif Bu==3&Bl==1
           B(1,1)=-(G*L^2/(2*K)+Qu*L/K+Tl);
           A(1,1)=-3;
       elseif Bu==1&Bl==3   
           B(1,1)=-(G*L^2/(2*K)+Ql*L/K+Tu);
           A(1,1)=-3;
       elseif Bu==2&Bl==3
           B(1,1)=-(G*L^2/(2*K)+Ql*L/K+H*L*Tu/K);
           A(1,1)=-(2+H*L/K);
       elseif Bu==3&Bl==2
           B(1,1)=-(G*L^2/(2*K)+Qu*L/K+H*L*Tl/K);
           A(1,1)=-(2+H*L/K);
       end
                
    % the second angle B(r,1)
   
       if Bd==1&Bl==1     %both constant Temp
          B(r,1)=-(G*L^2/(2*K)+(Td+Tl));
          A(r,r)=-4;
       elseif Bd==2&Bl==2  %both convection
           B(r,1)=-(G*L^2/(2*K)+H*L*(Td+Tl)/K);
           A(r,r)=-(2+2*H*L/K);
       elseif Bd==3&Bl==3  %both heat flux or insulated
           B(r,1)=-(G*L^2/(2*K)+2*(Qd+Ql)*L/K);
           A(r,r)=-2;
       elseif Bd==1&Bl==2 
           B(r,1)=-(G*L^2/(2*K)+H*L*Tl/K+Td);
           A(r,r)=-(3+H*L);
       elseif Bd==2&Bl==1
           B(r,1)=-(G*L^2/(2*K)+H*L*Td/K+Tl);
           A(r,r)=-(3+H*L);
       elseif Bd==3&Bl==1
           B(r,1)=-(G*L^2/(2*K)+Qd*L/K+Tl);
           A(r,r)=-3;
       elseif Bd==1&Bl==3   
           B(r,1)=-(G*L^2/(2*K)+Ql*L/K+Td);
           A(r,r)=-3;
       elseif Bd==2&Bl==3
           B(r,1)=-(G*L^2/(2*K)+Ql*L/K+H*L*Td/K);
           A(r,r)=-(2+H*L/K);
       elseif Bd==3&Bl==2
           B(r,1)=-(G*L^2/(2*K)+Qd*L/K+H*L*Tl/K);
           A(r,r)=-(2+H*L/K);
       end
   
       % the third angle B(r*(c-1)+1,1)
       if Bu==1&Br==1     %both constant Temp
          B(r*(c-1)+1,1)=-(G*L^2/(2*K)+(Tu+Tr));
          A(r*(c-1)+1,r*(c-1)+1)=-4;
       elseif Bu==2&Br==2  %both convection
           B(r*(c-1)+1,1)=-(G*L^2/(2*K)+H*L*(Tu+Tr)/K);
           A(r*(c-1)+1,r*(c-1)+1)=-(2+2*H*L/K);
       elseif Bu==3&Br==3  %both heat flux or insulated
           B(r*(c-1)+1,1)=-(G*L^2/(2*K)+2*(Qu+Qr)*L/K);
           A(r*(c-1)+1,r*(c-1)+1)=-2;
       elseif Bu==1&Br==2 
           B(r*(c-1)+1,1)=-(G*L^2/(2*K)+H*L*Tr/K+Tu);
           A(r*(c-1)+1,r*(c-1)+1)=-(3+H*L);
       elseif Bu==2&Br==1
           B(r*(c-1)+1,1)=-(G*L^2/(2*K)+H*L*Tu/K+Tr);
           A(r*(c-1)+1,r*(c-1)+1)=-(3+H*L);
       elseif Bu==3&Br==1
           B(r*(c-1)+1,1)=-(G*L^2/(2*K)+Qu*L/K+Tr);
           A(r*(c-1)+1,r*(c-1)+1)=-3;
       elseif Bu==1&Br==3   
           B(r*(c-1)+1,1)=-(G*L^2/(2*K)+Qr*L/K+Tu);
           A(r*(c-1)+1,r*(c-1)+1)=-3;
       elseif Bu==2&Br==3
           B(r*(c-1)+1,1)=-(G*L^2/(2*K)+Qr*L/K+H*L*Tu/K);
           A(r*(c-1)+1,r*(c-1)+1)=-(2+H*L/K);
       elseif Bu==3&Br==2
           B(r*(c-1)+1,1)=-(G*L^2/(2*K)+Qu*L/K+H*L*Tr/K);
           A(r*(c-1)+1,r*(c-1)+1)=-(2+H*L/K);
       end
  
       % the fourth angle B(r*c,1)
       if Bd==1&Br==1     %both constant Temp
           B(r*c,1)=-(G*L^2/(2*K)+(Td+Tr));
           A(r*c,r*c)=-4;
       elseif Bd==2&Br==2  %both convection
           B(r*c,1)=-(G*L^2/(2*K)+H*L*(Td+Tr)/K);
           A(r*c,r*c)=-(2+2*H*L/K);
       elseif Bd==3&Br==3  %both heat flux or insulated
           B(r*c,1)=-(G*L^2/(2*K)+2*(Qd+Qr)*L/K);
           A(r*c,r*c)=-2;
       elseif Bd==1&Br==2 
           B(r*c,1)=-(G*L^2/(2*K)+H*L*Tr/K+Td);
           A(r*c,r*c)=-(3+H*L);
       elseif Bd==2&Br==1
           B(r*c,1)=-(G*L^2/(2*K)+H*L*Td/K+Tr);
           A(r*c,r*c)=-(3+H*L);
       elseif Bd==3&Br==1
           B(r*c,1)=-(G*L^2/(2*K)+Qd*L/K+Tr);
           A(r*c,r*c)=-3;
       elseif Bd==1&Br==3   
           B(r*c,1)=-(G*L^2/(2*K)+Qr*L/K+Td);
           A(r*c,r*c)=-3;
       elseif Bd==2&Br==3
           B(r*c,1)=-(G*L^2/(2*K)+Qr*L/K+H*L*Td/K);
           A(r*c,r*c)=-(2+H*L/K);
       elseif Bd==3&Br==2
           B(r*c,1)=-(G*L^2/(2*K)+Qd*L/K+H*L*Tr/K);
           A(r*c,r*c)=-(2+H*L/K);
       end
            
        %calculate and plot the distribution of temp
        X=inv(A)*B;
        Q=reshape(transpose(X),r,c);
        figure(1)
        mesh(Q);
        xlabel('Top boundary')
        ylabel('Left boundary')
        zlabel('Temperature(℃)')
        figure(2)
        surf(Q);
        xlabel('Top boundary')
        ylabel('Left boundary')
        zlabel('Temperature(℃)')