我在求解方程组时遇到问题。要理解代码,
术语:(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```
重写程序的第一部分,使其看起来像这样
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
)
删除内部循环并将其转换为附加方程。另外,将整个程序转换为使用对数刻度;否则你肯定会遇到稳定性和准确性问题。我将在稍后的编辑中演示如何执行此操作。