我正在嘗試撰寫一個函式，該函式回傳一個用于在 C 中計算重心有理插值的函式。
注意：在作為“重心插值”結果回傳的函式本身內重新計算權重系數 W、?? = 1、2、...、?? 根本不是一個明智的想法。即這樣，每次呼叫這樣的回傳函式時，都會重新計算這些系數，如果它很大，這是不可接受的浪費時間。相反，這些系數應該在“Barricentric Interpolation”函式中計算，放置在某個地方，并“困”在作為結果回傳的函式中。這樣，每次呼叫回傳的函式都會使用預先計算的系數，這樣可以節省大量時間
公式：
例子：
對于節點：{{1, 3}, {2, 5}, {4, 4}, {5, 2}, {7, 1}} 和 row=2
std::cout<<f(2.5)將列印 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）
uj5u.com熱心網友回復：
干得好。請注意，我使用以前編輯版本中的自己的代碼重寫了代碼的基本部分。此外，澄清一下：您實際上正在尋找一種名為 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 級別上，請注意我沒有使用 a std::function，而是直接使用 lambda，它不會通過間接增加任何開銷。此外，我使用了 a std::pair<double, double>，因為在進行近似時我們通常處于浮點狀態。
結果是，按要求
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
}
演示
uj5u.com熱心網友回復：
您忘記了 lambda 正文中的 return 關鍵字return ((sum_f(1, n)...
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/477888.html