疎行列を使って一気に解く
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が減るので、計算時間が１/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
