我正在求解耦合方程的非线性系统，并计算了离散系统的雅可比行列式。结果真的很复杂，下面是（仅！）a 的前3 列矩阵，
（之所以会出现复杂性，部分原因是数值方案需要指数拟合才能实现稳定性。）
我对使用雅可比行列式实现数字代码有一个相当普遍的问题。
我可以继续在代码中实现这个矩阵。但是我的直觉告诉我，由于绝对的复杂性和不可避免的引入错误，我预计会进行几天（可能是几周！）繁琐的调试。在数字代码中如何应对这样的复杂性，这似乎是不可避免的？！您是否使用符号包的自动代码生成（然后手动调整代码）？
首先，我计划使用有限差分近似来调试解析雅可比行列式，我应该注意任何陷阱吗？您如何处理代码中的类似问题？
更新
一个词：模块化。
您的雅可比行列式中有很多重复的表达式可以写成它们自己的函数。您没有理由多次编写相同的操作，这将使调试更容易；如果你只写一次，那么只有一个地方出错（理论上）。
模块化代码也将使测试更容易；您可以为雅可比矩阵的每个组件编写测试，而不是尝试测试整个矩阵。例如，如果您以模块化方式编写函数 am()，您可以轻松地为它编写健全性测试，检查您是否正确区分它等等。
另一个建议是查看用于组装雅可比行列式的自动微分库。不能保证它们没有错误，但可能会比编写自己的调试/错误更少。这里有一些你可能想看的：
抱歉，刚刚看到您使用的是 python。ScientificPython 支持 AD。
有几种策略需要考虑：
使用 CAS 以符号形式查找导数，然后导出用于计算导数的代码。
使用自动微分 (AD) 工具生成代码，从代码计算导数以计算函数。
使用有限差分逼近来逼近雅可比。
自动微分可以生成更有效的代码来计算整个雅可比矩阵，然后使用符号计算为矩阵中的每个条目生成公式。有限差分是仔细检查衍生品的好方法。
让我在这里提几句谨慎的话，以一个故事开头。很久以前，我刚开始的时候和一个同事一起工作。他有一个优化问题要解决，目标相当混乱。他的解决方案是生成分析导数以进行优化。
我看到的问题是这些衍生品很讨厌。使用 Macsyma 生成，转换为 fortran 代码，它们每个都有几十个延续语句。事实上，Fortran 编译器对此感到不安，因为它超过了连续语句的最大数量。虽然我们找到了一个允许我们解决该问题的标志，但还有其他问题。
在 CA 系统通常生成的长表达式中，存在大量减法抵消的风险。计算大量大数，却发现它们都相互抵消以产生一个小数。
通常，解析生成的导数实际上比使用有限差分数值生成的导数更昂贵。n 个变量的梯度可能需要 n 倍于评估目标函数的成本。（您可能可以节省一些时间，因为许多术语可以在各种衍生品中重复使用，但这也将迫使您进行仔细的手工编码，而不是使用计算机生成的表达式。并且任何时候您编写讨厌的数学表达式，错误的概率不是微不足道的。确保验证这些导数的准确性。）
我的故事的重点是这些 CA 生成的表达式有它们自己的问题。有趣的是，我的同事实际上为问题的复杂性感到自豪，他显然正在解决一个非常困难的问题，因为代数是如此讨厌。我认为他没有考虑的是该代数是否真的在计算正确的东西，它是否如此准确，并且它是否如此有效。
如果我当时是这个项目的资深人士，我会读他的暴动行为。他的骄傲使他使用了一个可能过于复杂的解决方案，甚至没有检查基于有限差分的梯度是否足够。我敢打赌，我们可能花了一个人周的时间来运行这个优化。至少，我会建议他仔细测试产生的梯度。准确吗？与有限差分导数相比，它有多准确？事实上，今天有一些工具也可以在他们的导数预测中返回一个误差估计。这对于我用 MATLAB 编写的自适应微分代码（派生）当然是正确的。
测试代码。验证导数。
但在您执行任何此操作之前，请考虑是否可以选择其他更好的优化方案。例如，如果您正在进行指数拟合，那么您很有可能可以使用分区非线性最小二乘法（有时称为可分离最小二乘法。我认为这是 Seber 和 Wild 在他们的书中使用的术语。）这个想法是将参数集分解为本质上线性和本质上非线性的集合。使用仅适用于非线性参数的优化。鉴于这些参数是“已知的”，则可以使用简单的线性最小二乘法估计本质上的线性参数。该方案将减少优化中的参数空间。它使问题更加稳健，因为您不需要找到线性参数的起始值。它降低了搜索空间的维度，从而使问题运行得更快。我再次提供了用于此目的的工具，但仅在 MATLAB 中。
如果您确实使用分析导数，请将它们编码以重用术语。这可以节省大量时间，并且实际上可以减少错误，节省您自己的时间。但是然后检查这些数字！
这是我们在一个代码中使用 Sacado 使用自动微分的示例：http: //www.dealii.org/developer/doxygen/deal.II/step_33.html