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

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

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

简单介绍

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

所谓Vandermonde方程组,指的是如下方程组:

Vz=b
其中V是Vandermonde矩阵,
V=V(x0,,xn)=1x0xn01x1xn11xnxnn
而Toeplitz方程组指的是:
Tnx=b

其中Tn是Toeplitz矩阵,

Tn=1γ1γn1γ11γ1γn1γ11

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

基本思想

Vandermonde方程组

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

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

VTv=y
这里
v=(a0,a1,an)T
.

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

v=ULy

结合前式,可以得到VTy=ULy,有y的任意性,我们就得到了VT的一个形如VT=UL的分解,这就得到了我们想要的。

Teoplitz方程组

在说Teoplitz方程组前,可以先提一下Yule-Walker方程组,Yule-Walker方程组指的是Teoplitz方程组b=(γ1,,γn1,γn)Tγn***)的特殊情况。简单地说Yule-Walker方程组的求解,是通过矩阵分块的方法,找到ykyk+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;
}

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

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;
}

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

猜你喜欢