LMCAS
Lamina Computer Algebra System
载入中...
搜索中...
未找到
LMCAS::ode_root_detail 命名空间参考

命名空间

namespace  anonymous_namespace{ode_characteristic_roots.cpp}

类

struct  CharRoot
struct  Complex

类型定义

using RealPolynomial = std::vector<double>

函数

RealPolynomial approximate_polynomial_gcd (RealPolynomial lhs, RealPolynomial rhs)
std::pair< RealPolynomial, RealPolynomial > divide_polynomials (const RealPolynomial &dividend, const RealPolynomial &divisor)
Result< std::vector< CharRoot > > find_characteristic_roots (const std::vector< double > &coeffs, const std::string &operation)
 使用数值方法求解特征多项式的所有根.
bool is_exact_zero_polynomial_root (const RealPolynomial &coefficients, const Complex &root)
RealPolynomial monic_polynomial (RealPolynomial polynomial)
RealPolynomial polynomial_derivative (const RealPolynomial &polynomial)
bool polynomial_remainder_is_small (const RealPolynomial &remainder, const RealPolynomial &dividend)
double polynomial_root_backward_error (const RealPolynomial &coefficients, const Complex &root)
double polynomial_scale (const RealPolynomial &polynomial)
void trim_polynomial (RealPolynomial &polynomial)

类说明

◆ LMCAS::ode_root_detail::CharRoot

struct LMCAS::ode_root_detail::CharRoot

在文件 ode_polynomial_utils.hpp 第 12 行定义.

LMCAS::ode_root_detail::CharRoot 的协作图:
类成员
double imag_part
bool is_complex
int multiplicity
double real_part

类型定义说明

◆ RealPolynomial

using LMCAS::ode_root_detail::RealPolynomial = std::vector<double>

在文件 ode_polynomial_utils.hpp 第 10 行定义.

函数说明

◆ approximate_polynomial_gcd()

RealPolynomial LMCAS::ode_root_detail::approximate_polynomial_gcd ( RealPolynomial lhs,
RealPolynomial rhs )
inline

在文件 ode_polynomial_utils.hpp 第 150 行定义.

151 {
152 lhs = monic_polynomial(std::move(lhs));
153 rhs = monic_polynomial(std::move(rhs));
154 for (int iteration = 0; iteration < 8 && rhs.size() > 1; ++iteration) {
155 auto division = divide_polynomials(lhs, rhs);
156 RealPolynomial remainder = std::move(division.second);
157 if (remainder.size() == 1 && remainder.front() == 0.0) {
158 return rhs;
159 }
160 lhs = std::move(rhs);
161 rhs = monic_polynomial(std::move(remainder));
162 }
163 return {1.0};
164}
std::pair< RealPolynomial, RealPolynomial > divide_polynomials(const RealPolynomial &dividend, const RealPolynomial &divisor)
RealPolynomial monic_polynomial(RealPolynomial polynomial)
const ExprPtr & rhs

引用了 approximate_polynomial_gcd(), divide_polynomials(), monic_polynomial() , 以及 LMCAS::rhs.

被这些函数引用 approximate_polynomial_gcd() , 以及 LMCAS::ode_root_detail::anonymous_namespace{ode_characteristic_roots.cpp}::repeated_characteristic_roots().

函数调用图:
这是这个函数的调用关系图:

◆ divide_polynomials()

std::pair< RealPolynomial, RealPolynomial > LMCAS::ode_root_detail::divide_polynomials ( const RealPolynomial & dividend,
const RealPolynomial & divisor )
inline

在文件 ode_polynomial_utils.hpp 第 98 行定义.

