Searching...
No Matches
finitediff.cpp
1
3#include <symengine/mul.h>
5
6namespace SymEngine
7{
8
9vec_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.