数值计算中的搜索、插值与系数求解

# 数值计算中的搜索、插值与系数求解 ## 1. 有序表搜索方法 在计算函数 $f(x)$ 时，若采用特定插值方案（如四阶多项式插值），需从一组制表的 $x_i$ 和 $f_i$ 中快速定位 $x$ 在 $x_i$ 表中的位置。这一问题虽不属于数值分析范畴，但在实际中极为常见。 ### 1.1 二分查找法 给定一个横坐标数组 $xx[j]$（$j = 1, 2, \cdots, n$），其元素单调递增或递减，以及一个数 $x$，要找到一个整数 $j$，使 $x$ 位于 $xx[j]$ 和 $xx[j + 1]$ 之间。可定义虚拟数组元素 $xx[0]$ 和 $xx[n + 1]$ 为正负无穷，这样 $j$ 始终在 $0$ 到 $n$ 之间。二分查找法通常是较好的选择，大约经过 $\log_2 n$ 次尝试就能找到正确位置。以下是二分查找的代码实现： ```c void locate(float xx[], unsigned long n, float x, unsigned long *j) { unsigned long ju,jm,jl; int ascnd; jl=0; ju=n+1; ascnd=(xx[n] >= xx[1]); while (ju-jl > 1) { jm=(ju+jl) >> 1; if (x >= xx[jm] == ascnd) jl=jm; else ju=jm; } if (x == xx[1]) *j=1; else if(x == xx[n]) *j=n-1; else *j=jl; } ``` 若使用零偏移数组，需从 $xx$ 的地址和返回值 $j$ 中减去 $1$。 ### 1.2 相关值搜索法 当多次搜索大表，且连续搜索的横坐标几乎相同时，每次都进行完整的二分查找会很浪费。此时可采用以下方法：先从表中的一个猜测位置开始，以 $1, 2, 4$ 等增量“搜索”，直到将目标值括起来，然后在括起来的区间内进行二分查找。以下是代码实现： ```c void hunt(float xx[], unsigned long n, float x, unsigned long *jlo) { unsigned long jm,jhi,inc; int ascnd; ascnd=(xx[n] >= xx[1]); if (*jlo <= 0 || *jlo > n) { *jlo=0; jhi=n+1; } else { inc=1; if (x >= xx[*jlo] == ascnd) { if (*jlo == n) return; jhi=(*jlo)+1; while (x >= xx[jhi] == ascnd) { *jlo=jhi; inc += inc; jhi=(*jlo)+inc; if (jhi > n) { jhi=n+1; break; } } } else { if (*jlo == 1) { *jlo=0; return; } jhi=(*jlo)--; while (x < xx[*jlo] == ascnd) { jhi=(*jlo); inc <<= 1; if (inc >= jhi) { *jlo=0; break; } else *jlo=jhi-inc; } } } while (jhi-(*jlo) != 1) { jm=(jhi+(*jlo)) >> 1; if (x >= xx[jm] == ascnd) *jlo=jm; else jhi=jm; } if (x == xx[n]) *jlo=n-1; if (x == xx[1]) *jlo=1; } ``` ### 1.3 搜索后的处理 `locate` 和 `hunt` 例程返回一个索引 $j$，使目标值位于表项 $xx[j]$ 和 $xx[j + 1]$ 之间。若要使用 `polint` 或 `ratint` 等例程获得 $m$ 点插值，需提供长度为 $m$ 的 $xx$ 和 $yy$ 数组。可通过以下公式计算： $k = IMIN(IMAX(j - (m - 1) / 2, 1), n + 1 - m)$ 然后使用偏移 $k$ 的数组地址调用插值例程，如 `polint(&xx[k - 1], &yy[k - 1], m, ...)`。 ## 2. 插值多项式的系数 有时需要知道通过少量点的插值多项式的系数，而非其值。例如，计算函数及其多个导数的同时插值值，或对制表函数的一段与其他函数进行卷积。但需注意，插值多项式的系数通常不如其在所需横坐标处的值准确，因此仅为计算插值值而确定系数并非明智之举。 ### 2.1 求解系数的线性方程 若插值多项式表示为 $y = c_0 + c_1x + c_2x^2 + \cdots + c_Nx^N$，则系数 $c_i$ 需满足线性方程： $\begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^N \\ 1 & x_1 & x_1^2 & \cdots & x_1^N \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_N & x_N^2 & \cdots & x_N^N \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_N \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_N \end{bmatrix}$ 这是一个范德蒙德矩阵，可使用特殊方法求解，比一般线性方程求解技术更高效。 ### 2.2 求解系数的代码实现 以下是求解系数的代码： ```c #include "nrutil.h" void polcoe(float x[], float y[], int n, float cof[]) { int k,j,i; float phi,ff,b,*s; s=vector(0,n); for (i=0;i<=n;i++) s[i]=cof[i]=0.0; s[n] = -x[0]; for (i=1;i<=n;i++) { for (j=n-i;j<=n-1;j++) s[j] -= x[i]*s[j+1]; s[n] -= x[i]; } for (j=0;j<=n;j++) { phi=n+1; for (k=n;k>=1;k--) phi=k*s[k]+x[j]*phi; ff=y[j]/phi; b=1.0; for (k=n;k>=0;k--) { cof[k] += b*ff; b=s[k]+x[j]*b; } } free_vector(s,0,n); } ``` ### 2.3 另一种求解方法 另一种技术是利用已有的函数值插值例程 `polint`。先插值或外推求插值多项式在 $x = 0$ 处的值，该值即为 $c_0$。然后从 $y_i$ 中减去 $c_0$ 并除以相应的 $x_i$，去掉一个点后重复此过程求 $c_1$ 等。以下是代码实现： ```c #include <math.h> #include "nrutil.h" void polcof(float xa[], float ya[], int n, float cof[]) { void polint(float xa[], float ya[], int n, float x, float *y, float *dy); int k,j,i; float xmin,dy,*x,*y; x=vector(0,n); y=vector(0,n); for (j=0;j<=n;j++) { x[j]=xa[j]; y[j]=ya[j]; } for (j=0;j<=n;j++) { polint(x-1,y-1,n+1-j,0.0,&cof[j],&dy); xmin=1.0e38; k = -1; for (i=0;i<=n-j;i++) { if (fabs(x[i]) < xmin) { xmin=fabs(x[i]); k=i; } if (x[i]) y[i]=(y[i]-cof[j])/x[i]; } for (i=k+1;i<=n-j;i++) { y[i-1]=y[i]; x[i-1]=x[i]; } } free_vector(y,0,n); free_vector(x,0,n); } ``` ### 2.4 注意事项 若 $x = 0$ 不在制表 $x_i$ 的范围内或离其较远，插值多项式的系数通常会变得很大，导致病态问题。此外，对光滑函数进行过高阶插值时，插值多项式可能会出现振荡。 ## 3. 二维或多维插值 ### 3.1 二维插值的基本概念 在二维插值中，给定一个函数值矩阵 $ya[1..m][1..n]$，以及数组 $x1a[1..m]$ 和 $x2a[1..n]$，满足 $ya[j][k] = y(x1a[j], x2a[k])$。要估计未制表点 $(x_1, x_2)$ 处的函数值 $y$。重要的概念是 $