VBA中Newton-Raphson法计算n元酸滴定曲线不收敛问题

代码中的核心错误分析 #

• B数组赋值错误：输入解离常数B值时，代码错误地将所有值存入B(n)而非对应索引B(j)，导致后续计算使用的解离常数完全失效。
• 变量未初始化：S1、S2、S3、S4及导数相关的Sd1-Sd4未在循环前初始化为0，VBA中未初始化变量为Empty，参与运算后会产生无效值。
• S系列计算位置错误：原代码仅在滴定前计算一次S值，但每次滴定剂加入后，H⁺浓度、酸浓度（Caci）和OH⁻浓度（M）均会改变，S值需随当前H值动态计算，而非固定初始值。
• 初始值设置不合理：初始Hini设为12（强碱性），对于未滴定的酸性溶液，初始值与真实值偏差极大，直接导致牛顿法发散。应根据滴定阶段设置合理初始值，且每次滴定迭代时用上一次的收敛值作为新初始值。
• 迭代循环逻辑错误：Y的初始值仅在全局设置一次，后续滴定循环中未重置，导致只有第一次滴定会进入迭代，后续直接跳过。需在每次滴定步骤中重置迭代条件。

修正后的VBA代码 #

Sub Calculo_Curva_Titulacion()
    '变量定义
    Dim n As Integer
    Dim j As Integer, t As Integer, el As Integer
    Dim Kw As Double, Mini As Double, M As Double
    Dim VM As Double, Cacini As Double, Caci As Double
    Dim Vali As Double, Vf As Double, VMtit As Double
    
    '用户输入参数
    n = InputBox("输入酸的质子数（n元酸）", "质子数", 1)
    Kw = 10 ^ -14
    
    ReDim B(n) As Double
    t = 1
    For j = 1 To n
        '修正：将值存入对应索引B(j)
        B(j) = InputBox("输入第" & j & "个解离常数的B值", "B" & j, 2)
        Cells(t, 1) = j
        Cells(t, 2) = B(j)
        t = t + 1
    Next j
    
    Cacini = InputBox("输入酸的初始浓度（M）", "酸浓度", 0.1)
    Vali = InputBox("输入酸溶液体积（mL）", "酸体积", 10)
    Mini = InputBox("输入NaOH浓度（M）", "NaOH浓度", 0.1)
    VMtit = InputBox("输入每次滴定加入的NaOH体积（mL）", "滴定步长", 0.5)
    
    '滴定终点体积（扩展10%避免提前终止）
    Vf = (n * Cacini * Vali) / Mini * 1.1
    el = 10 '结果起始行
    Cells(el, 1) = "NaOH体积(mL)"
    Cells(el, 2) = "H+浓度"
    Cells(el, 3) = "pH"
    el = el + 1
    
    Dim S1, S2, S3, S4 As Double
    Dim Sd1, Sd2, Sd3, Sd4 As Double
    Dim Hini As Double, Hnew As Double, Y As Double
    Dim Tol As Double
    Tol = 1e-10 '收敛阈值
    
    '初始H+浓度估算（未滴定的酸溶液）
    Hini = Sqr(Cacini * 10 ^ -B(1)) '一元酸近似，多元酸可调整
    Hnew = Hini
    Y = Tol + 1 '确保首次进入迭代
    
    For VM = 0 To Vf Step VMtit
        Caci = (Cacini * Vali) / (Vali + VM)
        M = (Mini * VM) / (Vali + VM)
        
        '重置迭代条件
        Y = Tol + 1
        '用上一次收敛值作为本次初始值，加快收敛
        If VM > 0 Then Hini = Hnew
        
        '牛顿-拉夫逊迭代
        While Y > Tol
            '每次迭代重新计算S系列值（基于当前Hini）
            S1 = 0: S2 = 0: S3 = 0: S4 = 0
            For j = 1 To n
                S1 = S1 + B(j) * (Hini ^ (j + 2))
                S2 = S2 + B(j) * (Hini ^ (j + 1))
                S3 = S3 + B(j) * (Hini ^ j)
            Next j
            For j = 1 To n - 1
                S4 = S4 + (n - j) * B(j) * (Hini ^ (j + 1))
            Next j
            
            '计算导数项Sd系列
            Sd1 = 0: Sd2 = 0: Sd3 = 0: Sd4 = 0
            For j = 1 To n
                Sd1 = Sd1 + (j + 2) * B(j) * (Hini ^ (j + 1))
                Sd2 = Sd2 + (j + 1) * B(j) * (Hini ^ j)
                Sd3 = Sd3 + j * B(j) * (Hini ^ (j - 1))
            Next j
            For j = 1 To n - 1
                Sd4 = Sd4 + (n - j) * (j + 1) * B(j) * (Hini ^ j)
            Next j
            
            '迭代计算新的H+浓度
            Hnew = Hini - (Polinomio(S1, S2, S3, S4, Kw, Hini, n, Caci, M) / DevPolinomio(Sd1, Sd2, Sd3, Sd4, Kw, Hini, n, Caci, M))
            '防止出现负值或不合理值
            If Hnew <= 0 Then Hnew = 1e-14
            Y = Abs(Hnew - Hini)
            Hini = Hnew
        Wend
        
        '写入结果
        Cells(el, 1) = VM
        Cells(el, 2) = Hnew
        Cells(el, 3) = -Log(Hnew) / Log(10)
        el = el + 1
    Next VM
End Sub

Function Polinomio(S1, S2, S3, S4, Kw, H, n, Caci, M) As Double
    '质子平衡方程
    Polinomio = S1 + M * S2 - Caci * S4 - Kw * S3 + H ^ 2 + (M - n * Caci) * H - Kw
End Function

Function DevPolinomio(Sd1, Sd2, Sd3, Sd4, Kw, H, n, Caci, M) As Double
    '质子平衡方程的导数
    DevPolinomio = Sd1 + M * Sd2 - Caci * Sd4 - Kw * Sd3 + 2 * H + (M - n * Caci)
End Function

关键修正说明 #

1. B数组赋值修正：将B(n)改为B(j)，确保每个解离常数存入对应索引。
2. 变量初始化：每次计算S和Sd前，都将其重置为0，避免累计错误值。
3. 动态计算S系列：将S和Sd的计算放入牛顿迭代循环内，确保每次迭代都基于当前H值重新计算。
4. 合理初始值：初始Hini基于酸的近似H+浓度估算，后续滴定用上一次收敛值作为初始值，大幅提升收敛速度和稳定性。
5. 迭代逻辑重置：在每个滴定步骤中重置Y为大于阈值的值，确保每次滴定都会进入迭代计算。
6. 异常值处理：加入If Hnew <=0 Then Hnew=1e-14，避免出现无效的H+浓度值。

内容的提问来源于stack exchange，提问作者Andrés Felipe Perdomo P.