假设我有一个方程组:
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;
}
型
2条答案
按热度按时间ldfqzlk81#
尽管你没有真正解释代码中“尴尬”的地方。
我看了一下,并能够提出一些优化和简化(例如使用C++17语言功能)。
注意,为了清晰(和安全),我使用命名空间别名:
字符串
优化是为了避免向量/矩阵的无限数组存储,因为它们是固定大小的:
型
我选择了
std::array
,因为它很简单,并且对constexpr友好,并且(与U::bounded_array<>
相反)允许构造初始化器列表。这将允许我们更自然的赋值:型
另一部分是使用结构化绑定,而不是
get<n>
或[n]
舞蹈。从效率的Angular 来看,我将验证编译器是否消除了矩阵初始化器中的公共子表达式。
对于
f0
/f1
函数,我让类型演绎来完成工作:型
主驱动程序没有太大变化,尽管我重新塑造了循环,使其更传统:
Live On Coliru
型
印刷
型
摘要
我希望这些更改能使代码稍微快一点(看到没有I/O的微基准测试会很有趣)。
我不确定这些是否解决了你在问题中提到的一些“尴尬”。
ltskdhd12#
我想你可以不用升压,自己卷。
方程采用泛函建立,雅可比矩阵采用有限差分法求解,但还有待改进。
字符串
输出量:
型