一、前言

参考知乎文章:LU分解、LUP分解、Cholesky分解中的LUP分解，LUP分解相比LU分解更加稳定，在矩阵求逆、求解解线性方程组有广泛运用。笔者修复了文章所述LUP算法一个BUG，当待分解的矩阵某一列均为0时，例如矩阵[0 1;0 1]，无法进行LUP分解。增加基于LUP分解求解线性方程组代码、在矩阵方程AX=B中，当B为单位阵时，利用LUP分解可以快速稳定地求解矩阵的逆。

二、LUP分解MATLAB代码

function [L, U, P] = lup(A)
if size(A, 1) ~= size(A, 2)error('矩阵行数与列数不等')
endU = A;
L = eye(size(A));
P = eye(size(A));
n = size(A, 2);
for j = 1 : n - 1%Select i(>=j) that maximizes |U(i,j)|index = -1;maxValue = 0.0;for i = j : ntemp = abs(U(i,j));if temp > maxValuemaxValue = temp;index = i;endendif index == -1continue;end%Interchange rows of U: U(j, j : n) <-> U(i, j : n)for k = j : ntemp = U(j, k);U(j, k) = U(index, k);U(index, k) = temp;end%Interchange rows of L: L(j, 1 : j - 1) <-> L(i, 1 : j - 1)for k = 1 : j - 1temp = L(j, k);L(j, k) = L(index, k);L(index, k) = temp;end%Interchange rows of P: P(j, 1 : n) <-> P(i, 1 : n)for k = 1 : ntemp = P(j, k);P(j, k) = P(index, k);P(index, k) = temp;endfor i = j + 1 : nL(i, j) = U(i, j) / U(j, j);for k = j : nU(i, k) = U(i, k) - L(i, j) * U(j, k);endend
end
end

三、LUP分解求解方程组MATLAB代码

function X = solve_matrix_equation_by_lup(A, B)
if size(A, 1) ~= size(A, 2)error('矩阵行数与列数不等')
endU = A;
L = eye(size(A));
n = size(A, 2);
m = size(B, 2);
for j = 1 : n - 1%Select i(>=j) that maximizes |U(i,j)|index = -1;maxValue = 0.0;for i = j : ntemp = abs(U(i,j));if temp > maxValuemaxValue = temp;index = i;endendif index == -1continue;end%Interchange rows of U: U(j, j : n) <-> U(i, j : n)for k = j : ntemp = U(j, k);U(j, k) = U(index, k);U(index, k) = temp;end%Interchange rows of L: L(j, 1 : j - 1) <-> L(i, 1 : j - 1)for k = 1 : j - 1temp = L(j, k);L(j, k) = L(index, k);L(index, k) = temp;end%Interchange rows of P: P(j, 1 : n) <-> P(i, 1 : n), C = P * B，等价于对B交换行for k = 1 : mtemp = B(j, k);B(j, k) = B(index, k);B(index, k) = temp;endfor i = j + 1 : nL(i, j) = U(i, j) / U(j, j);for k = j : nU(i, k) = U(i, k) - L(i, j) * U(j, k);endend
endfor i = 1 : nif abs(U(i, i)) < 1.0e-20disp('方程无解')return;end
end%L * y = C
y = zeros(n, m);
for j = 1 : mfor i = 1 : nsum = 0.0;for k = 1 : i - 1sum = sum + L(i, k) * y(k, j);endy(i, j) = B(i, j) - sum;end
end%U * x = y
X = zeros(n, m);
for j = 1 : mfor i = n :-1 : 1sum = 0.0;for k = i + 1 : nsum = sum + U(i, k) * X(k, j);endX(i, j) = (y(i, j) - sum) / U(i, i);end
endend

四、LUP分解、矩阵求逆、求解线性方程组测试例子

clc;
clear;
A = [1 4 5;2 3 2;3 4 2];
B = [1 0 0;0 1 0;0 0 1];[l,u,p] = lu(A)
[L,U,P] = lup(A)invA = inv(A)
x = A\B
X = solve_matrix_equation_by_lup(A, B)

LU分解、矩阵求逆与解线性方程组(matlab代码)相关推荐

ADMM,ISTA,FISTA算法步骤详解，MATLAB代码，求解LASSO优化问题
ADMM,ISTA,FISTA算法步骤详解,MATLAB代码,求解LASSO优化问题原创文章!转载需注明来源:©️ Sylvan Ding's Blog ❤️ 实验目的了解 ADMM, ISTA, ...
紧凑存储的杜利特尔分解法Doolittle(LU分解法)_解线性方程组的直接解法
紧凑存储的杜利特尔分解法Doolittle(LU分解法)_解线性方程组的直接解法标签:计算方法实验 /* 紧凑存储的杜利特尔分解法Doolittle:如果初始矩阵不要求保留的话,可以紧凑存储.因为每 ...
【机械臂优化】基于粒子群算法实现考虑关节限位约束下的冗余机械臂求逆解附Matlab代码）
1 简介 2 部分代码 %%%%%%%%%%%%%%%%%%采用PSO算法对运动学冗余机械臂求一组最优逆解%%%%%%%%%%%%%%%%%%% %该程序对一个具有四自由度的机械臂做位置控制,由操作空 ...
Dijkstra算法和Floyd算法详解（MATLAB代码）
一.Dijkstra算法 1.算法简介 Dijkstra算法是由E.W.Dijkstra于1959年提出,又叫迪杰斯特拉算法,它应用了贪心算法模式,是目前公认的最好的求解最短路径的方法.算法解决的是有 ...
小波变换对图像的分解与重构（含matlab代码）
01 小波变换原理所谓的小波的小是针对傅里叶波而 ...
数值计算——雅可比迭代法解线性方程组(附代码)
1.雅克比迭代法的计算过程: (1).取初始向量: (1) (2 ...
gammatone 滤波器详解及其MATLAB代码实现
一.GammaTone 滤波器详解定义: 外界语音信号进入耳蜗的基底膜后,将依据频率进行分解并产生行波震动,从而刺激听觉感受细胞[1].GammaTone 滤波器是一组用来模拟耳蜗频率分解特点的滤波 ...
鲁棒局部均值分解（RLMD）附Matlab代码
✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信.
蚁群算法原理详解和matlab代码
1原理: 蚂蚁在寻找食物源的时候,能在其走过的路径上释放一种叫信息素的激素,使一定范围内的其他蚂蚁能够察觉到.当一些路径上通过的蚂蚁越来越多时,信息素也就越来越多,蚂蚁们选择这条路径的概率也就越高,结 ...
lle算法的matlab实现,lle算法详解及matlab代码实现
LLE算法代码 % LLE ALGORITHM (using K nearest neighbors) % % [Y] = lle(X,K,dmax) % % X = data as D x N ma ...

LU分解、矩阵求逆与解线性方程组(matlab代码)

文章目录

一、前言

二、LUP分解MATLAB代码

三、LUP分解求解方程组MATLAB代码

四、LUP分解、矩阵求逆、求解线性方程组测试例子

LU分解、矩阵求逆与解线性方程组(matlab代码)相关推荐

最新文章

热门文章