SN .zip_SN_SN方法

# -*- coding: utf-8 -*- import math import numpy from matplotlib import pyplot as plt R = 145.5436 Sigma_t = 0.05 Sigma_s = 0.03 nv_Sigma_f = 0.0225 accuracy = 1e-5 # S4 mu = [-0.86113, -0.33998, 0.33998, 0.86113] w = [0.34785, 0.65214, 0.65214, 0.34785] # S6 #mu = [-0.93246, -0.66120, -0.23861, 0.23861, 0.66120, 0.93246] #w = [0.17132, 0.36076, 0.46791, 0.46791, 0.36076, 0.17132] # S12 #mu = [-0.9815606342, -0.9041172564, -0.7699026742, -0.5873179543, -0.3678314990, -0.1252334085, 0.1252334085, 0.3678314990, 0.5873179543, 0.7699026742, 0.9041172564, 0.9815606342] #w = [0.0471753364, 0.1069393260, 0.1600783285, 0.2031674267, 0.2334925365, 0.2491470458, 0.2491470458, 0.2334925365, 0.2031674267, 0.1600783285, 0.1069393260, 0.0471753364] # S16 #mu = [-0.98940, -0.94457, -0.86563, -0.75540, -0.61787, -0.45801, -0.28160, -0.09601, 0.09601, 0.28160, 0.45801, \ # 0.61787, 0.75540, 0.86563, 0.94457, 0.98940] #w = [0.02715, 0.06225, 0.09515, 0.12462, 0.14959, 0.16915, 0.18260, 0.18945, 0.18945, 0.18260, 0.16915, 0.14959, \ # 0.12462, 0.09515, 0.06225, 0.02715] K = 10000 # 段数 array_number = K + 1 #数组维数 delta_r = R / K r = [] r_i = 0 # 定义初始量 for i in range(array_number): r.append(r_i) r_i += delta_r i = 0 A = numpy.zeros([array_number]) while i < K: i += 1 A[i] = r[i] * r[i] i = 0 V = numpy.zeros([K]) while i < K: V[i] = (r[i + 1]**3 - r[i]**3) / 3 i += 1 alpha = numpy.zeros([len(mu) + 1]) for i in range(len(mu)): alpha[i + 1] = - (mu[i] * w[i]) + alpha[i] k_avg = [0] print r print A print V print alpha[4] S_his = [] def main(): a = numpy.ones([2 * len(mu) + 1,array_number + K]) S = numpy.ones([K]) phi = numpy.ones([2 * len(mu) + 1,array_number + K]) cycle = 0 S_his.append(S[1]) while 1: # 算源 k = 0 while k < K: sum = 0 j = 2 while j < (2 * len(mu) + 1): sum = sum + phi[j - 1][2 * k + 1] * w[j / 2 - 1] j += 2 S[k] = ((Sigma_s + nv_Sigma_f) / 2) * sum k += 1 # 算通量 i = 0 k = K phi[0][2 * K] = 0 while k > 0: k -= 1 phi[0][2 * k + 1] = (V[k] * S[k] + (A[k + 1] + A[k]) * phi[0][2 * k + 2]) / (2 * A[k] + Sigma_t * V[k]) phi[0][2 * k] = 2 * phi[0][2 * k + 1] - phi[0][2 * k + 2] for i in range(2,(2 * len(mu) + 1),2): if (mu[i / 2 - 1] < 0): k = K while k > 1: k -= 1 phi[i - 1][2 * k + 1] = (V[k] * S[k] + abs(mu[i / 2 - 1]) * (A[k + 1] + A[k]) * phi[i - 1][2 * k + 2] + ((A[k + 1] - A[k]) * (alpha[i / 2] + alpha[i / 2 - 1]) * phi[i - 2][2 * k + 1] / w[i / 2 - 1]))\ / (abs(mu[i / 2 - 1]) * (A[k + 1] + A[k]) + ((A[k + 1] - A[k]) * (alpha[i / 2] + alpha[i / 2 - 1]) / w[i / 2 - 1]) + Sigma_t * V[k]) phi[i - 1][2 * k] = 2 * phi[i - 1][2 * k + 1] - phi[i - 1][2 * k + 2] phi[i][2 * k + 1] = 2 * phi[i - 1][2 * k + 1] - phi[i - 2][2 * k + 1] if (mu[i / 2 - 1] > 0): k = 0 phi[i - 1][0] = phi[2 * len(mu) - (i - 1)][0] while k < K: phi[i - 1][2 * k + 1] = (V[k] * S[k] + abs(mu[i / 2 - 1]) * (A[k + 1] + A[k]) * phi[i - 1][2 * k] + ((A[k + 1] - A[k]) * (alpha[i / 2] + alpha[i / 2 - 1]) * phi[i - 2][2 * k + 1] / w[i / 2 - 1]))\ / (abs(mu[i / 2 - 1]) * (A[k + 1] + A[k]) + ((A[k + 1] - A[k]) * (alpha[i / 2] + alpha[i / 2 - 1]) / w[i / 2 - 1]) + Sigma_t * V[k]) phi[i - 1][2 * k + 2] = 2 * phi[i - 1][2 * k + 1] - phi[i - 1][2 * k] phi[i][2 * k + 1] = 2 * phi[i - 1][2 * k + 1] - phi[i - 2][2 * k + 1] k += 1 cycle = cycle + 1 S_his.append(S[1]) print phi[3][K - 1] k_avg.append(abs(S_his[cycle]/S_his[cycle - 1])) if abs(S_his[cycle] - S_his[cycle - 1]) < accuracy: print k_avg[cycle] break i = 0 #while i < (2 * len(mu) + 1): ''' phi_new = numpy.zeros([2 * len(mu) + 1, K]) r_new = numpy.zeros([K]) for i in range(0, 2 * len(mu) + 1): phi_new[] for index in range(len(r) - 1): r_new[index] = ((r[index] + r[index + 1]) / 2) phi_new[i][index] = phi[i][2 * index + 1] plt.plot(r_new,phi_new[i]) # i += 2 ''' phi_new = numpy.zeros([K]) r_new = numpy.zeros([K]) for j in range(0, K): phi_new[j] = phi[5][2 * j + 1] r_new[j] = (r[j] + r[j + 1]) / 2 plt.plot(r_new, phi_new) plt.show() if __name__ == '__main__': main()