寻求适用于欠定非线性方程组的高效稳健Binodal相图求解方法

针对复刻Deviri and Safran (2021)中Fig3/4双节线相图的需求，核心是求解含4个未知量、3个方程的欠定非线性方程组。以下是基于numpy/scipy生态的高效稳健方案，适配任意可推导化学势与渗透压的自由能函数$F$：

核心思路 #

欠定系统的破局点是固定一个自由度，将问题转化为正定方程组求解，再通过**连续法（路径跟踪）**遍历参数空间——既避免暴力搜索的低效，也解决单一优化的收敛稳定性问题：

• 固定其中一相的某一组分浓度（比如相1的$\phi_a^{{(1)}$），将剩余3个未知量$\phi_b}{(1)}, \phi_a^{(2)}, \phi_b^{(2)}$作为待求解变量
• 从一个已知的可靠平衡点（比如临界点、稀释极限解）出发，逐步调整固定的浓度参数，以上一步的解作为初始值调用非线性求解器，确保收敛连续性

实现步骤 #

1. 定义自由能衍生函数 #

首先从自由能$F(\phi_a, \phi_b)$推导三个关键函数（优先用解析形式，数值求导为备选）：

• 组分a的化学势：$\mu_a = \frac{\partial F}{\partial \phi_a}$
• 组分b的化学势：$\mu_b = \frac{\partial F}{\partial \phi_b}$
• 渗透压：$\Pi = \phi_a \frac{\partial F}{\partial \phi_a} + \phi_b \frac{\partial F}{\partial \phi_b} - F$

示例代码（以Flory-Huggins型自由能为例，可替换为目标论文的自由能形式）：

import numpy as np
from scipy.optimize import root

def free_energy(phi_a, phi_b, chi_ab, chi_as, chi_bs):
    phi_s = 1 - phi_a - phi_b
    if phi_s <= 0 or phi_a <=0 or phi_b <=0:
        return np.inf
    term_a = phi_a * np.log(phi_a)
    term_b = phi_b * np.log(phi_b)
    term_s = phi_s * np.log(phi_s)
    interaction = chi_ab * phi_a * phi_b + chi_as * phi_a * phi_s + chi_bs * phi_b * phi_s
    return term_a + term_b + term_s + interaction

def mu_a(phi_a, phi_b, *params):
    # 解析导数优先，此处用数值求导示例
    h = 1e-8
    return (free_energy(phi_a+h, phi_b, *params) - free_energy(phi_a-h, phi_b, *params))/(2*h)

def mu_b(phi_a, phi_b, *params):
    h = 1e-8
    return (free_energy(phi_a, phi_b+h, *params) - free_energy(phi_a, phi_b-h, *params))/(2*h)

def pi(phi_a, phi_b, *params):
    return phi_a * mu_a(phi_a, phi_b, *params) + phi_b * mu_b(phi_a, phi_b, *params) - free_energy(phi_a, phi_b, *params)

2. 构建方程组与求解器 #

针对固定的$\phi_a^{(1)}$，定义残差函数，目标是让三个平衡方程的残差为0：

def residual(x, phi_a1, params):
    phi_b1, phi_a2, phi_b2 = x
    # 三个平衡条件的残差
    res1 = mu_a(phi_a1, phi_b1, *params) - mu_a(phi_a2, phi_b2, *params)
    res2 = mu_b(phi_a1, phi_b1, *params) - mu_b(phi_a2, phi_b2, *params)
    res3 = pi(phi_a1, phi_b1, *params) - pi(phi_a2, phi_b2, *params)
    return [res1, res2, res3]

3. 连续法遍历参数空间 #

从初始平衡点开始，逐步扫描$\phi_a^{(1)}$的取值范围，每次用前一次的解作为初始值，大幅提升收敛概率：

def solve_binodal(params, phi_a1_range, initial_guess):
    binodal_phase1 = []
    binodal_phase2 = []
    prev_sol = initial_guess  # 初始解：临界点或已知两相点

    for phi_a1 in np.linspace(*phi_a1_range, 100):
        sol = root(residual, prev_sol, args=(phi_a1, params), method='lm')
        if sol.success:
            phi_b1, phi_a2, phi_b2 = sol.x
            # 验证相组成的物理合理性
            if all([phi_b1 >=0, phi_a2 >=0, phi_b2 >=0, 
                   phi_a1+phi_b1 <=1, phi_a2+phi_b2 <=1]):
                binodal_phase1.append((phi_a1, phi_b1))
                binodal_phase2.append((phi_a2, phi_b2))
                prev_sol = sol.x  # 更新初始值为当前有效解
        else:
            # 收敛失败时可尝试缩小步长或微调初始值
            continue
    return np.array(binodal_phase1), np.array(binodal_phase2)

4. 调用与可视化 #

# 替换为Deviri and Safran 2021的参数
params = (0.5, 0.3, 0.2)  # chi_ab, chi_as, chi_bs
# 初始猜测：可从临界点或稀释极限解获取
initial_guess = [0.1, 0.2, 0.3]
# 设定phi_a1的扫描范围
phi_a1_range = (0.01, 0.4)

phase1, phase2 = solve_binodal(params, phi_a1_range, initial_guess)

# 绘制双节线
import matplotlib.pyplot as plt
plt.scatter(phase1[:,0], phase1[:,1], label='Phase 1')
plt.scatter(phase2[:,0], phase2[:,1], label='Phase 2')
plt.xlabel(r'$\phi_a$')
plt.ylabel(r'$\phi_b$')
plt.legend()
plt.show()

关键优化点 #

• 优先使用解析导数：数值求导易引入误差，解析形式能大幅提升求解稳定性与速度
• 初始值精准选择：若未知初始平衡点，可先求解临界点（通过二阶偏导为0的条件），从临界点出发扫描能覆盖完整双节线
• 求解器适配：scipy的root方法中，lm（Levenberg-Marquardt）适合非线性最小二乘问题；遇刚性问题可尝试hybr方法
• 动态步长调整：某点收敛失败时，可缩小该区域的扫描步长，用插值后的初始值重新尝试

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