如何用Python求解含多个未知向量的矩阵方程组？

首先通过代数消元减少未知变量数量，再把剩余未知向量展平成一维数组，适配scipy.optimize.fsolve的输入要求。

步骤1：代数消元简化方程组 #

从第一个方程直接解出X3，减少一个未知量：

X2 - C1*X3 = X1 → X3 = (X2 - X1)/C1 （前提是C1≠0，若C1=0则X2=X1，直接代入第二个方程即可）

将X3代入第二个方程，结合第三个线性方程后，未知量只剩X2（N维）和X4（N维），总共有2N个未知参数。

步骤2：适配fsolve的一维输入要求 #

fsolve仅接受一维数组作为输入和输出，因此需要：

• 把X2和X4拼接成一个长度为2N的一维数组（例如先存放X2的所有元素，再存放X4的所有元素）
• 在目标函数内部，将一维数组重新拆分为X2和X4，计算各方程的残差，最后把残差拼接成一维数组返回。

步骤3：编写目标函数和求解代码 #

以下是具体的Python代码示例：

import numpy as np
from scipy.optimize import fsolve
from scipy.linalg import solve

# 假设已知参数已定义：C1, k, C2, C3, X1, b, A，N为向量维度
N = len(X1)

# 计算X4的初始猜测（解线性方程A@X4=C3*b）
X4_init = solve(A, C3 * b)
# X2的初始猜测可设为X1加小扰动，或直接用X1
X2_init = X1 + 0.1 * np.ones(N)

# 将初始猜测展平为一维数组
x0 = np.concatenate([X2_init, X4_init])

def objective(x):
    # 拆分一维数组为X2和X4
    X2 = x[:N]
    X4 = x[N:]
    
    # 计算方程2的残差：(X2 - X1)/C1 - k*(X4 - C2)*sqrt(X2*(1-X2))
    res2 = (X2 - X1)/C1 - k * (X4 - C2) * np.sqrt(X2 * (1 - X2))
    
    # 计算方程3的残差：A@X4 - C3*b
    res3 = A @ X4 - C3 * b
    
    # 拼接残差并展平为一维数组返回
    return np.concatenate([res2, res3])

# 调用fsolve求解
solution = fsolve(objective, x0)

# 拆分得到最终解
X2_sol = solution[:N]
X4_sol = solution[N:]
X3_sol = (X2_sol - X1)/C1

注意事项 #

• 若C1=0，直接将X2=X1代入第二个方程，此时X4可通过解线性方程A@X4=C3*b直接得到，再代入计算X3即可，无需使用fsolve。
• 初始猜测值对收敛性影响很大，需根据问题的物理背景调整X2和X4的初始值，避免离真实解过远。
• 需保证X2*(1-X2)非负（即X2每个元素在[0,1]区间内），否则开平方会得到NaN导致求解失败。若X2必须在该区间，可在目标函数中加入约束，或改用带约束的优化器（如scipy.optimize.root的特定方法、scipy.optimize.minimize结合残差平方和）。

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