Scipy方程组求解遇死循环：无法找到符合精度的delta_Ni_mod

问题详情 #

• 原表达式(0.999 - Ni_oct + delta_Ni_mod)曾作为独立变量，为移除未展示的方程将其直接代入现有方程组
• 需求：找到使diff = j - Q满足精度要求（≤1e-20）的delta_Ni_mod，其中j是循环中的迭代量，Q由方程组的解计算得出
• 当前问题：
  • 代码陷入无限循环，始终无法找到符合精度的解
  • j的取值范围极小（10^-170 到 10^-21），强行提高计算精度会得到不符合物理逻辑的解

测试代码 #

import numpy as np
from scipy.optimize import root


dH_3 = 5.02*1.602e-19
N_A=6.022e23
R = 8.314 
T_values = np.arange(0,1300, 50)
dH_5_neutr = 0.708*1.602e-19
p_H2 = np.logspace(-6,1, num=100)
dH_6_neutr = 6.151*1.602e-19

p_H2 = np.logspace(-6,1, num=100)
K_def3 = np.exp(-(dH_3*N_A)/(R*T_values))
Kp_def3=K_def3
K_def5_neutr = np.exp(-(dH_5_neutr*N_A)/(R*T_values))
Kp_def5_neutr=K_def5_neutr
K_def6_neutr = np.exp((205.15/2)/R)*np.exp(-(dH_6_neutr*N_A)/(R*T_values))
Kp_def6_neutr = [[K_def6_neutr[i] * (p_H2[j]/10e-5) for j in range(len(p_H2))] for i in range(len(K_def6_neutr))]



# Initialize variables
delta_Ni_mod_max = np.float64(-0.002)
delta_Ni_mod_min = np.float64(0.0)
tolerance = 1e-20
diff_list = []
delta_Ni_mod_values = []
resultlist_def1_3_4_5_6 = []
Q_list = []
delta_Ni_mod = (delta_Ni_mod_min + delta_Ni_mod_max) / 2
for j, l, Kp_def6_subarray in zip(Kp_def3, Kp_def5_neutr, Kp_def6_neutr):
    for m in Kp_def6_subarray:
        def equations(vars):
            Ni_oct, V_oct, O_o, V_o, Al_tet, Al_oct, V_tet = vars

            eq1 = Ni_oct + V_oct + (0.999 - Ni_oct + delta_Ni_mod) + Al_oct - 2.0  # Oktaederplätze
            eq2 = O_o + V_o - 4.0  # Sauerstoffplätze
            eq3 = -1 * Ni_oct - 3 * V_oct + 2 * V_o + Al_tet - 2 * V_tet  # Ladungsbilanz
            eq7 = Al_tet + V_tet - 1  # Tetraederplätze
            eq8 = Al_tet + Al_oct - 2.0  # Massebilanz Al
            eq10 = (Al_oct * V_tet) / (V_oct * Al_tet) - l  # Kp_def5_neutr Site exchange
            eq11 = ((V_o * (Ni_oct ** 2)) / (O_o * ((0.999 - Ni_oct + delta_Ni_mod) ** 2))) - m  # K Bildung Sauerstoffvakanzen

            return [eq1, eq2, eq3, eq7, eq8, eq10, eq11]

        initial_guess = [0.999,  # Ni_oct # 0.999
                         0.0001,  # V_oct
                         3.999,  # O_o
                         0.0001,  # V_o
                         1.,  # Al_tet
                         1.,  # Al_oct
                         0.001]  # V_tet
        #delta_Ni_mod = (delta_Ni_mod_min + delta_Ni_mod_max) / 2

        while True:
            
            sol = root(equations, initial_guess, method='lm', tol=1e-40)
            Ni_oct, V_oct, O_o, V_o, Al_tet, Al_oct, V_tet = sol.x

            Q = (V_oct * (0.999 - Ni_oct + delta_Ni_mod) ** 2) / (Ni_oct ** 3)
            diff = j - Q

            if abs(diff) <= tolerance:
                diff_list.append(diff)
                delta_Ni_mod_values.append(delta_Ni_mod)
                resultlist_def1_3_4_5_6.append(sol.x)
                Q_list.append(Q)
                break

            if diff > 0:
                # Q muss positiver gemacht werden, also delta_Ni muss richtung negative Zahl gebtacht werden

                delta_Ni_mod_min = delta_Ni_mod
            else:
                # Hier gilt K < Q, also muss Q negativer gemacht werden
                # Q muss negativer gemacht werden, also delta_Ni Ri 0 gebracht werden
                delta_Ni_mod_max = delta_Ni_mod

            delta_Ni_mod = (delta_Ni_mod_min + delta_Ni_mod_max) / 2

内容的提问来源于stack exchange，提问作者NND