finitediff.cpp
1 
2 #include <symengine/add.h>
3 #include <symengine/mul.h>
4 #include <symengine/constants.h>
5 
6 namespace SymEngine
7 {
8 
9 vec_basic generate_fdiff_weights_vector(const vec_basic &grid,
10  const unsigned max_deriv,
11  const RCP<const Basic> around)
12 {
13  // Parameters
14  // ----------
15  // grid: grid point locations
16  // max_deriv: highest derivative.
17  // around: location where approximations are to be accurate
18  //
19  // Returns
20  // -------
21  // weights[grid_index, deriv_order]: weights of order 0 to max_deriv (column
22  // major order)
23  //
24  // References
25  // ----------
26  // Generation of Finite Difference Formulas on Arbitrarily Spaced Grids
27  // Bengt Fornberg, Mathematics of compuation, 51, 184, 1988, 699-706
28  //
29  const unsigned len_g = numeric_cast<unsigned>(grid.size());
30  const unsigned len_w = len_g * (max_deriv + 1);
31  RCP<const Basic> c1, c4, c5;
32  c1 = one;
33  c4 = sub(grid[0], around);
34  vec_basic weights(len_w);
35  weights[0] = one;
36  for (unsigned idx = 1; idx < len_w; ++idx)
37  weights[idx] = zero; // clear weights
38  for (unsigned i = 1; i < len_g; ++i) {
39  const int mn = (i < max_deriv) ? i : max_deriv; // min(i, max_deriv)
40  RCP<const Basic> c2 = one;
41  c5 = c4;
42  c4 = sub(grid[i], around);
43  for (unsigned j = 0; j < i; ++j) {
44  const RCP<const Basic> c3 = sub(grid[i], grid[j]);
45  c2 = mul(c2, c3);
46  if (j == i - 1) {
47  for (int k = mn; k >= 1; --k) {
48  weights[i + k * len_g]
49  = div(mul(c1, sub(mul(integer(k),
50  weights[i - 1 + (k - 1) * len_g]),
51  mul(c5, weights[i - 1 + k * len_g]))),
52  c2);
53  }
54  weights[i]
55  = mul(minus_one, div(mul(c1, mul(c5, weights[i - 1])), c2));
56  }
57  for (int k = mn; k >= 1; --k) {
58  weights[j + k * len_g]
59  = div(sub(mul(c4, weights[j + k * len_g]),
60  mul(integer(k), weights[j + (k - 1) * len_g])),
61  c3);
62  }
63  weights[j] = div(mul(c4, weights[j]), c3);
64  }
65  c1 = c2;
66  }
67  return weights;
68 }
69 } // namespace SymEngine
Classes and functions relating to the binary operation of addition.
Main namespace for SymEngine package.
Definition: add.cpp:19
RCP< const Basic > div(const RCP< const Basic > &a, const RCP< const Basic > &b)
Division.
Definition: mul.cpp:426
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
Definition: integer.h:200
RCP< const Basic > sub(const RCP< const Basic > &a, const RCP< const Basic > &b)
Substracts b from a.
Definition: add.cpp:495
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:347