python - 生成方程组

Question

我正在尝试创建一个可以生成方程组的函数，该方程组可以在单独的程序中求解。这些方程式是从同位素衰变树生成的，但为简单起见，我有以下树:

所以这可以制成 2 个可能的衰变链:

[(A,0,1,5), (B,1,.4,4), (C,0,.4,0)]

[(A,0,1,5), (B,1,.6,6), (C,0,.6,0)]

格式是(物种、数量、衰变概率、半衰期)。我正在尝试创建一个函数，该函数会自动为一棵腐烂树创建一个方程组，这可能比这更复杂。规则对于任何树都是相同的:

对于某些具有 parent Y_1、Y_2、...、Y_n 的物种 X:

X_final = 每个亲本物种的总和(Y_n 的衰变概率 -> X * Y_n 的数量/Y_n 的半衰期)- X 的数量/X 的半衰期，这可能表示为:

并且链中的每个物种都有自己的方程式，稍后需要求解。因此，为此，我需要以下方程组:

A_f = - A_i/5
B1_f = .4 * A_i/5 - B1_i / 4
B2_f = .6 * A_i/5 - Β2_i / 6
C = B1_i / 4 + B2_i / 6

此外，如果半衰期为 0，则表示它是稳定的。目前，我正在通过制作字符串字典来生成方程组，但我认为有更好的方法来做到这一点。我计划在使用字符串创建系统后稍后将字符串转换为变量。这是我的代码:

A = 'A'
B = 'B'
C = 'C'
D = 'D'

chain1 = [(A,0,1,5),(B,1,.4,4),(C,0,.4,0),(D,0,.4,0)]
chain2 = [(A,0,1,5),(B,2,.6,6),(C,0,.6,0),(D,0,.6,0)]
master_chain = [chain1, chain2]

def equation_builder(master_chain):
    master_equations = {}
    m = 0
    for chain in master_chain:
        n = 0
        for item in chain:
            if item == chain[0]:
                equation = {str(item[0]) + str(item[1]) + 'f' :\
                '-' +  str(item[0]) + str(item[1])  + 'i/' +  str(item[3])}
                master_equations.update(equation)
            elif str(item[0])+str(item[1])+'f' not in master_equations:
                equation = {str(item[0]) + str(item[1]) + 'f' :\
                str(item[2]/chain[n-1][2])+str(chain[n-1][0]) + 
                str(chain[n-1][1])+'i/' + str(chain[n-1][3])+\
                '-'+str(item[0])+str(item[1])+'i/'+str(item[3])}
                master_equations.update(equation)
            elif str(item[0])+str(item[1])+'f' in master_equations \
            and master_chain[m-1][n-1] != master_chain[m][n-1]:
                old_equation = master_equations[str(item[0])+str(item[1])+'f']
                new_equation = old_equation + '+' +\
                str(item[2]/chain[n-1][2])+str(chain[n-1][0]) +\
                str(chain[n-1][1])+'i/' + str(chain[n-1][3])
                equation = {str(item[0])+str(item[1])+'f' : new_equation}
                master_equations.update(equation)
            n += 1
        m += 1

    return master_equations

if __name__ == '__main__':
    print equation_builder(master_chain)

Answer 1

使用SymPy . SymPy 是一个符号计算工具箱，非常适合这种用例。您可以使用 A = sympy.Symbol("A") 创建符号，然后像使用任何变量一样在表达式中使用 A。例如，如果 A 和 B 是符号，那么如果你写 C=A*exp(B)，print C 将输出 A*exp(B)。使用表达式的 args 属性，您还可以访问任何表达式的语法树表示，如果您想进一步处理方程式，这可能会很有用。

这是一个使用你的图表的例子(我不太明白你是如何得出结果的，所以这可能需要一些调整，但它应该足以理解这个想法):

import sympy as sp

A, B1, B2, C = sp.symbols("A, B1, B2, C")

chain1 = [(A,0,1,5),(B1,1,.4,4),(C,0,0.4,0)]
chain2 = [(A,0,1,5),(B2,2,.6,6),(C,0,0.6,0)]
master_chain = [chain1, chain2]

finals = {}
for subchain in master_chain:
    for i, (species, number, decay_prob, half_life) in enumerate(subchain):
        input_species = sp.Symbol(str(species) + "_i")
        if species not in finals:
            finals[species] = -input_species / half_life if half_life else 0
        if i < len(subchain) - 1:
            (other_species, other_number, other_decay_prob, other_half_life) = subchain[i+1]
            if other_species not in finals:
                finals[other_species] = -sp.Symbol(str(other_species) + "_i") / other_half_life if other_half_life else 0
            finals[other_species] += input_species * decay_prob / half_life

print finals

输出是

{C: 0.1*B1_i + 0.1*B2_i, B2: A_i/5 - B2_i/6, A: -A_i/5, B1: A_i/5 - B1_i/4}

请注意 Symbol("x") == Symbol("x")，例如，符号由其字符串表示形式标识，因此您可以在每次需要时安全地重新创建符号。