重心有理插值

Question

我正在尝试编写一个函数，该函数 returns 是一个用于在 C++ 中计算重心有理插值的函数。

注意：重新计算权重系数根本不是一个明智的想法 W, = 1, 2, ..., 在函数本身内 returns 作为“重心插值”的结果。也就是说，在这种方式下，每次调用这样的返回函数时，都会重新计算这些系数，如果它很大，这是不可接受的时间浪费。相反，这些系数应该在“Barricentric Interpolation”函数中计算，放置在某处，并“陷入”函数 returns 作为结果。这样每次调用返回函数都会使用预先计算好的系数，非常节省时间

公式：

示例：

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

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

代码：

#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;
}

你能帮我把这个函数写好吗？

我的输出是：-nan（正确的是5.425）

Answer 1

您在 lambda 正文中忘记了 return 关键字 return ((sum_f(1, n)...

Answer 2

给你。请注意，我使用以前编辑版本中我自己的代码重写了代码的基本部分。此外，澄清一下：您实际上是在寻找一个名为 Floater-Hormann 近似值的东西，它基本上是插值区域中没有极点的有理函数的近似值（通常使用重心公式计算）。

#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); };
}

此外，在 C++ 级别，请注意我没有使用 std::function，而是直接使用 lambda，这不会通过间接方式增加任何开销。此外，我使用了 std::pair<double, double>，因为我们通常在 floating-point 范围内进行近似。

结果是，按要求

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

重心有理插值

Barycentric rational interpolation

c++

generics

lambda