Demo entry 6658583

NumericalAnalysis1

   

Submitted by anonymous on Nov 08, 2017 at 02:23
Language: C. Code size: 10.3 kB.

#include <stdio.h>
#include <math.h>

#define n 501       //矩阵A的维数
#define s 2         //上半带宽
#define r 2         //下半带宽
#define error 1E-12 //精度水平

double element[s+r+1][n] = {{0},{0}};   //压缩存储的矩阵A的带内元素
double u0[n] = {0};                     //迭代的初始向量
double b = 0.16, c = -0.064;
double yk[n] = {0};                     //迭代过程中特征向量序列
double uk[n] = {0};                     //迭代过程中特征向量序列(单位化)
double b0 = 0, b1 = 0;                  //迭代过程中特征值
double e = 1;                           //相对误差

/****************************自定义函数********************************/

max(int a, int b, int c);                 //三值中求取最大值
min(int a, int b, int c);                 //三值中求取最小值
printFunc(double x[n]);                   //控制台打印向量函数
printMatrix(double aa[r+s+1][n]);         //控制台打印矩阵函数
length(double x[n]);                      //计算向量的模
unit(double x[n]);                        //向量单位化函数
mtx_vtr2(double aa[r+s+1][n],double x[n]);//幂法迭代过程中A*y->u
vtr_vtr(double x1[n],double x2[n]);       //1*n向量与n*1向量相乘
powerMethod(double x[n], double aa[r+s+1][n],double e0);//幂法迭代
powerMethod_Translation(double x[n], double aa[r+s+1][n],double p,double e0);//带原点平移的幂法迭代
doolittle(double aa[r+s+1][n]);           //对矩阵做doolittle分解(大型带状矩阵,用压缩矩阵C)
doolittle_func(double aa[r+s+1][n],double b[n]);//doolittle分解后求解方程(大型带状矩阵,用压缩矩阵C)
inPowerMethod(double x[n], double aa[r+s+1][n],double e0);//反幂法迭代
inPowerMethod_Translation(double x[n], double aa[r+s+1][n],double p,double e0);//带原点平移的反幂法迭代
det(double aa[r+s+1][n]);//计算A的行列式
init(void);//初始化函数

/*****************************主函数**********************************/
int main()
{
    double tz1=0,tz501=0,tzs=0,bm1=0,bm2=0;
    double tzi[39]={0}; //与dk最接近的特值
    double d;
    double conda; //条件数
    double deta; //行列式
    int k=0;
    init();
    bm1=powerMethod(u0,element,error); //幂法求出按模最大特征值
    tzs=inPowerMethod(u0,element,error); //反幂法求出按模最小特征值
    bm2=inPowerMethod_Translation(u0,element,-bm1,error); //带原点平移的反幂法求出最大/最小特征值

    //判断bm1,bm2大小
    if(bm1>bm2)
    {
        tz1=bm2;
        tz501=bm1;
    }
    else
    {
        tz1=bm1;
        tz501=bm2;
    }

    conda=fabs(bm1/tzs); //A的条件数
    deta=det(element); //A的行列式

    //输出最小、最大、按模最小特征值
    printf("tz1 = %12.12e\n",tz1);
    printf("tz501 = %12.12e\n",tz501);
    printf("tzs = %12.12e\n",tzs);
    printf("\n");

    //输出条件数
    printf("cond A = %12.12e\n",conda);
    printf("det A = %12.12e\n",deta);
    printf("\n");

    //计算与dk最接近的特征值并输出
    for(k=1;k<=39;k++)
    {
        d=tz501-tz1;
        d=tz1+k*d/40;
        tzi[k-1]=inPowerMethod_Translation(u0,element,d,error);
        printf("tzi%d = %12.12e\n",k,tzi[k-1]);
    }

    return 0;
}



/**************************自定义函数详细说明************************/

int max(int a, int b, int c)                 //三值中求取最大值
{
    int m;
    if(a>=b) m = a;
    else m = b;
    if(c>m) m = c;
    return m;
}

int min(int a, int b, int c)                 //三值中求取最小值
{
    int m;
    if(a<=b) m = a;
    else m = b;
    if(c<m) m = c;
    return m;
}

void printFunc(double x[n])               //控制台打印向量函数
{
    int i;
    for(i=0; i<n; i++)
    {
        printf("%12.12e\t",x[i]);
    }
    printf("\n");
}

