范德蒙德和Teoplitz方程组的解法

计算数学与科学工程计算研究所 陆嵩

简单介绍

工程中的很多实际问题的处理，比如说图像处理的某些情况，最后往往归结为比较容易处理的Vandermonde方程组和Teoplitz方程组的求解问题，因此对于这类特殊问题的快速求解方法就显得尤为重要。

所谓Vandermonde方程组，指的是如下方程组：

Vz=b

其中V是Vandermonde矩阵，

V=V(x0,…,xn)=⎛⎝⎜⎜⎜⎜1x0⋮xn01x1⋮xn1⋯⋯⋯1xn⋮xnn⎞⎠⎟⎟⎟⎟

而Toeplitz方程组指的是：

Tnx=b

其中Tn是Toeplitz矩阵，

Tn=⎛⎝⎜⎜⎜⎜⎜⎜1γ1⋮γn−1γ11⋮⋯⋯⋯⋱γ1γn−1⋮γ11⎞⎠⎟⎟⎟⎟⎟⎟

一般矩阵的求解有各种分解，那么对于这类特殊的方程组，有么有什么特殊的解法呢？下面做一个简单的介绍。

基本思想

Vandermonde方程组

假设我们已经得到了范德蒙德逆转置的的一个分解V−T=UL，那么V−1=LTUT，可以得到方程组的解为z=V−1b=LTUTb。
所以我们主要把力量放在的到范德蒙德逆转置的一个分解。

首先考虑这样的一个方程组求解:

VTv=y

这里

v=(a0,a1,…an)T

.

注意这里的VT是带转置的，并不是原来的Vandermonde方程组的求解问题了。因为范德蒙德矩阵的特殊性，不难想到，这其实是以v为稀系数的多项式在(xi,yi)点的插值问题，如果用ci表示牛顿插值表达式的系数，其实就是k阶差商，我们可以通过差商之间的递推关系，找到ci
和yi的一个关系，进而，通过将牛顿插值多项式以某种巧妙的方式展开，找到ci和插值多项式系数ai之间的关系，换言之，就是找到了v和y之间的一个关系。这个过程比较繁琐，但是细细品味，也是很有意思的。这里考虑到期末考试将至，就不再细述了。给出一个结果，同过一步步推导，最后可以得到v和y之间的形如下面这样的一个关系：

v=ULy

结合前式，可以得到V−Ty=ULy，有y的任意性，我们就得到了V−T的一个形如V−T=UL的分解，这就得到了我们想要的。

Teoplitz方程组

在说Teoplitz方程组前，可以先提一下Yule-Walker方程组，Yule-Walker方程组指的是Teoplitz方程组b=−(γ1,…,γn−1,γn)T（γn*）的特殊情况。简单地说Yule-Walker方程组的求解，是通过矩阵分块的方法，找到yk和yk+1之间的关系，进而可以从一阶解y1出发，递推得到方程组的解，这一个过程所需的运算量为3n2/2。类似于Yule-Walker方程组的求解过程，我们也可以递推地得到一般Teoplitz方程组的解，由于时间关系，就不再细述。

程序和结果

下面直接给出程序和简单的算例和结果，我这里用的环境是VS2013：

Vandermonde方程组

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 6

int main(int argc, char** argv) {
    double b[N] = { 1, 2, 3, 4, 5, 6 };
    double x[N] = { 5, 6, 7, 8, 9, 10 };
    int i, k, k0;
    int n = N - 1;

    printf("要求解的方程组为：\n x = \n");
    for (i = 0; i < N; i++)
    {
        printf("%f ", x[i]);

    }

    printf("\n b = \n");
    for (i = 0; i < N; i++)
    {
        printf("%f\n ", b[i]);

    }


    for (k = 0; k <= n - 1; k++)
    {
        for (i = n; i > k; i--)
        {
            b[i] = b[i] - x[k] * b[i - 1];
        }
    }

    for (k0 = k; k >= 0; k--)
    {
        for (i = k + 1; i <= n; i++)
        {
            b[i] = b[i] / (x[i] - x[i - k - 1]);

        }
        for (i = k; i <= n - 1; i++)
        {
            b[i] = b[i] - b[i + 1];
        }
    }
    printf("\nthe solution is :\n");
    for (i = 0; i <= n; i++)
    {
        printf("%f \n", b[i]);
    }
    getchar();

    return 0;
}

Toeplitz方程组

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 3

int main(int argc, char** argv)

{
    double gamma[N - 1] = { 0.3, 0.6 };
    double b[N] = { 5, 6, 7 };
    int n = N;
    int i, k;
    double s = 0;
    double zeta[N], xi[N], delta, alpha, mu, v[N], sigma[N];
    zeta[0] = -gamma[0];
    xi[0] = b[0];
    delta = 1.0;
    alpha = -gamma[0];
    k = 1;

    printf("gamma （from 1 to %d） is\n", N - 1);
    for (i = 0; i <= N - 2; i++)
    {
        printf("%f ", gamma[i]);
    }
    printf("\n\nb  is \n");
    for (i = 0; i <= N - 1; i++)
    {
        printf("%f ", b[i]);
    }

    while (1)
        //while(k<n-1)
    {
        delta = (1.0 - alpha*alpha)*delta;
        for (i = 1; i <= k; i++) {
            s = s + gamma[i - 1] * xi[k - i];
        }
        mu = (b[k] - s) / delta;
        s = 0.0;
        for (i = 1; i <= k; i++) {
            v[i - 1] = xi[i - 1] + mu*zeta[k - i];
        }
        for (i = 1; i <= k; i++) {
            xi[i - 1] = v[i - 1];
        }
        xi[k] = mu;
        if (k >= n - 1) break;
        for (i = 1; i <= k; i++)
        {
            s = s + zeta[k - i] * gamma[i - 1];
        }
        alpha = -(gamma[k] + s) / delta;
        s = 0.0;

        for (i = 1; i <= k; i++) {
            sigma[i - 1] = zeta[i - 1] + alpha*zeta[k - i];
        }

        for (i = 1; i <= k; i++) {
            zeta[i - 1] = sigma[i - 1];

        }
        zeta[k] = alpha;
        k = k + 1;
    }

    printf("\n\nsolution:\n");
    for (i = 0; i <= n - 1; i++) {
        printf("%lf\n", xi[i]);
    }
    getchar();
    return 0;
}