解线性方程组的直接方法:LU分解法及其C语言算法实现
在上一篇博客里面,笔者介绍了解线性方程组的列主元Guass消元法,这篇将介绍LU分解法及其算法实现.
什么是LU分解?
对于一个线性方程组Ax=b,其中A是非奇异系数矩阵,b是线性方程组右端项,在列主元Guass消元法里面我们知道,最后的系数矩阵A将变成一个上三角矩阵,并且是通过一系列的行变换而来的,设最后得到的上三角矩阵为U,结合高等代数的知识,一个矩阵左乘一个初等矩阵,相当于进行一次行变换,因此设每一次A左乘的初等矩阵为Li(i=1,2,…,n),则有LnL(n-1)…L1A=U,由于Ln均为初等矩阵,且均为下三角单位矩阵(因为每次A进行消元所做的行变换均是从上面的行消去下面的行),所以设L=LnL(n-1)*…*L1,LA=U,A=L-1U,L-1也为单位下三角矩阵,然后得到了A=LU,我们可以设Ux=y,Ly=b,由于L为下三角矩阵,求解难度较小,因此通过求解y向量,再由Ux=y求解x,这便是LU分解全部步骤了,对于LU分解矩阵的详细计算过程,大家可以参考这个网站link
接下来话不多说,上代码
初始化矩阵
double** init_Matrix(int r, int c)
{double** p = new double* [r];int d = c + 1;for (int i = 0; i < r; i++){p[i] = new double[d];memset(p[i], 0, sizeof(double) * d);}cout << "请输入线性方程组对应的增广矩阵:" << endl;for (int i = 0; i < r; i++){for (int j = 0; j < d; j++){cin >> p[i][j];}}return p;
}
进行LU分解
void LU_Position(double**a,int r ,int c)
{double num1 = 0,num2=0;for (int i=0;i<r;i++){if (i==0)//第一行不做处理,单独算第一列{for (int j = 1; j < c; j++){a[j][i] = a[j][i] / a[0][0];}}else{for (int j = i; j < c; j++){num1 = 0;for (int k = 0; k < i; k++){num1 += a[i][k] * a[k][j];}a[i][j] = a[i][j] - num1;}for (int j = i+1; j < r; j++){num2 = 0;for (int k = 0; k < i; k++){num2 += a[j][k] * a[k][i];}a[j][i] = (a[j][i] - num2) / a[i][i];}}}cout << "所得的LU矩阵为:" << endl;for (int k = 0; k < r; k++){for (int n = 0; n < c; n++){printf("%f\t", a[k][n]);}cout << endl;}
}
注意:此时我们得到了LU,由于这两个矩阵均为稀疏矩阵,且拼接后大小与A矩阵相同,因此我们一直接将计算的两个矩阵储存在A矩阵中(A中最后一列为右端项,不进行处理),这部分计算大家特别要注意的是数组的行列下边关系
计算y向量及x向量
void calculate(double*y, double*x, double**a,int r)
{y[0] = a[0][r];double sum=0;for (int i = 1; i < r; i++){sum = 0;for (int k = 0; k < i; k++){sum += a[i][k] * y[k];}y[i] = a[i][r] - sum;}cout << "所求的向量Y为:" << endl;for (int i = 0; i < r; i++){printf("%f\t", y[i]);}cout << endl;for (int i = r-1; i >=0; i--){sum = 0;for (int j=i+1;j<r;j++){sum += a[i][j] * x[j];}x[i] = (y[i] - sum) / a[i][i];}cout << "所求线性方程组的解向量X为:" << endl;for (int i = 0; i < r; i++){printf("%f\t", x[i]);}cout << endl;
}
这里我们将得到的y向量储存在数组y中
程序完整代码
#include<iostream>
#include<Windows.h>
using namespace std;
/*
测试数据2 2
2 3 5
1 -1 03 3
3 2 -3 -2
1 1 1 6
1 2 -1 23 3
1 2 -3 1
2 -1 3 5
3 -2 2 14 4
4 -3 6 7 11
1 1 3 4 10
-2 9 -7 1 10
3 3 -4 20 255 5
28 -3 0 0 0 10
-3 38 -10 0 -5 0
0 -10 25 -15 0 0
0 0 -15 45 0 0
0 -5 0 0 30 0
*/
double** init_Matrix(int r, int c)
{double** p = new double* [r];int d = c + 1;for (int i = 0; i < r; i++){p[i] = new double[d];memset(p[i], 0, sizeof(double) * d);}cout << "请输入线性方程组对应的增广矩阵:" << endl;for (int i = 0; i < r; i++){for (int j = 0; j < d; j++){cin >> p[i][j];}}return p;
}
//直接用A来储存LU和方程组右边的常数项 进行LU分解时 右边常数项不做处理 即最后一列全部不处理
void LU_Position(double**a,int r ,int c)
{double num1 = 0,num2=0;for (int i=0;i<r;i++){if (i==0){for (int j = 1; j < c; j++){a[j][i] = a[j][i] / a[0][0];}}else{for (int j = i; j < c; j++){num1 = 0;for (int k = 0; k < i; k++){num1 += a[i][k] * a[k][j];}a[i][j] = a[i][j] - num1;}for (int j = i+1; j < r; j++){num2 = 0;for (int k = 0; k < i; k++){num2 += a[j][k] * a[k][i];}a[j][i] = (a[j][i] - num2) / a[i][i];}}}cout << "所得的LU矩阵为:" << endl;for (int k = 0; k < r; k++){for (int n = 0; n < c; n++){printf("%f\t", a[k][n]);}cout << endl;}
}
void calculate(double*y, double*x, double**a,int r)
{y[0] = a[0][r];double sum=0;for (int i = 1; i < r; i++){sum = 0;for (int k = 0; k < i; k++){sum += a[i][k] * y[k];}y[i] = a[i][r] - sum;}cout << "所求的向量Y为:" << endl;for (int i = 0; i < r; i++){printf("%f\t", y[i]);}cout << endl;for (int i = r-1; i >=0; i--){sum = 0;for (int j=i+1;j<r;j++){sum += a[i][j] * x[j];}x[i] = (y[i] - sum) / a[i][i];}cout << "所求线性方程组的解向量X为:" << endl;for (int i = 0; i < r; i++){printf("%f\t", x[i]);}cout << endl;
}
void LU_position_main() {cout << "输入矩阵的行列:" << endl;int i = 0, j = 0;cin >> i >> j;double** p = init_Matrix(i, j);LU_Position(p, i, j);double* a = new double[i];memset(a, 0, sizeof(double) * i);double* b = new double[i];memset(b, 0, sizeof(double) * i);calculate(a, b, p, i);delete[]a;delete[]b;for (int i = 0; i < j; i++){delete[]p[i];}delete[]p;
}
int main(void) {LU_position_main();system("pause");return 0;
}
欢迎交流探讨~~
解线性方程组的直接方法:LU分解法及其C语言算法实现相关推荐
- 数值分析3-解线性方程组的高斯消去法、LU分解法及列主元消去法的matlab程序和调试方法
对于形如Ax=b的线性方程组,在线性代数中是通过求逆的方式求解的,即x=A-1b,而在数值分析中,解线性方程组的方法是通过直接法或者迭代法来实现的,今天写的三个程序都属于直接法,分别为高斯消去法.LU ...
- matlab lu解线性方程,MATLAB报告用LU分解法求解线性方程组.doc
MATLAB报告用LU分解法求解线性方程组 <MATLAB>期末大报告 报告内容:用LU分解法求解线性方程组 学院(系): 专 业: 班 级: 学 号: 学生姓名: 2014 年 6 月 ...
- 紧凑存储的杜利特尔分解法Doolittle(LU分解法)_解线性方程组的直接解法
紧凑存储的杜利特尔分解法Doolittle(LU分解法)_解线性方程组的直接解法 标签:计算方法实验 /* 紧凑存储的杜利特尔分解法Doolittle:如果初始矩阵不要求保留的话,可以紧凑存储.因为每 ...
- 计算方法:列主元消去法,LU分解法, 雅可比迭代法,高斯塞德尔迭代法 解线性方程(C++)
Matrix.h包括矩阵类Matrix的定义,Matrix.cpp包括该类成员函数的实现,LinearEqu.h包括线性方程类LinearEqu的定义,继承自Matrix类,其中solve()方法为列 ...
- Matlab | Lab4——用LU 分解法、 Jacobi 迭代、 Gauss-Seidel 迭代 解线性病态方程组(系数矩阵为Hilbert矩阵)
1.要求 考虑线性方程组Hx=b,其中H为n阶Hilbert矩阵,即 通过先给定解(例如取x的各个分量为1),再计算出右端向量b的办法给出一个精确解已知的问题. (1)分别编写Doolittle LU ...
- 二、解线性方程组的直接方法
https://zhuanlan.zhihu.com/p/30485749 设 n n n阶线性方程组: { a 11 x 1 + a 12 x 2 + . . . + a 1 n x n = b 1 ...
- Doolittle分解法(LU分解法)的Python实现
在解一般的非奇异矩阵线性方程组的时候,或者在迭代改善算法中,需要使用LU分解法. 对于一个一般的非奇异矩阵A=(a11, a12,-,a1n,a21,-ann),可分解为一个下三角矩阵L和一个上三角矩 ...
- Python02 雅克比迭代法 Gauss-Seidel迭代法 列选主元法 LU分解法(附代码)
1. 实验结果 (1)在定义的矩阵类中设置需要求解的方程为: (2)在 test.py 中选择雅克比迭代法求解: 输入:最大容许迭代次数和精度要求: 输出:根据谱半径判断方法是否收敛,收敛时得到满足精 ...
- matlab解方程组方法,第二章解线性方程组的直接方法matlab用法
第二章解线性方程组的直接方法matlab用法 第二章 解线性方程组的直接方法的 MATLAB 程序24. 在这章中我们要学习线性方程组的直接法,特别是适 合用数学软件在计算机上求解的方法. 2.1 方 ...
最新文章
- Spring StateMachine,教你快速实现一个状态机
- Golang map输出排序
- 脱了马甲我也认识你: 聊聊 Android 中类的真实形态
- JAVA中的适配器应用_Java适配器模式应用之电源适配器功能详解
- Playfab开发(一)如何调用PlayFab接口
- flask中关于endpoint端点、url_map映射、view_func视图函数,view_functions、及视图函数名是否何以相同的问题?
- SpringBoot集成MyBatis-Plus自定义SQL
- ipc原理linux,传统的Linux中IPC通信原理
- python引入模块教程_python导入模块--案例
- jQuery入门笔记
- 如何用计算机制作思维导向图,电脑怎样制作思维导图,手把手教你绘制思维导图简单方法...
- android 客户端和服务端cookie,如何在Android客户端注入及清除Cookie教程
- 我国传统-基本能力常识
- 常见的数字证书格式与格式转换
- ARM开发初级-Windows环境下的STM32开发环境搭建(包含missing compiler version 5的解决方法)-学习笔记02
- 【程序员的自我修养】[动态图文] 超详解函数栈帧
- Linux的拓扑结构,linux底下的makefile框架拓扑结构分析
- 高仿富途牛牛-组件化(一)-支持页签拖拽、增删、小工具
- 基于NIOS II的1553B总线开发板
- 批量用title的内容命名html文件,使用批处理批量复制文件并重命名
热门文章
- matlab的imshow, image, imagesc区别
- RouterOS搭建一台SSTP Server用于远程办公
- AD 2020 入门教程
- 明解C语言 7-6例题分析
- 《你好,放大器》----学习记录(四)
- 【计算机网络】-【考研复试面试】-整合
- linux找不到镜像文件,为什么启动u盘找不到镜像文件_u盘启动找不到映像文件的解决方法...
- ORACLE EBS出现In order to access this application, you must install the J2SE Plugin version 1.6.0_07
- dux修改index.php,[mcj]Dux大前端主题增加网站顶端公告模块
- 第一课:Python变量