这是一个我已经挣扎了一周的问题，在浪费了几个小时后回来只是为了放弃。。。
我应该找到以下拉盖尔多项式的系数：
P0(x) = 1
P1(x) = 1 - x
Pn(x) = ((2n - 1 - x) / n) * P(n-1) - ((n - 1) / n) * P(n-2)
我相信我的实现中有一个错误，因为出于某种原因，我得到的系数似乎太大了。这是该程序生成的输出：
a1 = -190.234
a2 = -295.833
a3 = 378.283
a4 = -939.537
a5 = 774.861
a6 = -400.612
代码描述(如下)：
如果你把代码向下滚动到我声明数组的部分，你会发现给定的x和y。
函数多项式只是用某个x的多项式值填充一个数组。它是一个递归函数。我相信它运行良好，因为我已经检查了输出值。
高斯函数通过对输出阵列进行高斯消去来求出系数。我认为这就是问题的开始。我想知道，这段代码中是否有错误，或者我的计算结果的方法可能不好？我正试图这样验证它们：
-190.234 * 1.5 ^ 5 - 295.833 * 1.5 ^ 4 ... - 400.612 = -3017,817625 =/= 2
代码：
#include "stdafx.h"
#include <conio.h>
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
double polynomial(int i, int j, double **tab)
{
double n = i;
double **array = tab;
double x = array[j][0];
if (i == 0) {
return 1;
} else if (i == 1) {
return 1 - x;
} else {
double minusone = polynomial(i - 1, j, array);
double minustwo = polynomial(i - 2, j, array);
double result = (((2.0 * n) - 1 - x) / n) * minusone - ((n - 1.0) / n) * minustwo;
return result;
}
}
int gauss(int n, double tab[6][7], double results[7])
{
double multiplier, divider;
for (int m = 0; m <= n; m++)
{
for (int i = m + 1; i <= n; i++)
{
multiplier = tab[i][m];
divider = tab[m][m];
if (divider == 0) {
return 1;
}
for (int j = m; j <= n; j++)
{
if (i == n) {
break;
}
tab[i][j] = (tab[m][j] * multiplier / divider) - tab[i][j];
}
for (int j = m; j <= n; j++) {
tab[i - 1][j] = tab[i - 1][j] / divider;
}
}
}
double s = 0;
results[n - 1] = tab[n - 1][n];
int y = 0;
for (int i = n-2; i >= 0; i--)
{
s = 0;
y++;
for (int x = 0; x < n; x++)
{
s = s + (tab[i][n - 1 - x] * results[n-(x + 1)]);
if (y == x + 1) {
break;
}
}
results[i] = tab[i][n] - s;
}
}
int _tmain(int argc, _TCHAR* argv[])
{
int num;
double **array;
array = new double*[5];
for (int i = 0; i <= 5; i++)
{
array[i] = new double[2];
}
//i 0 1 2 3 4 5
array[0][0] = 1.5; //xi 1.5 2 2.5 3.5 3.8 4.1
array[0][1] = 2; //yi 2 5 -1 0.5 3 7
array[1][0] = 2;
array[1][1] = 5;
array[2][0] = 2.5;
array[2][1] = -1;
array[3][0] = 3.5;
array[3][1] = 0.5;
array[4][0] = 3.8;
array[4][1] = 3;
array[5][0] = 4.1;
array[5][1] = 7;
double W[6][7]; //n + 1
for (int i = 0; i <= 5; i++)
{
for (int j = 0; j <= 5; j++)
{
W[i][j] = polynomial(j, i, array);
}
W[i][6] = array[i][1];
}
for (int i = 0; i <= 5; i++)
{
for (int j = 0; j <= 6; j++)
{
cout << W[i][j] << "t";
}
cout << endl;
}
double results[6];
gauss(6, W, results);
for (int i = 0; i < 6; i++) {
cout << "a" << i + 1 << " = " << results[i] << endl;
}
_getch();
return 0;
}
我认为你对递归多项式生成的解释要么需要修改，要么对我来说有点太聪明了。
给定p[0][5]＝{1,0,0,0,0，…}；P[1][5]={1，-1,0,0,0，…}
则P[2]是*P[0]+卷积(P[1]，{c，d})
其中a=-((n-1)/n)c=(2n-1)/n和d=-1/n
这可以推广为：p[n]=a*p[n-2]+conv(p[n-1]，{c，d})；在每一步中，都涉及到与(c+d*x)的多项式乘法，它将次数增加一(仅增加一…)，并与乘以标量a的P[n-1]相加。
那么插值因子x很可能在[0..1]的范围内。
(卷积意味着，你应该实现多项式乘法，幸运的是这很容易…)
[a,b,c,d]
* [e,f]
------------------
af,bf,cf,df +
ae,be,ce,de, 0 +
--------------------------
(= coefficients of the final polynomial)
P1(x) = x - 1的定义没有按规定实现。计算中有1 - x。
我没有再看下去。