99 {
100 if (divisor.empty() || divisor.front() == 0.0 ||
101 dividend.size() < divisor.size()) {
102 return {{0.0}, dividend};
103 }
104
105 RealPolynomial remainder = dividend;
106 RealPolynomial quotient(
107 dividend.size() - divisor.size() + 1, 0.0);
108 for (size_t i = 0; i < quotient.size(); ++i) {
109 const double factor = remainder[i] / divisor.front();
110 quotient[i] = factor;
111 for (size_t j = 0; j < divisor.size(); ++j) {
112 remainder[i + j] -= factor * divisor[j];
113 }
114 }
115 remainder.erase(
116 remainder.begin(), remainder.begin() + quotient.size());
117 if (remainder.empty()) {
118 remainder.push_back(0.0);
119 }
120 trim_polynomial(quotient);
121 trim_polynomial(remainder);
122 return {std::move(quotient), std::move(remainder)};
123}
void trim_polynomial(RealPolynomial &polynomial)

引用了 divide_polynomials() , 以及 trim_polynomial().

被这些函数引用 approximate_polynomial_gcd(), divide_polynomials(), LMCAS::ode_root_detail::anonymous_namespace{ode_characteristic_roots.cpp}::recover_root_multiplicities() , 以及 LMCAS::ode_root_detail::anonymous_namespace{ode_characteristic_roots.cpp}::repeated_characteristic_roots().

函数调用图:
这是这个函数的调用关系图:

◆ find_characteristic_roots()

Result< std::vector< CharRoot > > LMCAS::ode_root_detail::find_characteristic_roots ( const std::vector< double > & coeffs,
const std::string & operation )

使用数值方法求解特征多项式的所有根.

对于一、二次多项式使用解析公式;对于三至六次多项式使用 Durand-Kerner 迭代。迭代前通过多项式与其导数的近似 GCD 识别并提取重根。

在文件 ode_characteristic_roots.cpp 第 338 行定义.

340{
341 int n = static_cast<int>(coeffs.size()) - 1;
342 if (n <= 0) {
343 return Result<std::vector<CharRoot>>::failure(
345 "characteristic polynomial has no roots", operation);
346 }
347
348 std::vector<CharRoot> roots;
349 double leading = coeffs[0];
350 if (leading == 0.0) {
351 return Result<std::vector<CharRoot>>::failure(
353 "characteristic polynomial has no usable leading coefficient",
354 operation);
355 }
356 std::vector<double> norm_coeffs(coeffs.size());
357 for (size_t i = 0; i < coeffs.size(); ++i) {
358 norm_coeffs[i] = coeffs[i] / leading;
359 if (!std::isfinite(norm_coeffs[i])) {
360 return Result<std::vector<CharRoot>>::failure(
362 "characteristic polynomial normalization is non-finite",
363 operation);
364 }
365 }
366
367 if (n == 1) {
368 roots.push_back({-norm_coeffs[1], 0.0, 1, false});
369 return roots;
370 }
371 if (n == 2) {
372 return quadratic_characteristic_roots(norm_coeffs, operation);
373 }
374 auto scale_result = characteristic_root_scale(norm_coeffs, operation);
375 if (!scale_result) {
376 return Result<std::vector<CharRoot>>::failure(scale_result.error());
377 }
378 const double root_scale = scale_result.value();
379 std::vector<double> solver_coeffs = norm_coeffs;
380 for (int i = 1; i <= n; ++i) {
381 for (int power = 0; power < i; ++power) {
382 solver_coeffs[i] /= root_scale;
383 }
384 }
385
386 auto repeated = repeated_characteristic_roots(solver_coeffs, operation);
387 if (!repeated) {
388 return repeated;
389 }
390 if (!repeated.value().empty()) {
391 return restore_root_scale(std::move(repeated.value()), root_scale, operation);
392 }
393 auto candidates = iterate_characteristic_roots(solver_coeffs, operation);
394 if (!candidates) {
395 return Result<std::vector<CharRoot>>::failure(candidates.error());
396 }
397 auto classified = classify_characteristic_roots(
398 solver_coeffs, candidates.value(), operation);
399 if (!classified) {
400 return classified;
401 }
402 return restore_root_scale(std::move(classified.value()), root_scale, operation);
403}
Result< std::vector< CharRoot > > quadratic_characteristic_roots(const RealPolynomial &norm_coeffs, const std::string &operation)
Result< std::vector< CharRoot > > classify_characteristic_roots(const RealPolynomial &solver_coeffs, const std::vector< Complex > &z, const std::string &operation)
Result< std::vector< Complex > > iterate_characteristic_roots(const RealPolynomial &solver_coeffs, const std::string &operation)
Result< double > characteristic_root_scale(const RealPolynomial &norm_coeffs, const std::string &operation)
Result< std::vector< CharRoot > > repeated_characteristic_roots(const RealPolynomial &solver_coeffs, const std::string &operation)
Result< std::vector< CharRoot > > restore_root_scale(std::vector< CharRoot > scaled_roots, double root_scale, const std::string &operation)
ExprSetResult roots(const ExprPtr &expression, const std::string &variable, ComputationContext &context, const SolveOptions &options)

