多项式外推法实现代码的并行化.PDF

多项式外推法实现代码的并行化 Clay P. Breshears 问题描述 多项式插值法是指对于一个有 N 个点的集合 （x_i,y_i），找出一个多项式函数 P，该函数 经过集合中的每一个点。即，y_i = P（x_i）,1 = i = N。而多项式外推法是指沿着 X 轴 的方向，在集合外寻找函数 P 的一个解，即 X 点对应的函数值。本文描述的并行代码主要 是利用外推法对插值多项式中的点进行求值。 利用拉格朗日公式，可以对多项式 P 进行插值得到 P(X)。这种方法需要计算一个 N 项相 加的 n – 1 阶多项式。每一项的计算包含 2n – 2 次减法和 2n – 1 次乘法。这种公式 中有很多独立求解的部分，但没有办法估算误差，同时为了确保每一项得到正确的值，需要 编写大量复杂的代码。 另一种是内维尔算法。这种公式使用有限差分方法，通过逐步求精高阶多项式优化函数值。 此处提供的并行解决方案将从内维尔算法的串行代码实现开始。该串行代码参见“Numerical Recipes”，Flannery, Teukolsky ，Vetterling（剑桥大学出版社，1986年）。一些算法的 基本数学细节不在本文的讨论范围内。 串行代码 内维尔算法首先在给定数据集中 N个点 X 的零阶多项式（即常量函数）的函数值 X。假定 P1 是 P(X) = y_1 中 X 的值，并且 P2 是 P(X) = y_2 中 X 的值，从 P3 到 PN 以此类推。 一旦我们算出零阶多项式的值，就可以算出 1 阶多项式的值。假定 P12 是通过 (x_1,y_1) 和 (x_2,y_2) 两点的一阶多项式中 X 的值 ，P23 是通过 (x_2,y_2) 和 (x_3,y_3) 两点一阶多项式中 X 的值，以此类推 P34 ... P（N-1）N。用二阶多项式中的 X 值求三 阶多项式中的 X 值，再用三阶多项式中的 X 值求四阶多项式中的 X 值，以此类推，用 N-2 阶多项式中的 X 值来求 N-1 阶多项式中的 X 值。从这一点来看，它和采用拉格朗日公式 的算法差别并不大，但它要应用公式 O(N^2) 次！ 下表显示了由 4 个初始值推导后续多项式的过程（从左至右）。 x1: y1 = P1 P12 x2: y2 = P2 P123 P23 P1234 x3: y3 = P3 P234 P34 x4: y4 = P4 该表就像上面所说的那样从左向右推导。所以，每次的计算结果都可以作为下一次计算的初 始值。内维尔算法的关键是有一个父子递归关系。这种关系使用 X 点处的两个 k 阶多项式 外推的小差分直接求出相关的 K+1 阶多项式外推。这意味着，可以使用一个相对简单的公 式来利用 P1 和 P2 的值计算出 P12 的值。（如果您对推导计算感兴趣，请参阅 “Numerical Recipes” 第 3.1 节。 ） 串行 polint()函数有3 个输入参数，使用两个数组（xa和 ya）作为定义多项式的初始点集， N(n)代表数组长度，X(x)作为输入值。有两个输出参数，分别返回插值多项式(Y)在 X 点处 的近似函数值和近似误差估计(dy)。函数 vector()和 free_vector()负责为指定长度的数组 分配和释放内存。 void polint (double xa[], double ya[], int n, double x, double *y, double *dy) { int i, m, ns=1; double den,dif,dift,ho,hp,w; double *c, *d; dif = fabs(