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