Demo entry 6769992

test

   

Submitted by anonymous on Nov 11, 2018 at 08:14
Language: Python. Code size: 3.2 kB.

# -*- coding: utf-8 -*-
import math
import numpy
from matplotlib import pyplot as plt
# 给初始条件
a = 66.00526
Sigma_t = 0.05
Sigma_s = 0.03
nv_Sigma_f = 0.0225
# S4
mu = [-0.8611363116, -0.3399810436, 0.3399810436, 0.8611363116]
w = [0.3478548451, 0.6521451549, 0.6521451549, 0.3478548451]
accuracy = 1e-5 #精确度
# 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 = 200 # 段数

array_number = K + 1 #数组维数

delta_x = a / K

x = []
x_i = 0
for i in range(array_number):
    x.append(x_i)
    x_i += delta_x
k_avg = [0]
def main():
    S = numpy.ones([array_number])
    phi = numpy.zeros([len(mu),array_number])
    S_his = []
    cycle = 0
    S_his.append(S[K - 1])
    while 1:
        # 算通量
        for m in range(len(mu)):
            if  mu[m] < 0:
                k = K
                phi[m][K] = 0
                while (k > 0):
                    k = k - 2
                    phi[m][k] = (2 + (Sigma_t * 2 * delta_x / mu[m]))/(2 - (Sigma_t * 2 * delta_x / mu[m])) * phi[m][k + 2] - 2 * (Sigma_t * 2 * delta_x / mu[m])/(2 - (Sigma_t * 2 * delta_x / mu[m])) * S[k + 1] / Sigma_t
                    phi[m][k + 1] = (phi[m][k] + phi[m][k + 2]) / 2
            if mu[m] > 0:
                k = 0
                phi[m][0] = phi[len(mu)-m-1][0]
                while (k < K):
                    k = k + 2
                    phi[m][k] = (2 - (Sigma_t * 2 * delta_x / mu[m]))/(2 + (Sigma_t * 2 * delta_x / mu[m])) * phi[m][k - 2] + 2 * (Sigma_t * 2 * delta_x / mu[m])/(2 + (Sigma_t * 2* delta_x / mu[m])) * S[k - 1] / Sigma_t
                    phi[m][k - 1] = (phi[m][k] + phi[m][k - 2]) / 2

        # 更新S
        k = 0
        while (k < K):
            k = k + 2
            sum = 0
            for m in range(len(mu)):
                sum = sum + w[m] * (phi[m][k - 2] + phi[m][k])
            S[k - 1] = ((Sigma_s + nv_Sigma_f) / 4) * sum
        
        cycle = cycle + 1
        #记录每一次迭代的S值
        S_his.append(S[K - 1]) 
        #计算k
        k_avg.append(S_his[cycle] / S_his[cycle - 1])
        #打印K值
        #判断是否达到精确度
        if abs((S_his[cycle] - S_his[cycle - 1])) < accuracy:
            print k_avg[cycle]
            break

    # 画分布
    for i in range(len(mu)):
        plt.plot(x, phi[i])
        plt.title('phi')

    plt.show()
if __name__ == '__main__':
    main()

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).