c++ Boost库,雅可比矩阵和系统设置

7xzttuei  于 11个月前  发布在  其他
关注(0)|答案(2)|浏览(129)

假设我有一个方程组：

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>

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

template<typename X, typename Y>
promote<X, Y> f0(const X &x, const Y &y) {
    return x * x + y + y - 1;
}

template<typename X, typename Y>
promote<X, Y> f1(const X &x, const Y &y) {
    return x - y;
}

vector<double> f(const vector<double> &var) {
    vector<double> f(2, 0);
    f[0] = f0(var[0], var[1]);
    f[1] = f1(var[0], var[1]);
    return f;
}

matrix<double> J(const vector<double> &var) {
    matrix<double> J(2, 2, 0);
    
    auto const variables = make_ftuple<double, 1, 1>(var[0], var[1]);
    
    auto const &x = std::get<0>(variables);
    auto const &y = std::get<1>(variables);
    
    J(0, 0) = f0(x, y).derivative(1, 0);
    J(0, 1) = f0(x, y).derivative(0, 1);
    J(1, 0) = f1(x, y).derivative(1, 0);
    J(1, 1) = f1(x, y).derivative(0, 1);
    
    return J;
}

int main() {
    vector<double> x(2, 0);
    x[0] = 0.5;
    x[1] = 0.5;
    double eps = 1.e-6;
    int iter = 0;
    
    vector<double> fx, delta_x;
    matrix<double> Jx;
    
    while (iter < 1.e3) {
        Jx = J(x);
        fx = f(x);
        delta_x = -fx;
        // Solve the linear system J(x) * delta_x = -f(x)
        
        permutation_matrix<std::size_t> pm(Jx.size1());
        lu_factorize(Jx, pm);
        lu_substitute(Jx, pm, delta_x);
        
        x += delta_x;
        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 0;
}

型

c++

2条答案

按热度按时间
ldfqzlk8

ldfqzlk81#

尽管你没有真正解释代码中“尴尬”的地方。
我看了一下，并能够提出一些优化和简化（例如使用C++17语言功能）。
注意，为了清晰（和安全），我使用命名空间别名：

namespace U = boost::numeric::ublas;
namespace D = boost::math::differentiation;

字符串
优化是为了避免向量/矩阵的无限数组存储，因为它们是固定大小的：

template <size_t N, typename T = double> using FixedVector = U::vector<T, std::array<T, N>>;
using Vec = FixedVector<2>;
using Mat = U::matrix<double, U::row_major, std::array<double, 4>>;

型
我选择了std::array，因为它很简单，并且对constexpr友好，并且（与U::bounded_array<>相反）允许构造初始化器列表。这将允许我们更自然的赋值：

Vec f(Vec const& var) {
    auto& [x, y] = var.data();
    return Vec(2, {f0(x, y), f1(x, y)});
}

Mat J(Vec const& var) {
    auto const& [x, y] = D::make_ftuple<double, 1, 1>(var[0], var[1]);
    return Mat{2, 2,
        {
            f0(x, y).derivative(1, 0),
            f0(x, y).derivative(0, 1),
            f1(x, y).derivative(1, 0),
            f1(x, y).derivative(0, 1),
        }};
}

型
另一部分是使用结构化绑定，而不是get<n>或[n]舞蹈。
从效率的Angular 来看，我将验证编译器是否消除了矩阵初始化器中的公共子表达式。
对于f0/f1函数，我让类型演绎来完成工作：

constexpr auto f0(auto const& x, auto const& y) { return x * x + y + y - 1; }
constexpr auto f1(auto const& x, auto const& y) { return x - y; }

型
主驱动程序没有太大变化，尽管我重新塑造了循环，使其更传统：

Live On Coliru

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

constexpr auto f0(auto const& x, auto const& y) { return x * x + y + y - 1; }
constexpr auto f1(auto const& x, auto const& y) { return x - y; }

namespace U = boost::numeric::ublas;
namespace D = boost::math::differentiation;

template <size_t N, typename T = double> using FixedVector = U::vector<T, std::array<T, N>>;
using Vec = FixedVector<2>;
using Mat = U::matrix<double, U::row_major, std::array<double, 4>>;

Vec f(Vec const& var) {
    auto& [x, y] = var.data();
    return Vec(2, {f0(x, y), f1(x, y)});
}

Mat J(Vec const& var) {
    auto const& [x, y] = D::make_ftuple<double, 1, 1>(var[0], var[1]);
    return Mat{2, 2,
        {
            f0(x, y).derivative(1, 0),
            f0(x, y).derivative(0, 1),
            f1(x, y).derivative(1, 0),
            f1(x, y).derivative(0, 1),
        }};
}

int main() {
    constexpr double eps  = 1.e-6;
    Vec x(2, {0.5, 0.5});

    U::vector<double> fx, delta_x;
    U::matrix<double> Jx;

    for (auto iteration = 0; iteration < 1'000; ++iteration) {
        Jx      = J(x);
        fx      = f(x);
        delta_x = -fx;
        // Solve the linear system J(x) * delta_x = -f(x)

        U::permutation_matrix<std::size_t> pm(Jx.size1());
        lu_factorize(Jx, pm);
        lu_substitute(Jx, pm, delta_x);
        
        x += delta_x;
        std::cout << "[" << iteration << "] Solution: " << x << std::endl;

        if (norm_2(delta_x) < eps)
            break;
    }
    std::cout << "Solution: " << x << std::endl;
}

型
印刷

[0] Solution: [2](0.416667,0.416667)
[1] Solution: [2](0.414216,0.414216)
[2] Solution: [2](0.414214,0.414214)
[3] Solution: [2](0.414214,0.414214)
Solution: [2](0.414214,0.414214)

型

摘要

我希望这些更改能使代码稍微快一点（看到没有I/O的微基准测试会很有趣）。
我不确定这些是否解决了你在问题中提到的一些“尴尬”。

赞(0）分享  回复(0）举报  11个月前
ltskdhd1

ltskdhd12#

我想你可以不用升压，自己卷。
方程采用泛函建立，雅可比矩阵采用有限差分法求解，但还有待改进。

#include <iostream>
#include <vector>
#include <functional>
#include <utility>
#include <cmath>
using namespace std;

using FUNC = function< double( vector<double> ) >;         // double-valued function with vector argument

bool solve( vector<FUNC> eqns, vector<double> &x );
vector<double> operator -( vector<double> a, const vector<double> &b );
vector<double> getFunction( vector<FUNC> eqns, const vector<double> &x );
vector< vector<double> > getJacobian( vector<FUNC> eqns, const vector<double> &x );
double getResiduals( const vector<double> &f );
bool GaussElimination( vector< vector<double> > A, vector<double> B , vector<double> &X );

//======================================================================

int main()
{
   vector<FUNC> eqns = {                                                  // equations to be solved
      []( vector<double> x ){ return x[0] * x[0] + x[1] * x[1] - 1; },
      []( vector<double> x ){ return x[0] - x[1];                   }
                       };
   vector<double> x = { 1.0, 1.0 };                                       // initial guess (gives positive root)
// vector<double> x = {-1.0,-1.0 };                                       // alternative (gives negative root)

   if ( solve( eqns, x ) )
   {
      for ( auto e : x ) cout << e << '\n';
   }
   else
   {
      cout << "Failed\n";
   }
}

//======================================================================

bool solve( vector<FUNC> eqns, vector<double> &x )
{
   const int MAXITER = 1000;
   const double TOL = 1.0e-6;
   double residuals = 1.0;
   for ( int iter = 1; iter <= MAXITER && residuals > TOL; iter++ )
   {
      auto f = getFunction( eqns, x );                       // compute f(x*)
      if ( residuals < TOL ) break;                          // if converged, finish

      auto J = getJacobian( eqns, x );                       // computed J(x*)
      vector<double> deltax( x.size() );
      if ( !GaussElimination( J, f, deltax ) ) break;        // update x -> x* - J^-1.f(x*)
      x = x - deltax;                                        //

      residuals = getResiduals( f );                         // test for convergence
   }
   return residuals <= TOL;
}

//======================================================================

vector<double> operator -( vector<double> a, const vector<double> &b )         // utility routine to subtract vectors
{
   for ( int i = 0; i < b.size(); i++ ) a[i] -= b[i];
   return a;
}

//======================================================================

vector<double> getFunction( vector<FUNC> eqns, const vector<double> &x )       // fill in current function values
{
   vector<double> f;
   for ( FUNC e : eqns ) f.push_back( e( x ) );
   return f;
}

//======================================================================

vector< vector<double> > getJacobian( vector<FUNC> eqns, const vector<double> &x )   // create Jacobian from FINITE DIFFERENCES
{
   const double SMALL = 1.0e-6;                                                      // needs improving for generality
   int N = x.size();
   vector< vector<double> > J( N, vector<double>( N ) );
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ ) 
      {
         auto xplus  = x;   xplus [j] = x[j] + 0.5 * SMALL;
         auto xminus = x;   xminus[j] = x[j] - 0.5 * SMALL;
         J[i][j] = ( eqns[i]( xplus ) - eqns[i]( xminus ) ) / SMALL;                 // finite-difference approximation for df_i/dx_j
      }
   }
   return J;
}

//======================================================================

double getResiduals( const vector<double> &x )                                       // L2 norm of vector x
{
   double res = 0.0;
   for ( auto e : x ) res += e * e;
   return sqrt( res );
}

//======================================================================

bool GaussElimination( vector< vector<double> > A, vector<double> B , vector<double> &X )
//-------------------------------------------------------------------
// Solve AX = B by Gaussian elimination
// Note: copies made of A and B
//-------------------------------------------------------------------
{
   const double SMALL = 1.0e-20;
   int n = A.size();

   // Row operations for i = 0, ,,,, n - 2 (n-1 not needed)
   for ( int i = 0; i < n - 1; i++ )
   {
      // Pivot: find row r below with largest element in column i
      int r = i;
      double maxA = abs( A[i][i] );
      for ( int k = i + 1; k < n; k++ )
      {
         double val = abs( A[k][i] );
         if ( val > maxA )
         {
            r = k;
            maxA = val;
         }
      }
      if ( r != i )
      {
         for ( int j = i; j < n; j++ ) swap( A[i][j], A[r][j] );
         swap( B[i], B[r] );
      }

      // Row operations to make upper-triangular
      double pivot = A[i][i];
      if ( abs( pivot ) < SMALL ) return false;            // Singular matrix

      for ( int r = i + 1; r < n; r++ )                    // On lower rows
      {
         double multiple = A[r][i] / pivot;                // Multiple of row i to clear element in ith column
         for ( int j = i; j < n; j++ ) A[r][j] -= multiple * A[i][j];
         B[r] -= multiple * B[i];
      }
   }

   // Check last row
   if ( abs( A[n-1][n-1] ) < SMALL ) return false;         // Singular matrix

   // Back-substitute
   X = B;
   for ( int i = n - 1; i >= 0; i-- )
   {
      for ( int j = i + 1; j < n; j++ )  X[i] -= A[i][j] * X[j];
      X[i] /= A[i][i];
   }

   return true;
}

字符串
输出量：

0.707107
0.707107

型

赞(0）分享  回复(0）举报  11个月前
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com