线性方程组求解器 [Linear Equation Solver]

线性方程组求解器 [Linear Equation Solver]

求解线性方程Ax=B,其中A是N*N的矩阵,B是有N个元素的列向量。若可解则输出唯一解,或输出”multiple”,无解时输出”inconsistent”。

题目链接:
https://open.kattis.com/problems/equationsolver

如果必然有解,用行列式法很容易。行列式的定义就符合递归,非常简洁地就可以实现对行列式的求解.例如用Python用二维列表a表示行列式,则求行列式的值det(a)就可以设计成:

def det(a):
    if len(a)==1:
        return a[0][0]
    res = 0
    sign = 1
    for i in range(len(a)):
        subMat = []
        for j in range(1,len(a)):
            subLine = []
            for k in range(len(a)):
                if k!=i:
                    subLine.append(a[j][k])
            subMat.append(subLine)
        res += sign*a[0][i]*det(subMat)
        sign = -sign
    return res

但是此题可能出现行列式为0的情况,当行列式为0时,可能无解,可能有多解,这就需要处理矩阵A,在矩阵变换过程中,若出现某一行全为0,若此时这一行对应的增广矩阵的对应元素不是0,则必然无解,对应元素是0则是多解(如果有多行都是全0,只要有一行对应增广矩阵不是0就是无解的)

几乎是无可避免的要处理矩阵A,按照如下方法:将A通过交换与行变换形成上三角矩阵(在此过程中可以判断是否是多解或无解,只有唯一解才能化成三角矩阵),然后将矩阵单位化。

main()函数可以写为:

#include <cstdio>
#include <cmath>
#define MAXN 107
double Mat[MAXN*MAXN];
double BMat[MAXN];
int main(void)
{
    int N;
    char zeroColumnFlag;
    while(1)
    {
        scanf("%d",&N);
        if(N==0)
            break;
        readMat(N);
        zeroColumnFlag = solveMat(N);
        if(zeroColumnFlag)
        {
            InvalidStatus(N);
        }
        else
        {
            identityMatrix(N);
            outputResult(N);
        }
    }
    return 0;
}

然后逐步把各函数按照数学方法设计完成即可。

void readMat(int N)
{
    int i;
    for(i=0;i<N*N;i++)
    {
        scanf("%lf",Mat+i);
    }
    for(i=0;i<N;i++)
    {
        scanf("%lf",BMat+i);
    }
}
void outputMat(int N)
{
    int i,j;
    for(i=0;i<N;i++)
    {
        for(j=0;j<N;j++)
        {
            printf("%9.3f",Mat[i*N+j]);
        }
        printf("\t%9.3f\n",BMat[i]);
    }
    printf("\n");
}
int isZero(double x)
{
    const double epsilon = 0.00000001;
    if(fabs(x)<epsilon)
        return 1;
    else
        return 0;
}
void swapElem(double *a,int i,int j)
{
    double tmp = a[i];
    a[i]=a[j];
    a[j]=tmp;
}
void swapLine(int N,int k1,int k2)
{
    int i;
    for(i=0;i<N;i++)
    {
        swapElem(Mat,k1*N+i,k2*N+i);
    }
}
void upperTriangularMat(int i,int N)
{
    int k1,k2;
    for (k1=i+1;k1<N;k1++)
    {
        double rate = Mat[k1*N+i]/Mat[i*N+i];
        BMat[k1] -= rate *BMat[i];
        for(k2=i;k2<N;k2++)
        {
            Mat[k1*N+k2] -= rate*Mat[i*N+k2];
        }
    }
}
void identityMatrix(int N)
{
    int i,j;
    for ( i=N-1;i>=0;i--)
    {
        double rate = 1/Mat[i*N+i];
        BMat[i] = rate*BMat[i];
        Mat[i*N+i] = 1;
        for (j=0;j<i;j++)
        {
            BMat[j] -= BMat[i]*Mat[j*N+i];
            Mat[j*N+i]=0;
        }
    }
}
void InvalidStatus(int N)
{
    int i,j;
    double s,*p;
    for (i=0;i<N;i++)
    {
        s=0;
        p = Mat+i*N;
        for(j=0;j<N;j++)
        {
            s+=p[j]*p[j];
        }
        if(isZero(s) && !isZero(BMat[i]))
        {
            printf("inconsistent\n");
            return ;
        }
    }
    printf("multiple\n");
    return ;
}
int solveMat(int N)
{
    int i,k1;
    int zeroColumnFlag=0;
    for(i=0;i<N;i++)
    {
        zeroColumnFlag = 1;
        for (k1 = i;k1<N;k1++)
        {
            if(!isZero(Mat[k1*N+i]))
            {
                swapLine(N,k1,i);
                swapElem(BMat,k1,i);
                zeroColumnFlag = 0;
                break;
            }
        }
        if(zeroColumnFlag)
        {
            break;
        }
        else
        {
            upperTriangularMat(i,N);
        }
    }
    return zeroColumnFlag;
}
void outputResult(int N)
{
    for (int i=0;i<N;i++)
    {
        if(i>0)
            printf(" ");
        printf("%.5lf",BMat[i]);
    }
    printf("\n");
}

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注