如何获取Scilab的balanc()函数实现逻辑以在Maxima中复现该功能

核心逻辑说明 #

Scilab的balanc()底层直接调用LAPACK标准的xGEBAL算法，无需查找Scilab源码，按照该公开算法逻辑实现即可完全对齐输出结果，算法目标是通过相似变换将矩阵调整为行、列∞范数近似相等的形式，降低后续数值计算的误差。

分步实现步骤 #

• 第一步：置换矩阵预处理
  遍历输入矩阵A的所有行和列，筛选出满足以下条件的索引：该行除对角元外所有元素绝对值小于阈值（默认取机器精度的10倍），或该列除对角元外所有元素绝对值小于阈值。通过置换矩阵P将这类行列移动到矩阵边缘，分离出孤立的特征值，得到第一次变换后的矩阵P^T · A · P。
• 第二步：对角缩放迭代
  1. 初始化对角缩放矩阵D为单位矩阵，设置迭代容差为2（即行范数和对应列范数的比值不超过2），最大迭代次数设为10次（该算法收敛速度极快，通常5次以内即可收敛）
  2. 每次迭代中逐行计算当前矩阵的行∞范数r_i、列∞范数c_i
  3. 对每个索引i，寻找2的整数次幂s，使得调整后的范数比值s*r_i / (c_i/s)尽可能接近1，选择2的整数次幂是为了避免浮点舍入误差
  4. 若s≠1，则更新D[i,i] = D[i,i] * s，同时将当前矩阵的第i行乘以s、第i列除以s
  5. 所有索引遍历完成后，检查所有行列范数比值是否都在容差范围内，是则终止迭代，否则进入下一轮迭代
• 第三步：合并输出结果
  最终的相似变换矩阵X = P · D，平衡后的矩阵Ab = X^(-1) · A · X，和Scilabbalanc()的双返回值输出完全对应。

适配wxMaxima的注意事项 #

• wxMaxima自带的相似变换函数未采用2次幂缩放、固定迭代规则和LAPACK标准的置换顺序，是输出结果和Scilab不一致的核心原因，严格按照上述规则实现即可对齐
• 不要采用符号变量求解方程组的方式实现，该算法本身是数值迭代逻辑，符号场景下将缩放因子表示为2^k（k为整数）的形式即可保留精确性
• 实现tf2ss时，先完成传递函数到可控标准型状态空间的转换，再调用上述balanc实现处理状态矩阵，即可完全复现Scilab的tf2ss输出逻辑，后续的离散化、z域传递函数转换步骤也能直接对齐。

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