Python求解非线性PDE隐式格式遇IndexError及残差数组构造错误

越界报错的直接成因 #

• 报错来源为残差数组f的维度设计完全偏离牛顿迭代的求解逻辑。使用隐式格式求解当前时间层的非线性方程时，残差是仅对应空间求解节点的1维向量，原代码错误将f初始化为(M+1, M+1)的二维数组，还把时间循环变量q作为f的第二维索引；当q取值达到31时，轴1的合法索引范围为0~30，直接触发index out of bounds错误。
• 原代码将残差与全时空数组强绑定，本质是混淆了“当前时间层求解的临时变量”和“全时空存储的解数组U”的维度边界。

求解结果偏差的核心逻辑问题 #

• 残差赋值逻辑完全错位：没有按照预期构造仅含空间维度的残差结构，正确的残差构造逻辑应为：
  • f[0] = 0：左边界节点值已知，残差为0
  • f[1:M]：所有内部求解节点的残差逐点计算
  • f[M] = 0：右边界节点值已知，残差为0
    原调整后的代码用循环给f赋值时，直接把f覆盖成单个标量值，完全无法形成正确的残差向量。
• 雅可比矩阵构造错误：原代码在10次牛顿迭代的外层循环里，每次仅取单个索引p点的导数数值，用这个常数值填充整个31*31的雅可比矩阵，没有针对每个内部节点逐点计算残差对邻点的偏导数，导致线性修正项计算完全失真。
• 边界条件应用错误：原代码错误将内部节点U[1, q+1]赋值为0，本来x=0（索引0）和x=1（索引M）才是固定为0的Dirichlet边界，人为把第一个内部节点的解强制置0，直接缩小了求解域，自然和显式方法的参考结果出现明显偏差。
• 雅可比公式存在笔误：原代码写雅可比对角项时，将-4*r*U[p, q+1]*U[p+1, q+1]误写为-4*r*U[p+1, q+1]*U[p+1, q+1]，进一步放大了导数计算的误差。

修正后的核心求解逻辑 #

关键调整点均围绕“当前时间层的牛顿迭代是对空间内部节点的1维求解”展开，核心代码片段如下：

for q in range(N):
    # 用上一时刻的解作为当前时间层的迭代初值
    U[:, q+1] = U[:, q]
    # 初始化边界条件
    U[0, q+1] = 0
    U[-1, q+1] = 0

    # 最多执行10次牛顿迭代
    for _ in range(10):
        # 初始化1维残差向量，覆盖所有空间节点
        f = np.zeros(M+1)
        # 逐内部节点计算残差
        for i in range(1, M):
            f[i] = (
                U[i, q+1]
                - r*(U[i+1, q+1] + U[i-1, q+1])
                + 2*r*U[i, q+1]
                - r*U[i, q+1]*U[i+1, q+1]**2
                + r*U[i, q+1]**2 * U[i+1, q+1]
                - r*U[i, q+1]**3
            )

        # 逐点计算雅可比矩阵的三类对角项
        lowdiag = np.zeros(M+1)
        diag = np.zeros(M+1)
        updiag = np.zeros(M+1)
        for i in range(1, M):
            lowdiag[i] = -r  # 残差对左邻点U[i-1]的偏导
            diag[i] = (
                1 + 2*r
                + 3*r*U[i, q+1]**2
                - 4*r*U[i, q+1]*U[i+1, q+1]
                - 2*r*U[i, q+1]**2
            )  # 残差对自身节点U[i]的偏导
            updiag[i] = (
                -r
                - 2*r*U[i, q+1]*U[i+1, q+1]
                + 2*r*U[i, q+1]**2
            )  # 残差对右邻点U[i+1]的偏导

        # 构造仅针对内部节点的雅可比矩阵（大小为(M-1, M-1)）
        J_data = np.array([diag[1:M], lowdiag[1:M], updiag[1:M]])
        J = spdiags(J_data, [0, -1, 1], M-1, M-1).toarray()

        # 求解线性修正量 delta，满足 J*delta = -f
        delta = np.linalg.solve(J, -f[1:M])
        # 更新内部节点的解
        U[1:M, q+1] += delta

        # 重新固定边界条件
        U[0, q+1] = 0
        U[-1, q+1] = 0

修正后可彻底消除索引越界问题，残差与雅可比的构造逻辑和数值格式完全匹配，求解结果可与显式方法的参考结果保持一致。

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