基于Boost库的牛顿迭代求解：Jacobian计算与系统搭建优化问询

给定方程组：

$$x^2 + y^2 = 1$$
$$x - y = 0$$

希望通过牛顿迭代法求解，迭代公式为 J(x) * delta_x = -f(x)，其中J(x)是系统的雅可比矩阵。现有实现代码中方程组与系统搭建方式较为繁琐，以下是优化后的实现思路与代码：

现有代码的问题 #

原代码将每个方程单独定义为函数，导致方程逻辑分散，雅可比矩阵计算时需要重复调用每个方程，代码冗余且扩展性差。

优化方案与实现 #

核心优化思路 #

• 统一方程组的定义，让函数值计算与雅可比矩阵计算复用同一套方程逻辑
• 简化雅可比矩阵的生成流程，借助自动微分工具减少手动操作
• 封装通用迭代逻辑，提高代码复用性

优化后的完整代码 #

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/math/differentiation/autodiff.hpp>
#include <functional>
#include <tuple>

using namespace boost::numeric::ublas;
using namespace boost::math::differentiation;

// 统一定义所有方程组，返回包含各方程结果的tuple
template<typename T>
auto system_equations(const T& x, const T& y) {
    return std::make_tuple(
        x*x + y*y - 1,  // 对应方程x² + y² = 1
        x - y            // 对应方程x - y = 0
    );
}

// 计算函数值向量f(x)
vector<double> compute_f(const vector<double>& var) {
    vector<double> f(2, 0);
    auto [f0_val, f1_val] = system_equations(var[0], var[1]);
    f[0] = f0_val;
    f[1] = f1_val;
    return f;
}

// 计算雅可比矩阵J(x)
matrix<double> compute_J(const vector<double>& var) {
    matrix<double> J(2, 2, 0);
    // 创建支持1阶偏导的自动微分变量元组
    auto vars = make_ftuple<double, 1, 1>(var[0], var[1]);
    const auto& x = std::get<0>(vars);
    const auto& y = std::get<1>(vars);

    auto [f0, f1] = system_equations(x, y);
    // 提取偏导数填充雅可比矩阵
    J(0, 0) = f0.derivative(1, 0);  // ∂f0/∂x
    J(0, 1) = f0.derivative(0, 1);  // ∂f0/∂y
    J(1, 0) = f1.derivative(1, 0);  // ∂f1/∂x
    J(1, 1) = f1.derivative(0, 1);  // ∂f1/∂y

    return J;
}

// 封装通用牛顿迭代逻辑，支持N维变量的方程组求解
template<std::size_t N>
vector<double> newton_iteration(
    std::function<vector<double>(const vector<double>&)> f_func,
    std::function<matrix<double>(const vector<double>&)> j_func,
    vector<double> initial_guess,
    double eps = 1e-6,
    int max_iter = 1000
) {
    vector<double> x = std::move(initial_guess);
    int iter = 0;

    while (iter < max_iter) {
        matrix<double> Jx = j_func(x);
        vector<double> fx = f_func(x);
        vector<double> delta_x = -fx;

        // LU分解求解线性系统 Jx * delta_x = -fx
        permutation_matrix<std::size_t> pm(Jx.size1());
        lu_factorize(Jx, pm);
        lu_substitute(Jx, pm, delta_x);

        x += delta_x;
        iter++;

        if (norm_2(delta_x) < eps) break;

        if (iter % 10 == 0) {
            std::cout << "[" << iter << "] Solution: " << x << std::endl;
        }
    }

    std::cout << "Final answer:" << std::endl;
    std::cout << "[" << iter << "] Solution: " << x << std::endl;
    return x;
}

int main() {
    vector<double> initial_guess(2);
    initial_guess[0] = 0.5;
    initial_guess[1] = 0.5;

    newton_iteration<2>(compute_f, compute_J, initial_guess);
    return 0;
}

优化点说明 #

1. 统一方程定义：所有方程逻辑集中在system_equations函数中，compute_f和compute_J直接复用该函数，避免重复编写方程，降低维护成本。
2. 简化雅可比计算：通过C++17结构化绑定直接获取自动微分后的方程结果，无需逐个调用独立的方程函数，代码更紧凑。
3. 通用迭代封装：将牛顿迭代的核心逻辑封装为模板函数，后续扩展到更多变量或不同方程组时，只需修改方程定义部分，迭代逻辑无需改动。
4. 可读性提升：函数命名清晰，代码分层明确，逻辑流程一目了然，便于后续调试和扩展。

内容的提问来源于stack exchange，提问作者Фархад Оруджов