python——反应性方程图解法

本文专业性较强，涉核专业学生可能理解起来轻松点。使用python实现反应性方程的图解法，采用变步长思想，图解过程体现了刚性问题的本质，欢迎各位讨论留言。

代码如下：

#-*- coding:utf-8 -*-
#May 16th,2021 on Sun
#Luo·Y·T

#import module
import matplotlib.pyplot as plt
import numpy as np

#initial values
l = 1e-4 #prompt neutron life, unit:s
β= np.array([0.000247, 0.001385, 0.001222, 0.002645, 0.000832, 0.000169]) #delayed neutron fraction, unit:none
t = np.array([78.64, 31.51, 8.66, 3.22, 0.716, 0.258]) #6 delayed neutrons' life, unit:s
λ= [] #void for decay constant
for i in range(len(t)): 
    λ.append(1/t[i]) #decayed constant
λ= np.array(λ) #transform list into matrix

#formula of reactivity
def ρ_formula(ω): 
    SUM = 0
    Σ = []
    for i in range(len(λ)):
        Σ.append(β[i]*ω/(ω+λ[i]))
        SUM += Σ[i]
    ρ=l*ω/(1+l*ω)+(1/(1+l*ω))*SUM
    return ρ
x,y=list(range(11)),list(range(11))
for i in range(11):
    x[i],y[i]=[],[]
    
#change the step-size 
x[0]= np.linspace(-10000000,-1/l,10000)
x[1]= np.linspace(-1/l, -4, 1000)
x[2]= np.linspace(-4, -λ[5], 1000)
x[3]= np.linspace(-λ[5],-λ[4], 20000)
x[4]= np.linspace(-λ[4],-λ[3], 5000)
x[5]= np.linspace(-λ[3],-λ[2], 5000)
x[6]= np.linspace(-λ[2],-λ[1], 6000)
x[7]= np.linspace(-λ[1],-λ[0], 50000)
x[8]= np.linspace(-λ[0],-0.01, 3000)
x[9]= np.linspace(-0.01,0,3)
x[10]= np.linspace(0,1000000, 1000)

#get the corresponding ρ
for i in range(11):
    x[i]=x[i][1:-1]
    for j in range(len(x[i])):
        y[i].append(ρ_formula(x[i][j]))
    if i==10:
        x[i]=list(np.log10(abs(x[i]))) #logarithmic coordinate
    else:
        x[i]=list(-np.log10(abs(x[i]))) 

#combine domain for 8 lines
x[0],y[0]=x[0],y[0] #line1       
x[1],y[1]=x[1]+x[2],y[1]+y[2] #line2
x[2],y[2]=x[3],y[3] #line3
x[3],y[3]=x[4],y[4] #line4
x[4],y[4]=x[5],y[5] #line5
x[5],y[5]=x[6],y[6] #line6
x[6],y[6]=x[7],y[7] #line7
x[7],y[7]=x[8]+x[9]+x[10],y[8]+y[9]+y[10] #line8

#draw 8 lines
for i in range(8):
    plt.plot(x[i],y[i],color='blue',linewidth=1)
    
#line:ρ=+1
x=np.linspace(-8,10,2)
y=[1,1]
plt.plot(x,y,linestyle='--',color='black',linewidth=1)
#line:ρ=+0.5
y=[0.5,0.5]
plt.plot(x,y,linestyle='--',color='black',linewidth=1)
#line:ρ=-0.5
y=[-0.5,-0.5]
plt.plot(x,y,linestyle='--',color='black',linewidth=1)

#optimization        
plt.xlim(-6,6)
plt.ylim(-3,3)
plt.xlabel('ω[logarithmic coordinates]')
plt.ylabel('ρ')
plt.savefig('ρ_formula.tif',dpi=100)

结果如下：