Sparse matrix with diffusion equation

疎行列を使って一気に解く

 

def RDE():
    dtnum = 4801 #時間をいくつに分割して計算するか 一の位を1とした方が良い 時刻tのfor文の終わりの数字に関係あり
    dxnum = 101 #xとyをいくつに分割して計算するか
    dynum = 101

    thick = 10 #x方向の大きさ
    width = 10 #y方向の大きさ
    time_calc = 48 #計算する時間
    #beta = 0.1 #上の式中ではラムダ 温度拡散率 lambda = 熱伝導率/(密度*比熱)
    
    #constant for PDE
    D = 0.001
    
    """計算準備"""
    #空の解を用意
    solution_k1 = np.zeros([dxnum,dynum])
    solution_k2 = np.zeros([dxnum,dynum]) + 100 * numpy.random.random*1 #solution_k2の初期値が初期条件になる
    
    #境界条件
    #今回はNeumann boundary conditionを作りたい
    
    
    dt = time_calc / (dtnum - 1)
    dx = thick / (dxnum - 1)
    dy = width / (dynum - 1)
    r_x = (dt / (2 * dx ** 2))
    r_y = (dt / (2 * dy ** 2))
    
    
    """Ax=bのAを用意"""
    a_matrix_1 = np.identity(dxnum) * (1 + 2 * D * r_x) \
        + np.eye(dxnum, k=1) * (- D * r_x) \
            + np.eye(dxnum, k=-1) * (- D * r_x)
    a_matrix_1[0,0] /= 2
    a_matrix_1[-1,-1] /= 2
    
    a_matrix_2 = np.identity(dynum) * (1 + 2 * D * r_y) \
        + np.eye(dynum, k=1) * (- D * r_y) \
            + np.eye(dynum, k=-1) * (- D * r_y)
    a_matrix_2[0,0] /= 2
    a_matrix_2[-1,-1] /= 2
    
    #疎行列を格納 csr方式
    a_matrix_csr1 = csr_matrix(a_matrix_1)
    a_matrix_csr2 = csr_matrix(a_matrix_2)
    
    a_matrix_arr1 = (a_matrix_csr1 for i in range(dynum))
    a_matrix_arr2 = (a_matrix_csr2 for i in range(dxnum))
    
    #さらに大きな行列を作る
    a_matrix_csr_agg1 = csr_matrix(block_diag(a_matrix_arr1))
    a_matrix_csr_agg2 = csr_matrix(block_diag(a_matrix_arr2))
    
    #ADI法ではk+1時刻とk+2時刻を計算するので、for文の回数は半分になる
    #number = Decimal(str(dtnum/2)).quantize(Decimal('0'), rounding=ROUND_HALF_UP)
    number = int(dtnum/2)
    
    solution = np.zeros([dxnum,dynum,number]) #解を代入する行列を作成
    
    #ADI法の計算
    
    for k in range(number): #時刻tについて
        #大きな行列で一気に計算を済ませてしまう。
        #for loopが減るので、計算時間が1/3になった。
        #一回目の拡散
        template_1 = np.zeros*2
        template_1[1 : dxnum + 1, 1 : dynum + 1] = solution_k2
        template_1[:, 0] = template_1[:, 2]
        template_1[:, -1] = template_1[:, -3]
        
        b_array_x = solution_k2 \
            + D * r_y * (template_1[1: dxnum + 1, 0: dynum] \
                - 2 * solution_k2 \
                    + template_1[1: dxnum + 1, 2: dynum + 2]) \
                + dt * (1 / 2) * function(solution_k2)
        b_array_x[0, :] /= 2
        b_array_x[-1,:] /= 2
        
        #一気に解くために、1次元配列にする。
        #ただし、足し合わせる方向性を変えるために、一度転置行列にしている
        b_array_x_1d = b_array_x.T.reshape(1, dxnum * dynum)[0]
        temp_solve1 = spla.spsolve(a_matrix_csr_agg1, b_array_x_1d)
        #得られた解を2次元配列に戻す。
        solution_k1 = temp_solve1.reshape(dynum,dxnum).T
        
        
        template_2 = np.zeros*3
        template_2[1:dxnum+1, 1:dynum+1] = solution_k1
        template_2[0,:] = template_2[2,:]
        template_2[-1,:] = template_2[-3,:]
        
        b_array_y = solution_k1 \
            + D * r_x * (template_2[0: dxnum, 1: dynum + 1] \
                - 2 * solution_k1 \
                    + template_2[2: dxnum + 2, 1: dynum + 1]) \
                + dt * (1 / 2) * function(solution_k1)
        b_array_y[:, 0] /= 2
        b_array_y[:,-1] /= 2
        
        #一気に解くために1次元配列にする。
        
        b_array_y_1d = b_array_y.reshape(1 , dxnum * dynum)[0]
        temp_solve2 = spla.spsolve(a_matrix_csr_agg2, b_array_y_1d)
        #得られた解を2次元配列に戻す。
        solution_k2 = temp_solve2.reshape(dxnum,dynum)
        
        solution[:,:,k] = solution_k2
        
    return solution

 

 

 

 

 

*1:dxnum, dynum

*2:dxnum + 2, dynum + 2

*3:dxnum + 2, dynum + 2