引用了 find_characteristic_roots(), LMCAS::NumericFailure , 以及 LMCAS::roots().

被这些函数引用 find_characteristic_roots(), LMCAS::ode_root_detail::anonymous_namespace{ode_characteristic_roots.cpp}::repeated_characteristic_roots(), LMCAS::solve_euler_ode_impl(), LMCAS::solve_higher_order_ode_impl() , 以及 LMCAS::solve_linear2_homogeneous().

函数调用图:
这是这个函数的调用关系图:

◆ is_exact_zero_polynomial_root()

bool LMCAS::ode_root_detail::is_exact_zero_polynomial_root ( const RealPolynomial & coefficients,
const Complex & root )
inline

在文件 ode_polynomial_utils.hpp 第 49 行定义.

50 {
51 return root.re == 0.0 && root.im == 0.0 && coefficients.back() == 0.0;
52}

引用了 LMCAS::ode_root_detail::Complex::im, is_exact_zero_polynomial_root() , 以及 LMCAS::ode_root_detail::Complex::re.

被这些函数引用 is_exact_zero_polynomial_root() , 以及 polynomial_root_backward_error().

函数调用图:
这是这个函数的调用关系图:

◆ monic_polynomial()

RealPolynomial LMCAS::ode_root_detail::monic_polynomial ( RealPolynomial polynomial)
inline

在文件 ode_polynomial_utils.hpp 第 125 行定义.

125 {
126 trim_polynomial(polynomial);
127 if (polynomial.empty() || polynomial.front() == 0.0) {
128 return {0.0};
129 }
130 const double leading = polynomial.front();
131 for (double& coefficient : polynomial) coefficient /= leading;
132 return polynomial;
133}

引用了 monic_polynomial() , 以及 trim_polynomial().

被这些函数引用 approximate_polynomial_gcd() , 以及 monic_polynomial().

函数调用图:
这是这个函数的调用关系图:

◆ polynomial_derivative()

RealPolynomial LMCAS::ode_root_detail::polynomial_derivative ( const RealPolynomial & polynomial)
inline

在文件 ode_polynomial_utils.hpp 第 135 行定义.

136 {
137 if (polynomial.size() <= 1) {
138 return {0.0};
139 }
140 RealPolynomial derivative;
141 derivative.reserve(polynomial.size() - 1);
142 const size_t degree = polynomial.size() - 1;
143 for (size_t i = 0; i < degree; ++i) {
144 derivative.push_back(
145 polynomial[i] * static_cast<double>(degree - i));
146 }
147 return derivative;
148}

引用了 polynomial_derivative().

被这些函数引用 polynomial_derivative() , 以及 LMCAS::ode_root_detail::anonymous_namespace{ode_characteristic_roots.cpp}::repeated_characteristic_roots().

函数调用图:
这是这个函数的调用关系图:

◆ polynomial_remainder_is_small()

bool LMCAS::ode_root_detail::polynomial_remainder_is_small ( const RealPolynomial & remainder,
const RealPolynomial & dividend )
inline

在文件 ode_polynomial_utils.hpp 第 166 行定义.

167 {
168 return polynomial_scale(remainder) <=
169 std::max(1.0, polynomial_scale(dividend)) * 1e-7;
170}
double polynomial_scale(const RealPolynomial &polynomial)
ExprResult e()

