python 切比雪夫插值 切比雪夫迭代法

方程组求解的切比雪夫半迭代加速方法

背景介绍

解方程组的迭代算法有Jacobi迭代，SOR方法等，但是对于一般的矩阵，这类算法不一定收敛，即使收敛，也有可能收敛得很慢。所以我们试图找到一个方法，来加速迭代算法的收敛速度。

基本思想

考虑迭代方法如下，$G$为迭代矩阵 $x_{k+1} = Gx_k+c$
这里计算$x_{k+1}$，只用到了前一个值，我们设想，能不能综合利用前面已知的所有的值的信息，使得收敛的速度更快呢？我们考虑前k+1个值的某种加权平均来作为$x_{k+1}$迭代的初值。即令 $\widetilde{x}_k = \sum\limits_{j = 1}^k {a_j^k}x^j$
这里的$a_j^k$作为权重，和要为1。事实上，就是通过这种方式，在原迭代序列的基础上，构造了一个新的迭代序列，用新的迭代序列，去替代旧的迭代序列。假设前k+1步的值都是用前所述一般的迭代产生的，那么对k+1步的值使用加权平均修正后的误差为
$u_k = (\sum\limits_{j=0}^k a_j^kG^j)u_0\triangleq q_k(G)u_0$
为了使误差减小得尽可能快，我们应该选择合适的$a_k^j$使得$q_k(G)$的谱半径尽可能小。换言之，我们应该选择合适的系数，使得多项式$q_k(t)$在$G$的所有的特征值处取值后的模的最大值达到最小。因为G的特征值是未知的一些离散的点，所以这看起来并不是一件容易的事情。我们可以在包含G所有特征值的凸区域上考虑这个问题，问题就变成了一个连续的问题，即找一个多项式，使得其在这一块区域上的最大值达到最小。咋一看，这其实就是一个求解最小零偏差多项式的问题，最小零偏差多项式实际上就是寻求函数$f(t)=t^k$的k-1阶最佳逼近多项式，当然，这是一个$L_inf$意义下的最佳逼近问题。数值逼近方法，对于这类问题给我我们一些很好的结果。切比雪夫多项式，在某种意义下，就是最小零偏差多项式……

给出如下的结果:

• 当G的特征值为实数都小于1，那么上述的连续最佳一致逼近问题有一个和切比雪夫多项式有关的解。
• 当G的特征值都为实数，但有一个大于1（由相容性保证其不等于1），那么我们可以构造新的迭代矩阵，将问题变为第一种情况。
• 当G的特征值不全是实数，但是包含在某个椭圆里，那么在复平面内它对应的最佳逼近问题，也有一个与切比雪夫多项式有关的解。

程序和结果

这里只给出第一种情况的一个程序和示例，示例具体参数可看程序。

C代码:

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

void main()
{
    //函数声明
    double*  numTimesVec(double a, double *b, double *c);
    double* vectorOperator(double *a, double *b, double *sum);
    void iterativeProcess(double* x1, double* x0);
    int assignment(double* a, double* b);
    int convergenceOrNot(double* a, double* b);


    //变量声明
    double rho1, rho2, xi, nu;
    double x1[N], x2[N];
    double alpha = -2;
    double beta = 0.5;
    double x0[N] = { 1, 1, 1 };

    //给定初始值
    rho1 = 2.0;

    //计算
    iterativeProcess(x1, x0);
    xi = (2 - beta - alpha) / (beta - alpha);
    nu = 2.0 / (2 - beta - alpha);
    while (convergenceOrNot(x0, x1) == 0)
    {
        rho2 = 1.0 / (1 - rho1 / (4 * xi*xi));
        double x_temp[N];
        iterativeProcess(x_temp, x1);
        double temp1[N];
        numTimesVec(nu, x_temp, temp1);//c语言不能返回一个数组，要想返回，最好传入一个空数组用于存放
        int i;
        /*  printf("temp1=\n");
            for (i = 0; i < N; i++)
            printf("%f   ", temp1[i]);*/


        double temp2[N];
        numTimesVec(1.0 - nu, x1, temp2);

        /*printf("temp2=\n");
        for (i = 0; i < N; i++)
        printf("%f   ", temp2[i]);*/


        double temp3[N];
        numTimesVec((1 - rho2), x0, temp3);

        /*  printf("temp3=\n");
            for (i = 0; i < N; i++)
            printf("%f   ", temp2[i]);*/


        //vectorOperator(temp1, temp2);
        double temp4[N];
        vectorOperator(temp1, temp2, temp4);
        double temp5[N];
        numTimesVec(rho2, temp4, temp5);

        vectorOperator(temp5, temp3, x2);
        rho1 = rho2;
        assignment(x0, x1);
        assignment(x1, x2);
    }
    int i;
    printf("迭代矩阵和过程见代码，迭代结果为：\n");
    for (i = 0; i < N; i++)
    {
        printf("%f\n", x1[i]);
    }

    getchar();

}




void iterativeProcess(double* x1, double* x0)
{
    int i, j;
    double G[N][N] = { { -1.0 / 2, 0.0, 0.0 },
    { 0.0, 0.0, 0.0 },
    { 0.0, 0.0, 1.0 / 2 }
    };
    double c[N] = { 1, 1, 1 };
    //  double x1[N] = { 0, 0, 0, 0, 0 };

    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {
            c[i] = c[i] + G[i][j] * x0[j];
        }
        x1[i] = c[i];
    }

    //return x1;
}
int assignment(double* a, double* b)
{
    int i;
    for (i = 0; i < N; i++)
    {
        a[i] = b[i];
    }
    return 0;

}
int convergenceOrNot(double* a, double* b)
{
    int i;
    double s = 0;
    double c[N];
    for (i = 0; i < N; i++)
    {
        c[i] = a[i] - b[i];
        s = s + c[i] * c[i];
    }
    s = sqrt(s / N);
    if (s < epsilon)

        return 1;

    else
        return 0;

}
double* vectorOperator(double *a, double *b, double *sum)
{
    int i;
    for (i = 0; i < N; i++)
    {
        sum[i] = a[i] + b[i];
    }
    return sum;
}
double*  numTimesVec(double a, double *b, double *c)
{
    int i;
    //  static double c[N];
    for (i = 0; i < N; i++)
    {
        c[i] = a*b[i];
    }
    return c;
}

MATLAB代码：

clc
clear
beta = 0.5;
alpha = -2;
epsilon = 0.001;
x0 = [1 1 1]';
rho1 = 2;
c = [1 1 1]';
G = diag([-1/2,0,1/2]);
x1 = G*x0 + c;
xi = (2-beta-alpha)/(beta-alpha);
nu = 2/(2-beta-alpha);



while(norm((x0-x1),2)>epsilon)
    rho2=1/(1-(rho1/4/xi^2));
    x2 = rho2*(nu.*(G*x1+c)+(1-nu).*x1)+(1-rho2)*x0;
    rho1 = rho2;
    x0 = x1;
    x1 = x2;
end
x1