function mtx_A = coefficient_A(n1, n2, n3, sigma, vec_k) cubic_length = 1; period_vec1 = cubic_length * [ 1; 0; 0 ]; period_vec2 = cubic_length * [ 0; 1; 0 ]; period_vec3 = cubic_length * [ 0; 0; 1 ]; length_x = cubic_length; length_y = cubic_length; length_z = cubic_length; delta = [ length_x / n1; length_y / n2; length_z / n3]; mtx_G = Sparse_mtx_G(n1,n2,n3,delta,period_vec1,period_vec2,period_vec3,vec_k); mtx_A = mtx_G' * mtx_G + sigma * speye(3*n1*n2*n3); end % =========================== function [mtx_G] = Sparse_mtx_G(n_1,n_2,n_3,delta,period_vec1,period_vec2,period_vec3,vec_k) n12 = n_1 * n_2; n123 = n12 * n_3; [i_row_1, j_col_1, val_1] = Gamma_1(n_1, n_2, n_3, period_vec1, vec_k); [i_row_2, j_col_2, val_2] = Gamma_2(n_1, n_2, n_3, period_vec2, vec_k); [i_row_3, j_col_3, val_3] = Gamma_3(n_1,n_2,n_3, period_vec3, vec_k); leng_1 = length(i_row_1); leng_2 = length(i_row_2); leng_3 = length(i_row_3); leng = 2 * ( leng_1 + leng_2 + leng_3 ); i_row = zeros(leng,1); j_col = zeros(leng,1); val = zeros(leng,1); val_1 = val_1 / delta(1); val_2 = val_2 / delta(2); val_3 = val_3 / delta(3); % mtx_C3 = sparse(i_row_3, j_col_3, val_3); % First block rows i_row(1:leng_3,1) = i_row_3; j_col(1:leng_3,1) = n123 + j_col_3; val (1:leng_3,1) = val_3; ii = leng_3 + leng_2; i_row(leng_3+1:ii,1) = i_row_2; j_col(leng_3+1:ii,1) = 2 * n123 + j_col_2; val (leng_3+1:ii,1) = - val_2; % Second block rows i2 = ii + leng_3; i_row(ii+1:i2,1) = n123 + i_row_3; j_col(ii+1:i2,1) = j_col_3; val (ii+1:i2,1) = - val_3; ii = i2 + leng_1; i_row(i2+1:ii,1) = n123 + i_row_1; j_col(i2+1:ii,1) = 2 * n123 + j_col_1; val (i2+1:ii,1) = val_1; % Third block rows i2 = ii + leng_2; i_row(ii+1:i2,1) = 2 * n123 + i_row_2; j_col(ii+1:i2,1) = j_col_2; val (ii+1:i2,1) = val_2; ii = i2 + leng_1; i_row(i2+1:ii,1) = 2 * n123 + i_row_1; j_col(i2+1:ii,1) = n123 + j_col_1; val (i2+1:ii,1) = - val_1; mtx_G = sparse(i_row, j_col, val); end % ========================================= function [i_row, j_col, val] = Gamma_1(n_1, n_2, n_3, period_vec1, vec_k) mtx_C_1 = gen_mtx_C_1(n_1,period_vec1,vec_k); [i_row_tmp, j_col_tmp, val_tmp] = find(mtx_C_1); leng = length(i_row_tmp); i_row = zeros(leng*n_2*n_3,1); j_col = zeros(leng*n_2*n_3,1); val = zeros(leng*n_2*n_3,1); for ii = 1:n_2*n_3 jj = ( ii - 1 ) * n_1; kk = ( ii - 1 ) * leng; i_row(kk+1:kk+leng,1) = jj + i_row_tmp; j_col(kk+1:kk+leng,1) = jj + j_col_tmp; val (kk+1:kk+leng,1) = val_tmp; end end % ======================================== function [mtx_C_1] = gen_mtx_C_1(n_1,period_vec1,vec_k) % % Generate matrix K_1^{*} % mtx_C_1 = -eye(n_1); for ii = 1:n_1-1 mtx_C_1(ii,ii+1) = 1; end tmp = 2 * pi * (vec_k' * period_vec1) * 1i; mtx_C_1(n_1,1) = exp(tmp); end % ======================================== function [i_row, j_col, val] = Gamma_2(n_1, n_2, n_3, period_vec2, vec_k) [i_row_tmp, j_col_tmp, val_tmp] = gen_mtx_C_2(n_1,n_2,period_vec2,vec_k); leng = length(i_row_tmp); i_row = zeros(leng*n_3,1); j_col = zeros(leng*n_3,1); val = zeros(leng*n_3,1); for ii = 1:n_3 jj = ( ii - 1 ) * n_1 * n_2; kk = ( ii - 1 ) * leng; i_row(kk+1:kk+leng,1) = jj + i_row_tmp; j_col(kk+1:kk+leng,1) = jj + j_col_tmp; val (kk+1:kk+leng,1) = val_tmp; end end % ================================================= function [i_row, j_col, val] = gen_mtx_C_2(n_1,n_2,period_vec2,vec_k) % % Generate matrix K_2^{*} % i_row = zeros(2*n_1*n_2,1); j_col = zeros(2*n_1*n_2,1); val = zeros(2*n_1*n_2,1); for jj = 1:n_2-1 ii = ( jj - 1 ) * n_1; kk = 2 * ii; i_row(kk+1:kk+n_1,1) = ii+1:ii+n_1; j_col(kk+1:kk+n_1,1) = ii+1:ii+n_1; val (kk+1:kk+n_1,1) = -ones(n_1,1); i_row(kk+n_1+1:kk+2*n_1,1) = ii+1:ii+n_1; j_col(kk+n_1+1:kk+2*n_1,1) = ii+n_1+1:ii+2*n_1; val (kk+n_1+1:kk+2*n_1,1) = ones(n_1,1); end ii = ( n_2 - 1 ) * n_1; kk = 2 * ii; i_row(kk+1:kk+n_1,1) = ii+1:ii+n_1; j_col(kk+1:kk+n_1,1) = ii+1:ii+n_1; val (kk+1:kk+n_1,1) = -ones(n_1,1); i_row(kk+n_1+1:kk+2*n_1,1) = ii+1:ii+n_1; j_col(kk+n_1+1:kk+2*n_1,1) = 1:n_1; tmp = 2 * pi * (vec_k' * period_vec2) * 1i; val (kk+n_1+1:kk+2*n_1,1) = exp(tmp) * ones(n_1,1); end % ================================================= function [i_row, j_col, val] = Gamma_3(n_1,n_2,n_3, period_vec3, vec_k) n12 = n_1 * n_2; i_row = zeros(2*n_1*n_2*n_3,1); j_col = zeros(2*n_1*n_2*n_3,1); val = zeros(2*n_1*n_2*n_3,1); for jj = 1:n_3-1 ii = ( jj - 1 ) * n12; kk = 2 * ii; i_row(kk+1:kk+n12,1) = ii+1:ii+n12; j_col(kk+1:kk+n12,1) = ii+1:ii+n12; val (kk+1:kk+n12,1) = -ones(n12,1); i_row(kk+n12+1:kk+2*n12,1) = ii+1:ii+n12; j_col(kk+n12+1:kk+2*n12,1) = ii+n12+1:ii+2*n12; val (kk+n12+1:kk+2*n12,1) = ones(n12,1); end ii = ( n_3 - 1 ) * n12; kk = 2 * ii; i_row(kk+1:kk+n12,1) = ii+1:ii+n12; j_col(kk+1:kk+n12,1) = ii+1:ii+n12; val (kk+1:kk+n12,1) = -ones(n12,1); i_row(kk+n12+1:kk+2*n12,1) = ii+1:ii+n12; j_col(kk+n12+1:kk+2*n12,1) = 1:n12; tmp = 2 * pi * (vec_k' * period_vec3) * 1i; val (kk+n12+1:kk+2*n12,1) = exp(tmp) * ones(n12,1); end % ===============================================