引用了 LMCAS::e(), polynomial_remainder_is_small() , 以及 polynomial_scale().

被这些函数引用 polynomial_remainder_is_small(), LMCAS::ode_root_detail::anonymous_namespace{ode_characteristic_roots.cpp}::recover_root_multiplicities() , 以及 LMCAS::ode_root_detail::anonymous_namespace{ode_characteristic_roots.cpp}::repeated_characteristic_roots().

函数调用图:
这是这个函数的调用关系图:

◆ polynomial_root_backward_error()

double LMCAS::ode_root_detail::polynomial_root_backward_error ( const RealPolynomial & coefficients,
const Complex & root )
inline

在文件 ode_polynomial_utils.hpp 第 54 行定义.

55 {
56 Complex value{coefficients.front(), 0.0};
57 double scale = std::abs(coefficients.front());
58 const double root_magnitude = root.magnitude();
59 for (size_t i = 1; i < coefficients.size(); ++i) {
60 value = value * root + Complex{coefficients[i], 0.0};
61 scale = scale * root_magnitude + std::abs(coefficients[i]);
62 }
63 if (!std::isfinite(value.re) || !std::isfinite(value.im) ||
64 !std::isfinite(scale)) {
65 return INFINITY;
66 }
67 if (scale == 0.0) {
68 // At an exact zero root with zero constant term, both the residual
69 // and its coefficientwise bound vanish: no perturbation is needed.
70 // A nonzero root can instead underflow the bound, which is not proof.
71 return is_exact_zero_polynomial_root(coefficients, root)
72 ? 0.0 : INFINITY;
73 }
74 return value.magnitude() / scale;
75}
bool is_exact_zero_polynomial_root(const RealPolynomial &coefficients, const Complex &root)
const std::string const ExprPtr & value

引用了 is_exact_zero_polynomial_root(), LMCAS::ode_root_detail::Complex::magnitude(), polynomial_root_backward_error() , 以及 LMCAS::value.

被这些函数引用 LMCAS::ode_root_detail::anonymous_namespace{ode_characteristic_roots.cpp}::classify_characteristic_roots(), LMCAS::ode_root_detail::anonymous_namespace{ode_characteristic_roots.cpp}::iterate_characteristic_roots(), polynomial_root_backward_error() , 以及 LMCAS::ode_root_detail::anonymous_namespace{ode_characteristic_roots.cpp}::quadratic_characteristic_roots().

函数调用图:
这是这个函数的调用关系图:

◆ polynomial_scale()

double LMCAS::ode_root_detail::polynomial_scale ( const RealPolynomial & polynomial)
inline

在文件 ode_polynomial_utils.hpp 第 77 行定义.

77 {
78 double scale = 0.0;
79 for (double coefficient : polynomial) {
80 scale = std::max(scale, std::abs(coefficient));
81 }
82 return scale;
83}

引用了 polynomial_scale().

被这些函数引用 polynomial_remainder_is_small(), polynomial_scale() , 以及 trim_polynomial().

函数调用图:
这是这个函数的调用关系图:

◆ trim_polynomial()

void LMCAS::ode_root_detail::trim_polynomial ( RealPolynomial & polynomial)
inline

在文件 ode_polynomial_utils.hpp 第 85 行定义.

85 {
86 const double scale = polynomial_scale(polynomial);
87 const double tolerance = std::max(1.0, scale) * 1e-10;
88 auto first = polynomial.begin();
89 while (first + 1 != polynomial.end() && std::abs(*first) <= tolerance) {
90 ++first;
91 }
92 polynomial.erase(polynomial.begin(), first);
93 if (polynomial.empty()) {
94 polynomial.push_back(0.0);
95 }
96}

引用了 LMCAS::e(), polynomial_scale() , 以及 trim_polynomial().

被这些函数引用 divide_polynomials(), monic_polynomial() , 以及 trim_polynomial().

函数调用图:
这是这个函数的调用关系图: