Lu分解法的C语言实现

将系数矩阵A转变成等价两个矩阵L和U的乘积
,其中L和U分别是单位下三角矩阵和上三角矩阵。当A的所有顺序主子式都不为0时,矩阵A可以分解为A=LU(所有顺序主子式不为0,矩阵不一定不可以进行LU分解)。其中L是下三角矩阵,U是上三角矩阵。

LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵,其变换矩阵就是一个单位下三角矩阵。这正是所谓的杜尔里特算法(Doolittle algorithm):从下至上地对矩阵A做初等行变换,将对角线左下方的元素变成零,然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵,这一系列单位下三角矩阵的乘积的逆就是L矩阵,它也是一个单位下三角矩阵。

这里写图片描述

这里写图片描述

这里写图片描述

构造矩阵A,b

                {1,2,3,4, }
    A=[         {1,4,2,-8,}
                {1,-1,4,1}   ]
                {1,3,5,2}

b = [14,-17,2,8]^T

代码实现:

#include <stdio.h>
#include <stdlib.h>
//LU分解法实现解线性方程组
//copyright @ Mryang

double sumU(double L[4][4] ,double U[4][4], int i, int j ){
    double sU = 0.0;
    for (int k = 1; k <= i-1 ; k++)
    {
        sU += L[i-1][k-1] * U[k-1][j-1];
    }

    return sU;
}//计算求和1

double sumL(double L[4][4] ,double U[4][4], int i, int j ){
    double sL = 0.0;
    for (int k = 0; k <= j-1; k++)
    {
        sL += L[i-1][k-1] * U[k-1][j-1];
    }

    return sL;

}//计算求和2

double sumY(double L[4][4] ,double y[4],int i){
    double sY=0.0;
    for (int k = 1; k <= i - 1; k++)
    {
        sY += L[i-1][k-1] * y[k-1];
    }
    return sY;
}//计算求和3

double sumX(double U[4][4] ,double x[4],int i ,int m){
    double sX = 0.0;
    for (int k = i+1; k <= m; k++)
    {
        sX += U[i-1][k-1] * x[k-1];
    }
    return sX;
}//计算求和4
int main(){

    double a[4][4] = { {1,2,3,1,},
                    {1,4,1,-1,},
                    {1,-1,-2,3,},
                    {1,3,-1,2}};//将系数存入二维数组
    double L[4][4] = {0};
    double U[4][4] = {0};//初始化部分
    double b[4] = {8,8,12,19};
    int n = 4;//n阶
    //输出[Ab]
    printf("[A]:\n");
    for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= n; j++)
        {
            printf("%f\t", a[i-1][j-1]);
        }
        printf("\n");
    }

    //计算L,U
    for (int i = 1; i <= n; i++)
    {
        L[i-1][i-1] = 1;//对角线元素为1
        for (int j = i; j <= n; j++)
        {
            //由于数组下标从0开始 所以i-1,j-1
            U[i-1][j-1] = a[i-1][j-1] - sumU(L,U,i,j);

            if(j+1 <= n) L[j][i-1] = (a[j][i-1] - sumL(L,U,j+1,i))/U[i-1][i-1];//i变j+1,j变i 
        }

    }





    //输出U
    printf("U:\n");
    for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= n; j++)
        {
            printf("%f\t",U[i-1][j-1]);
        }
        printf("\n");
    }
    //输出L
    printf("L:\n");
    for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= n; j++)
        {
            printf("%f\t",L[i-1][j-1]);
        }
        printf("\n");
    }

    //由Ly=b 求y
    double y[4] = {0.0};
    y[0] = b[0];//y(1) = b(1);

    for (int i = 2; i <= n; i++)
    {
        y[i-1] = b[i-1] - sumY(L,y,i);
    }

    //由 Ux=y 求x
    double x[4] = {0.0};

    for (int i = n; i >= 1; i--)
    {
        x[i-1] =( y[i-1] - sumX(U,x,i,n))/ U[i-1][i-1];
    }



    //输出y
    printf("y:\n");
    for (int i = 0; i < n; i++)
    {
        printf("%f\n",y[i]);
    }
    printf("\n");
    //输出x
    printf("x:\n");
    for (int i = 0; i < n; i++)
    {
        printf("%f\n",x[i]);
    }
    printf("\n");
    system("pause");
    return 0;
}

运行结果:

这里写图片描述

所以
x1=1
x2= 2
x3=-1
x4 = 3

  • 15
    点赞
  • 67
    收藏
    觉得还不错? 一键收藏
  • 7
    评论
### 回答1: 稀疏矩阵的LU分解是一种将稀疏矩阵分解成下三角矩阵L和上三角矩阵U的方LU分解可以在线性方程组、计算逆矩阵和计算行列式等方面起到重要作用。 在C语言中,可以通过以下步骤实现稀疏矩阵的LU分解。 1. 定义稀疏矩阵的数据结构 可以使用二维数组或链表等数据结构来表示稀疏矩阵,其中只存储非零元素的值和位置。 2. 实现稀疏矩阵的LU分解函数 首先,通过高斯消元将原始矩阵转化为上三角矩阵,同时记录下每一步的消元因子。然后,通过回代计算下三角矩阵。 3. 编写主程序 在主程序中读取输入的稀疏矩阵,并调用LU分解函数进行计算。最后,输出得到的下三角矩阵L和上三角矩阵U。 总结起来,稀疏矩阵的LU分解C语言中的实现需要定义合适的数据结构来表示稀疏矩阵,并编写相应的分解函数。该函数首先通过高斯消元将矩阵转化为上三角矩阵,再通过回代计算下三角矩阵。最后,主程序读取输入并调用分解函数,输出分解后的矩阵。 ### 回答2: 稀疏矩阵LU分解是一种用于稀疏矩阵的分解。在LU分解中,矩阵A被分解为一个下三角矩阵L和一个上三角矩阵U的乘积,即A = LU。 在C语言中,我们可以使用数组来表示稀疏矩阵。假设稀疏矩阵A的维度为n*n,并且我们已将稀疏矩阵A存储在一个名为matrix的二维数组中。我们可以使用以下代码来实现稀疏矩阵LU分解C语言函数: ```c #include <stdio.h> void sparseLU(int n, int matrix[n][n]) { int i, j, k; for (i = 0; i < n; i++) { for (k = i + 1; k < n; k++) { if (matrix[i][i] == 0) { printf("Error: Division by zero\n"); return; } matrix[k][i] = matrix[k][i] / matrix[i][i]; for (j = i + 1; j < n; j++) { matrix[k][j] = matrix[k][j] - matrix[k][i] * matrix[i][j]; } } } } int main() { int n = 3; // 稀疏矩阵的维度 int matrix[3][3] = {{2, 1, 1}, {4, 3, 3}, {8, 7, 9}}; sparseLU(n, matrix); // 输出LU分解后的矩阵 printf("L:\n"); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { if (i > j) { printf("%d ", matrix[i][j]); } else if (i == j) { printf("1 "); } else { printf("0 "); } } printf("\n"); } printf("U:\n"); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { if (i <= j) { printf("%d ", matrix[i][j]); } else { printf("0 "); } } printf("\n"); } return 0; } ``` 在以上代码中,我们使用两个for循环来遍历稀疏矩阵A,并进行LU分解。在每次迭代中,我们计算矩阵U的第i行和第k行之间的系数,并将其存储在matrix[k][i]中。然后,我们用这些系数更新矩阵A中第k行的元素。 最后,我们在main函数中打印出分解后的矩阵L和U。其中,矩阵L中的上三角元素都被设置为0,对角线元素为1;矩阵U中的下三角元素都被设置为0。 以上是一种用C语言实现稀疏矩阵LU分解的方。当然,具体的实现方式还需根据具体情况来进行调整。 ### 回答3: 稀疏矩阵LU分解是一种将稀疏矩阵分解为两个矩阵的方,其中一个矩阵是单位下三角矩阵,另一个是上三角矩阵,用于快速决稀疏线性方程组的问题。以下是用C语言实现稀疏矩阵LU分解的思路: 1. 首先定义一个结构体用于表示稀疏矩阵,包括矩阵的行列数和非零元素个数等属性,以及一个三元组的数组用于存储矩阵中非零元素的位置和值。 2. 创建函数lu_decomposition,该函数参数包括稀疏矩阵和LU分解后的矩阵。 3. 在函数内部,首先根据稀疏矩阵中非零元素的位置和值构建稀疏矩阵的三元组表示。 4. 然后从左上角开始遍历稀疏矩阵,按列将第一个非零元素作为主元。 5. 根据主元列位置,在主元下方的行中找到主元下面的非零元素,并进行消元操作,计算出相应的乘因子,将其保存在上三角矩阵中。 6. 在主元所在的列中,将主元上方的非零元素进行消元操作,并将乘因子保存在单位下三角矩阵中。 7. 遍历完所有的主元后,得到稀疏矩阵的LU分解结果。 8. 最后,将LU分解后的矩阵返回。 通过以上步骤,可以实现C语言中对稀疏矩阵进行LU分解的功能。在实际应用中,可以根据具体需对该算进行优化和改进,以提高分解效率和准确性。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 7
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值