void printMatrix(double aa[r+s+1][n])     //控制台打印矩阵函数
{
    int i,j;
    for(i=0; i<r+s+1; i++)
    {
        for(j=0; j<n; j++)
        {
            printf("%f\t",aa[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

double length(double x[n])                //计算向量的模
{
    double leng = 0;
    int i;
    for(i=0; i<n; i++)
    {
        leng = leng + pow(x[i],2);
    }
    leng = sqrt(leng);
    return leng;
}

void unit(double x[n])                    //向量单位化函数
{
    double lengthOfx;
    int i;
    for(i=0; i<n; i++)
    {
        yk[i]x[i]/lengthOfx;
    }
}

void mtx_vtr2(double aa[r+s+1][n],double x[n])//幂法迭代过程中A*y->u
{
    double mv[n] = {0};
    int i;
    mv[0] = aa[s][0]*x[0]+b*x[1]+c*x[2];
    mv[1] = aa[s][1]*x[1]+b*(x[0]+x[2])+c*x[3];

    for(i=2; i<n-2; i++)
    {
        mv[i] = aa[s][i]*x[i]+b*(x[i-1]+x[i+1])+c*(x[i-2]+x[i+2]);
    }

    mv[n-2] = aa[s][n-2]*x[1]+b*(x[n-3]+x[n-1])+c*x[n-4];
    mv[n-1] = aa[s][n-1]*[n-1]+b*x[n-2]+c*x[n-3];

    for(i=0; i<n; i++)
    {
        uk[i] = mv[i];
    }
}

double vtr_vtr(double x1[n],double x2[n])       //1*n向量与n*1向量相乘
{
    double vMultiV = 0;
    int i;
    for(i=0; i<n; i++)
    {
        vMultiV = vMultiV + x1[i]*x2[i];
    }
    return vMultiV;
}

double powerMethod(double x[n], double aa[r+s+1][n],double e0)//幂法迭代
{
    b0 = 0, b1 = 0, e = 1;
    int i,k = 0;
    for(i=0; i<n; i++)//初始化向量->u0
    {
        uk[i] = x[i];
    }
    do
    {
        k++;
        unit(uk);//单位化uk->yk
        mtx_vtr2(aa,yk);//迭代A*yk->u(k+1)
        b1 = vtr_vtr(yk,uk);//求特征向量b(k+1)=(yk,u(k+1))
        if(k>1)
        {
            e = (b1-b0)/b0;
            e = fabs(e);
        }
        b0 = b1;
    }while(e>e0);

    return b1;
}

double powerMethod_Translation(double x[n], double aa[r+s+1][n],double p,double e0)//带原点平移的幂法迭代
{
    b0 = 0,b1 = 0, e = 1;
    double apy[r+s+1][n] = {{0},{0}};
    double b2 = 0;
    int i,j,k = 0;

    for(i=0; i<n; i++)
    {
        uk[i] = x[i];
    }

    do
    {
        k++;
        unit(uk);//单位化uk->yk

        //主对角线元素-p,其余不变,存入apy
        for(i=0; i<r+s+1; i++)
        {
            if(i==s)
            {
                for(j=0; j<n; j++)
                {
                    apy[i][j] = aa[i][j]-p;
                }
            }
            else
            {
                for(j=0; j<n; j++)
                {
                    apy[i][j] = aa[i][j];
                }
            }
        }
        mtx_vtr2(apy,yk);//迭代A*yk->u(k+1)

        b1 = vtr_vtr(yk,uk);//求特征向量b(k+1)=(yk,uk(k+1))

        if(k>1)//判断是否满足精度要求
        {
            e = (b1-b0)/b0;//计算误差
            e = fabs(e);
        }
        b0 = b1

    }while(e>e0);
    b2 = b1 + p;

    return b2;
}

void doolittle(double aa[r+s+1][n])//对矩阵做doolittle分解(大型带状矩阵,用压缩矩阵C)
{
    double y[n] = {0};
    double x[n] = {0};
    int i,j,k,t=0;

    //Doolittle分解
    for(k=0; k<n; k++)
    {
        for(j=k; j<min(k+s+1,n,n);j++)//计算u第k行
        {
            for(t=max(0,k-r,j-s); t<k; t++)
            {
                aa[k-j+s][j] = aa[k-j+s][j]-aa[k-t+s][t]*aa[t-j+s][j];
            }
        }
        for(i=k+1; i<min(k+s+1,n,n);i++)//计算l第k行
        {
            for(t=max(0,k-r,i-s); t<k; t++)
            {
                aa[i-k+s][k] = aa[i-k+s][k]-aa[i-t+s][t]*aa[t-k+s][k];
            }
            aa[i-k+s][k] = aa[i-k+s][k]/aa[s][k];
        }
    }

}

void doolittle_func(double aa[r+s+1][n],double b[n])//doolittle分解后求解方程(大型带状矩阵,用压缩矩阵C)
{
    double y[n] = {0};
    double x[n] = {0};
    int i,j,k,t = 0;

    //解方程Ly=b
    y[0]=b[0];
    for(i=1;i<n;i++)
    {
        y[i]=b[i];
        for(t=max(0,0,i-r);t<i;t++)
        {
            y[i]=y[i]-aa[i-t+s][t]*y[t];
        }
    }

    //解方程Ux=y
    x[n-1]=y[n-1]/aa[s][n-1];
    for(i=n-2;i>=0;i--)
    {
        x[i]=y[i];
        for(t=i+1;t<min(i+s+1,n,n);t++)
        {
            x[i]=x[i]-aa[i-t+s][t]*x[t];
        }
        x[i]=x[i]/aa[s][i];
    }
    for(i=0;i<n;i++)
    {
        uk[i]=x[i];
    }
}

double inPowerMethod(double x[n], double aa[r+s+1][n],double e0)//反幂法迭代
{
    int i,j=0,k=0;
    double lu[r+s+1][n]={{0},{0}};
    b0=0,b1=0,e=1;

    for(j=0;j<n;j++) //初始向量->u0
    {
        uk[j]=x[j];
    }

    for(i=0;i<r+s+1;i++) //读取方程系数矩阵(A)
    {
        for(j=0;j<n;j++)
        {
            lu[i][j]=aa[i][j];
        }
    }
    doolittle(lu);//Doolittle分解
    do
    {
        k++;
        unit(uk); //单位化uk->yk

        doolittle_func(lu,yk); //解A*u(k+1)=yk->u(k+1)
        b1=1/vtr_vtr(yk,uk); //求特征向量1/b(k+1)=(yk,u(k+1))
        if(k>1) //判断是否满足精度要求
        {
            e=(b1-b0)/b0; //计算误差
            e=fabs(e);
        }
        b0=b1;
    }while(e>e0);

    return b1;
}

double inPowerMethod_Translation(double x[n], double aa[r+s+1][n],double p,double e0)//带原点平移的反幂法迭代
{
    int i,j=0,k=0;
    double b2=0;
    double lu[r+s+1][n]={{0},{0}};
    b0=0,b1=0,e=1;

    for(j=0;j<n;j++) //初始向量->u0
    {
        uk[j]=x[j];
    }

    //主对角线元素-p,其余不变,存入lu
    for(i=0;i<r+s+1;i++)
    {
        if(i==s)
        {
            for(j=0;j<n;j++)
            {
                lu[i][j]=aa[i][j]-p;
            }
        }
        else
        {
            for(j=0;j<n;j++)
            {
                lu[i][j]=aa[i][j];
            }
        }
    }

    doolittle(lu); //Doolittle分解

    do
    {
        k++;
        unit(uk); //单位化uk->yk

        doolittle_func(lu,yk); //解A*u(k+1)=yk->u(k+1)
        b1=1/vtr_vtr(yk,uk); //求特征向量1/b(k+1)=(yk,u(k+1))
        if(k>1) //判断是否满足精度要求
        {
            e=(b1-b0)/b0; //计算误差
            e=fabs(e);
        }
        b0=b1;
    }while(e>e0);
    b2=b1+p;

    return b2;
}

double det(double aa[r+s+1][n])//计算A的行列式
{
    double lu[r+s+1][n]={{0},{0}};
    double detaa=1;
    int i=0,j=0;

    for(i=0;i<r+s+1;i++)
    {
        if(i==s)
        {
            for(j=0;j<n;j++)
            {
                lu[i][j]=aa[i][j];
            }
        }
        else
        {
            for(j=0;j<n;j++)
            {
                lu[i][j]=aa[i][j];
            }
        }
    }

    doolittle(lu);

    for(i=0;i<n;i++)
    {
        detaa=detaa*lu[s][i];
    }

    return detaa;
}

void init(void) //初始化函数
{
    int i=1,j=0;

    //大型带状矩阵压缩
    for(i=1;i<=n;i++)
    {
        element[s][i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);
    }
    for(i=2;i<n;i++)
    {
        element[0][i]=c;
        element[4][n-1-i]=c;
    }
    for(i=1;i<n;i++)
    {
    element[1][i]=b;
    element[3][n-1-i]=b;
    }

    //选择初始向量u0
    for(i=0;i<n;i++)
    {
        if(i%2==0)
        {
            u0[i]=1;
        }
        else
        {
            u0[i]=0;
        }
    }
}

This snippet took 0.03 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).