Scipy 根方程求解

Question

我在求解方程组时遇到问题。要理解代码，

术语：(1. - Ni_oct + delta_Ni_mod) 曾经是一个变量。但我必须实现它才能删除此处未看到的一个方程。该方程是在变量 Q 中看到的。Q 应等于循环中的可迭代 j。目标是找到 diff 满足容差的 delta_Ni_mod 。但由于某种原因，代码找不到任何解决方案并且无限运行。增加容差不是一个选择，因为这样我得到的解决方案在逻辑上不满足适当的差异。顺便说一句，我应该提到可迭代 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```

Answer 1

重写程序的第一部分，使其看起来像这样

import numpy as np
from scipy.optimize import root

from scipy.constants import (
    e,    # atomic unit of charge (in C)
    N_A,  # Avogadro constant
    R,    # molar gas constant
)

delta_Ni_mod_max = -0.002
delta_Ni_mod_min = 0
delta_Ni_mod = (delta_Ni_mod_min + delta_Ni_mod_max)/2


def equations(vars: np.ndarray, l: float, m: float) -> tuple[float, ...]:
    Ni_oct, V_oct, O_o, V_o, Al_tet, Al_oct, V_tet = vars

    # Oktaederplätze / octahedral sites
    #    + Ni_oct              - Ni_oct
    eq1 =        + V_oct + .999        + delta_Ni_mod + Al_oct - 2

    # Sauerstoffplätze / oxygen sites
    eq2 = O_o + V_o - 4

    # Ladungsbilanz / charge balance
    eq3 = -Ni_oct - 3*V_oct + 2*V_o + Al_tet - 2*V_tet

    # Tetraederplätze / tetrahedral sites
    eq7 = Al_tet + V_tet - 1

    # Massebilanz Al / mass balance, aluminium
    eq8 = Al_tet + Al_oct - 2

    # # Kp_def5_neutr Site exchange
    eq10 = Al_oct * V_tet/(V_oct * Al_tet) - l

    # K Bildung Sauerstoffvakanzen / formation oxygen vacancies
    eq11 = V_o*Ni_oct**2/(O_o * (1 - Ni_oct + delta_Ni_mod)**2) - m

    return (eq1, eq2, eq3, eq7, eq8, eq10, eq11)


def main() -> None:
    # tolerance = 1e-21

    diff_list = []
    delta_Ni_mod_values = []
    resultlist_def1_3_4_5_6 = []

    dH_3 = 5.02*e
    dH_5_neutr = 0.708*e
    dH_6_neutr = 6.151*e

    T_values = np.arange(0, 1300, 50)
    N_A_RT = N_A/(R * T_values)

    p_H2 = np.logspace(-6, 1, num=100)
    K_def3 = np.exp(-dH_3 * N_A_RT)
    Kp_def3 = K_def3
    K_def5_neutr = np.exp(-dH_5_neutr * N_A_RT)
    Kp_def5_neutr = K_def5_neutr
    K_def6_neutr = np.exp(205.15/2/R) * np.exp(-dH_6_neutr * N_A_RT)
    Kp_def6_neutr = np.outer(K_def6_neutr, p_H2/10e-5)

    for j, l, Kp_def6_subarray in zip(Kp_def3, Kp_def5_neutr, Kp_def6_neutr):
        for m in Kp_def6_subarray:
            initial_guess = (
                0.999,   # Ni_oct
                0.0001,  # V_oct
                3.99,    # O_o
                0.0001,  # V_o
                1,       # Al_tet
                1,       # Al_oct
                0.001,   # V_tet
            )

删除内部循环并将其转换为附加方程。另外，将整个程序转换为使用对数刻度；否则你肯定会遇到稳定性和准确性问题。我将在稍后的编辑中演示如何执行此操作。