Barycentric rational interpolation-CodePudding

The

EXAMPLE:

for nodes: {{1, 3}, {2, 5}, {4, 4}, {5, 2}, {7, 1}}, and row=2

std::cout<<f(2.5) would print 5.425

Code:

#include <iostream> #include <stdexcept> #include <utility> #include <vector> #include <functional> #include <cmath> int max(int a, int b) { if (a > b)return a; return b; } int min(int a, int b) { if (a < b)return a; return b; } std::vector<double> ComputeWeights(const std::vector<std::pair<int, int>>& nodes, int d) { std::vector<double>w; int n = nodes.size(); for (int i = 0; i < n; i ) { double sum = 0; int k = max(1, i - d); double product = product_f(k, k d); for (int j = k; j < k d; j ) sum = (pow(-1, k - 1) * 1. / (nodes[i].first - nodes[k].first)); w.push_back(sum); } return w; } double ComputeF(double x, const std::vector<std::pair<int, int>>& nodes, const std::vector<double>& w) { double f, ff = 0; int n = nodes.size(); for (int i = 0; i < n; i ) { f = ((w[i] * nodes[i].second) / (x - nodes[i].first)) / (w[i] / (x - nodes[i].first)); ff = f; } return ff; } std::function<double(double)> formula(std::vector<std::pair<int, int>>nodes, int d) { if (d < 0 || d > nodes.size()) throw std::domain_error("Forbidden row"); int n = nodes.size(); for (int i = 0; i < n; i ) for (int j = i 1; j < n; j ) if (nodes[j].first == nodes[i].first) throw std::domain_error("Forbidden coordinates"); auto w = ComputeWeights(nodes, d); return [nodes, w](double x) { return ComputeF(x, nodes, w); }; } int main () { auto f = formula({{1, 3}, {2, 5}, {4, 4}, {5, 2}, {7, 1}}, 2); std::cout << f(2.5); return 0; }

Could you help me to write this function properly?

My output is: -nan (correct is 5.425)

CodePudding user response：

Here you go. Note that I rewrote basic parts of your code using my own code from the previous edit version. Further, to clarify: You are actually looking for a thing named Floater-Hormann approximation, which is basically an approximation by a rational function without poles in the interpolation region (and that is usually evaluated using the barycentric formula).

#include <iostream> #include <stdexcept> #include <utility> #include <vector> #include <functional> #include <cmath> std::vector<double> ComputeWeights(const std::vector<std::pair<double, double>>& nodes, int d) { int n = (int)nodes.size(); std::vector<double> w(n); for (int k = 0; k < n; k) { int imin = std::max(k - d, 0); int imax = std::min(n - d - 1,k); double temp = imin & 1 ? -1.0 : 1.0; double sum = 0.0; for (int i = imin; i <= imax; i) { int jmax = std::min(i d, n - 1); double term = 1.0; for (int j = i; j <= jmax; j) { if (j == k) continue; term *= (nodes[k].first - nodes[j].first); } term = temp / term; temp = -temp; sum = term; } w[k] = sum; } return w; } double ComputeF(double x, const std::vector<std::pair<double, double>>& nodes, const std::vector<double>& w) { double num = 0.0; double denom = 0.0; for (int i = 0; i < (int)nodes.size(); i) { if (x == nodes[i].first) { return nodes[i].second; } auto ad = w[i] / (x - nodes[i].first); num = ad * nodes[i].second; denom = ad; } return num / denom; } auto formula(std::vector<std::pair<double, double>>nodes, int d) { if (d < 0 || d > (int)nodes.size()) throw std::domain_error("Forbidden row"); int n = nodes.size(); for (int i = 0; i < n; i ) for (int j = i 1; j < n; j ) if (nodes[j].first == nodes[i].first) throw std::domain_error("Forbidden coordinates"); auto w = ComputeWeights(nodes, d); return [nodes, w](double x) { return ComputeF(x, nodes, w); }; }

Further, on the C level, note that I didn't use a std::function, but directly the lambda, which doesn't add any overhead by indirection. Moreover, I used a std::pair<double, double>, as we are usually in the floating-point regime when it comes to approximation.

The result is, as required

int main() { std::vector<std::pair<double, double> > nodes = {{1, 3}, {2, 5}, {4, 4}, {5, 2}, {7, 1}}; auto f = formula(nodes,2); std::cout<<f(2.5)<<std::endl; //Prints 5.425 }

DEMO

CodePudding user response：

You forgot return keyword inside lambda body return ((sum_f(1, n)...