CYCLUS
Loading...
Searching...
No Matches
uniform_taylor.cc
Go to the documentation of this file.
1// Implements the UniformTaylor class
2#include "uniform_taylor.h"
3
4#include <cmath>
5#include <string>
6
7#include "error.h"
8
9namespace cyclus {
10
11// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
13 const double t) {
14 int n = A.NumRows();
15
16 // checks if the dimensions of A and x_o are compatible for matrix-vector
17 // computations
18 if (x_o.NumRows() != n) {
19 std::string error = "Error: Matrix-Vector dimensions are not compatible: " +
20 boost::lexical_cast<std::string>(x_o.NumRows()) +
21 " rows vs " + boost::lexical_cast<std::string>(n) +
22 " nuclides.";
23 throw ValueError(error);
24 }
25
26 // step 1 of algorithm: calculates the largest diagonal element (alpha)
27 double alpha = MaxAbsDiag(A);
28
29 // step 2 of algorithm: creates the matrix B = A + alpha * I
30 Matrix B = identity(n);
31 B = alpha * B;
32 B += A;
33
34 // steps 3-7 of algorithm: computes the solution Vector x_t
35 double tol = 1e-3;
36 Vector x_t = x_o;
37
38 x_t = GetSolutionVector(B, x_o, alpha, t, tol);
39
40 return x_t;
41}
42
43// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44double UniformTaylor::MaxAbsDiag(const Matrix& A) {
45 int n = A.NumRows(); // stores the order of the matrix A
46 double a_ii = A(1, 1); // begins with the first diagonal element
47
48 // Initializes the maximum diagonal element to the absolute value of the
49 // first diagonal element a_ii
50 a_ii = fabs(a_ii);
51 double max_a_ii = a_ii;
52
53 // Searches the remaining diagonal elements for the largest absolute value
54 for (int i = 2; i <= n; ++i) {
55 a_ii = A(i, i);
56 a_ii = fabs(a_ii);
57
58 if (a_ii > max_a_ii) {
59 max_a_ii = a_ii;
60 }
61 }
62
63 return max_a_ii;
64}
65
66// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67Vector UniformTaylor::GetSolutionVector(const Matrix& B, const Vector& x_o,
68 double alpha, double t, double tol) {
69 // step 3 of algorithm: calculates exp( -alpha * t)
70 long double alpha_t = alpha * t;
71 long double expat = exp(-alpha_t);
72
73 if (expat == 0) {
74 std::string error =
75 "Error: exp(-alpha * t) exceeds the range of a long double.";
76 error += "\nThe Uniform Taylor method cannot solve the matrix exponential.";
77 throw ValueError(error);
78 }
79
80 // step 4 of algorithm: initializes the next term Ck, the total sum of Ck
81 // terms, and the previous term Ck-1
82 //
83 // NOTE: the exponential term is included at the beginning of the sum so
84 // that Ck_sum slowly gets larger as more terms are added, rather than
85 // simply multiplying a very small number by a very big one at the end
86 Vector C_next = expat * x_o;
87 Vector Ck_sum = C_next;
88 Vector C_prev = C_next;
89
90 // step 5 of algorithm: determines the maximum number of terms needed
91 int maxTerms = MaxNumTerms(alpha_t, tol);
92
93 // step 6 of algorithm: computes the sum of Ck terms until the maximum
94 // number of terms has been reached
95 for (int k = 1; k < maxTerms; ++k) {
96 // step 6a of algorithm: computes the next term in the series
97 C_next = (t / k) * B;
98 C_next *= C_prev;
99
100 // step 6b of algorithm: updates the solution Ck_sum
101 Ck_sum += C_next;
102
103 // step 6c of algorithm: resets the previous term for the next iteration
104 C_prev = C_next;
105 }
106
107 // step 7 of algorithm: returns the solution for x_t
108 return Ck_sum;
109}
110
111// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112int UniformTaylor::MaxNumTerms(long double alpha_t, double epsilon) {
113 long double nextTerm; // stores the next term in the series
114
115 // initializes the previous term and the sum of terms in the series
116 long double prevTerm = 1;
117 long double sumTerms = 1;
118
119 // calculates the lower bound of the series
120 long double lowerBound = exp(alpha_t);
121
122 // checks to see if exp(alpha * t) is infinite
123 if (lowerBound == HUGE_VAL) {
124 std::string error =
125 "Error: exp(alpha * t) exceeds the range of a long double";
126 error += "\nThe Uniform Taylor method cannot solve the matrix exponential.";
127 throw ValueError(error);
128 }
129
130 lowerBound = lowerBound * (1 - epsilon);
131
132 // computes the sum of terms until it is greater than the lower bound
133 int p = 1;
134 bool stopSum = false;
135
136 while (stopSum != true) {
137 // checks if the sum of terms is greater than the lower error bound
138 if (sumTerms >= lowerBound) {
139 stopSum = true;
140 } else {
141 // computes the next term in the series
142 nextTerm = (alpha_t / p) * prevTerm;
143
144 // adds the next term to the sum of terms and updates the total number
145 // of terms in the series
146 sumTerms += nextTerm;
147 ++p;
148
149 // resets the previous term for the next iteration of the loop
150 prevTerm = nextTerm;
151 }
152 }
153
154 return p;
155}
156
157} // namespace cyclus
int NumRows() const
Definition l_matrix.cc:58
static Vector MatrixExpSolver(const Matrix &A, const Vector &x_o, const double t)
Solves the matrix exponential problem:
For values that are too big, too small, etc.
Definition error.h:37
taken directly from OsiSolverInterface.cpp on 2/17/14 from https://projects.coin-or....
Definition agent.cc:14
LMatrix Matrix
LMatrix identity(int n)
Definition l_matrix.cc:284
LMatrix Vector
struct pyne::alpha alpha
a struct matching the '/decay/alphas' table in nuc_data.h5.