雲濤夢筆: 2006

繼上次的穩態後，這次繼續寫非穩態的情形
unsteady的情形的主要概念在於時間步的產生，從原先的位置網格分析後
這次多了時間網格分析，主要在於說位置網格分析是建立於時間網格分析上
，如此牽涉矩陣包矩陣的情形，故而採用了cell的指令去做回圈設計。

而數值法依然採用顯明法，在此介紹兩種方法。

Explicit method

Explicit method要求ITS（integration time step）的值非常小（否則不容易收斂），較適合應用在非常短暫的衝擊、高度非線性的分析。(由已知求未知)

Implicit method

Implicit method則可以容許較大的ITS值，但是對於短暫的衝擊、高度非線性的分析則常有收斂的困難。
（由未知求未知）

程式碼如下：


        disp('Now we will divide the 2-D plane for FDE')
        %to 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 ==>   ');
        dt=input('time step ==> ');
        lo=input('density==> ');
        K=input('conduction coefficient ==>  ');
        l=input('mesh length==> ');
        ts=input('total time==>  ');
        C=input('specific heat== >');
        Ti=input('initial temperature==> ');
        arfa=K/(lo*C);
        tt=[1:dt:ts];
        tau=(arfa/l^2)*dt;
        node=length(tt)

        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)=1-4*tau;
          if i+r              A(i,i+r)=tau;
          end
          if i-r>0
              A(i,i-r)=tau;
          end
          if i+1           A(i+1,i)=tau;
           A(i,i+1)=tau;
          end
        end
        for j=1:c-1
            A(j*r,j*r+1)=0;
            A(j*r+1,j*r)=0;
        end
        
        
        
        B=zeros(r*c,1);
        %Top boundary
        Tu=input('Temperature of top boundary>> ');
        for j=1:c
            B((j-1)*r+1,1)=Tu; 
        end

        %left boundary
        Tl=input('Temperature of left boundary>> ');
        for k=1:r
            B(k,1)=Tl;
        end

        %down boundary
        Td=input('Temperature of down boundary>> ');
        for h=1:c
            B(r*h,1)=Td;
        end

        %right boundary
        Tr=input('Temperature of right boundary>> ');
        for p=1:r
            B(r*(c-1)+p,1)=Tr;
        end

        %four corner's temp need another expression
        B(1,1)=(Tu+Tl);
        B(r,1)=(Td+Tl);
        B(r*(c-1)+1,1)=(Tu+Tr);
        B(r*c,1)=(Td+Tr);    
        B=tau*B;
        Tim=Ti*ones(r*c,1);
        H=cell(node,1);%%for every timematrix placing
        for i=1:node
            R=A*Tim+B;
            Tim=R;
            H(i,1)={R};
        end
        
        
        %%show the node and let the user decide the time they want to see
        %%the situation
        node
        nd=input('which node you want to see ==>');
        
        Q=reshape(transpose(H{nd,1}),r,c);
        figure(1)
        mesh(Q);
        xlabel('X index(△l)')
        ylabel('Y index(△l)')
        zlabel('Temperature(℃)')
        for i=1:size(Q,1)
        for j=1:size(Q,2)
                h=text(j, i, Q(i,j), num2str(V(i, j)));  
                % numbering (if use Q(i,j)to replace V(i,j) will show the value of every node)
                set(h, 'hori', 'center', 'vertical', 'bottom', 'color', 'r');  % h's position and color
           end
         end
        figure(2)
        surf(Q);
        xlabel('X index(△l)')
        ylabel('X index(△l)')
        zlabel('Temperature(℃)')
        
      %% to make a movie by pause  
        for i= 1:node
             Q=reshape(transpose(H{i,1}),r,c);
             surf(Q);
             xlabel('X index(△l)')
             ylabel('X index(△l)')
             zlabel('Temperature(℃)')
             pause(0.001)
        end

結果：
先輸入條件
Please input m row of the matrix you want ==> 8
Please input n column of the matrix you want ==> 7
time step ==> 0.02
density==> 750
conduction coefficient ==> 70
mesh length==> 0.05
total time==> 30
specific heat== >980
initial temperature==> 5