数值分析3-解线性方程组的高斯消去法、LU分解法及列主元消去法的matlab程序和调试方法

对于形如Ax=b的线性方程组，在线性代数中是通过求逆的方式求解的，即x=A^-1b,而在数值分析中，解线性方程组的方法是通过直接法或者迭代法来实现的，今天写的三个程序都属于直接法，分别为高斯消去法、LU分解法和列主元消去法。
还有全主元消去，平方根法，追赶法，其中平方根法，追赶法的思路都高斯消去法差不多，稍微修改一下程序就行，而全主元消去法相较于列主元消去法而言，判断大小和换行换列会更麻烦一些，用循环语句能用到吐，没啥需求的话我实在是懒得写；
高斯消去法，就是线性代数中通过把系数矩阵化为行阶梯矩阵然后求解的方法，而LU分解法和高斯消去法十分相似，只不过是回代过程有所不同，所以我把这两个程序放在一起写了。程序中用到了较多的循环语句，要仔细分析一下，不然写的过程中很容易出错。
首先是高斯消去法：
代码块1：高斯消去法

clear;
clc;
fprintf('高斯消去法解线性方程组:\n')
n=3;
A=[1 2 3;2 5 2;3 1 5];
b=[14 18 20];                 %注意，这里我用的S是行向量，计算时注意转置
%n=input('请输入系数矩阵A的阶数:n=');
%A=input('请输入系数矩阵A:A=');
%b=input('请输入结果向量b:b=');%b为行向量，计算时需要转置
%消元计算---->产生行阶梯矩阵：
A1=A;                         %为了保留原始参数，以下用A1存储产生的上三角矩阵
b1=b;
for j=1:n                     %从第一列到第n列for i=j+1:n               %每一行从主元素的下侧第一个元素开始if A1(j,j)~=0         %每次的主元素都不可为零，这是能使用高斯消去法的前提M(i,j)=A1(i,j)/A1(j,j);for k=j:n         %第i行j列开始计算到第i行n列A1(i,k)=A1(i,k)-M(i,j)*A1(j,k);endb1(i)=b1(i)-M(i,j)*b1(j);   else error('求解出错，高斯消去法不适用于此方程组')endend
end
%回代求解：x(n)=b1(n)/A1(n,n);    %先求出解向量的最后一个元素，用于回代i=n-1;                 %从n-1行开始，回代求解
while i>0               %从第n-1行到第一行，依次回代s=0;               %对s赋初值0，for j=i+1:ns=A1(i,j)*x(j)+s;%计算：（已知解向量×系数矩阵对应元素）的和endx(i)=(b1(i)-s)/A1(i,i);%回代求解，解向量中已知元素+1.用于回代求x(i-1)i=i-1;
end
fprintf('经过消元回代，此方程组的解为:\n')
disp(x);%可以用inv(A)*b'直接求解验证计算
%inv(A)*b'

对于有的线性方程组是无法使用高斯消去法求解的，能使用高斯消去法的前提是系数矩阵A的顺序主子式不为0，这也是系数矩阵能LU分解的充分不必要条件，还有，如果系数矩阵A在计算过程中的主元素绝对值过小的的话，得到的结果误差会很大，这也是数值计算中要避免的一种情况，即：绝对值较小的数做除数；这种线性方程组可以考虑使用列主元消去法或者全主元消去法来计算。（但是matlab计算精度很高，就硬算也没啥问题）

所谓LU分解法，其实就是高斯消去法的一种变形，对于顺序主子式不为零的n阶矩阵A，都可以分解成A=L*U的形式，其中L为对角元素为1的下三角矩阵，U为上三角矩阵；因此Ax=b可以等价为：Ux=y，Ly=b，分成两步来回代求解，其代码如下所示：
代码块2：LU分解法解线性方程组

clear;
clc;
fprintf('LU分解法解线性方程组:\n')
n=3;
A=[1 2 3;2 5 2;3 1 5];
b=[14 18 20];                 %注意，这里我用的是行向量，计算时注意转置
%n=input('请输入系数矩阵A的阶数:n=');
%A=input('请输入系数矩阵A:A=');
%b=input('请输入结果向量b:b=');%b为行向量，计算时需要转置
k=1;
%消元计算---->产生上三角矩阵：
A1=A;                         %为了保留原始参数，以下用A1存储产生的上三角矩阵
%LU分解：
for j=1:n                     %从第一列到第n列for i=j+1:n               %每一行从主元素的下侧第一个元素开始if A1(j,j)~=0         %每次的主元素都不可为零，这是能使用高斯消去法的前提M(i,j)=A1(i,j)/A1(j,j);for k=j:n         %第i行j列开始计算到第i行n列A1(i,k)=A1(i,k)-M(i,j)*A1(j,k);end  else error('求解出错，系数矩阵无法LU分解')endend
end
for i=1:nM(i,i)=1;%M即为所求的下三角矩阵L，上三角矩阵U=A1
end
L=M;
U=A1;
%回代求解：
y(1)=b(1);
for i=2:ns=0;for k=1:i-1s=s+L(i,k)*y(k);end y(i)=b(i)-s;
end
x(n)=y(n)/U(n,n);
i=n-1;
while i>0s=0;for k=i+1:ns=s+U(i,k)*x(k);endx(i)=(y(i)-s)/U(i,i);i=i-1;
end
fprintf('此方程组的解为:\n')
disp(x);
%用inv(A)*b'直接求解验证计算

列主元消去法是高斯消去法的一种改进方法，它可以避免高斯消去法中主元素过小导致的误差；有利必有弊，列主元消去法计算量要大一些，因为他要有一个比较大小和重新排序换行的问题，这里我写的比较麻烦，用了好多个循环，我自己看的都要晕了，代码块如下：
代码块3：列主元消去法解线性方程组：

%列主元消去法解线性方程组
clear;
clc;
fprintf('列主元消去法解线性方程组:\n')
n=3;
%A=[1 2 3;2 5 2;3 1 5];
%b=[14 18 20];                 %注意，这里我用的S是行向量，计算时注意转置
n=input('请输入系数矩阵A的阶数:n=');
A=input('请输入系数矩阵A:A=');
b=input('请输入结果向量b:b=');%b为行向量，计算时需要转置
%消元计算---->产生上三角矩阵：
A1=A;                         %为了保留原始参数，以下用A1存储产生的上三角矩阵
b1=b;
for j=1:n                     %从第一列到第n列for i=j+1:n               %每一行从主元素的下侧第一个元素开始for k=i:n%这个循环选出第j列最大的元素做主元if abs(A1(k,j))>abs(A1(j,j))%if语句，别忘了绝对值，ps：冒泡排序，真滴好用for k1=j:n              %用来交换第k行和第j行a1=A1(j,k1);a2=A1(k,k1);A1(j,k1)=a2;A1(k,k1)=a1;endt=b1(k);b1(k)=b1(j);b1(j)=t;%别忘了结果向量，也要相应的交换行endend
%上述增补的部分就是列主元消去的精髓，下面的部分和高斯消去一样if A1(j,j)~=0         %每次的主元素都不可为零，这是能使用高斯消去法的前提M(i,j)=A1(i,j)/A1(j,j);for k=j:n         %第i行j列开始计算到第i行n列A1(i,k)=A1(i,k)-M(i,j)*A1(j,k);endb1(i)=b1(i)-M(i,j)*b1(j);   else error('求解出错，高斯消去法不适用于此方程组')endend
end
%回代求解：x(n)=b1(n)/A1(n,n);    %先求出解向量的最后一个元素，用于回代i=n-1;                 %从n-1行开始，回代求解
while i>0               %从第n-1行到第一行，依次回代s=0;               %对s赋初值0，for j=i+1:ns=A1(i,j)*x(j)+s;%计算：（已知解向量×系数矩阵对应元素）的和endx(i)=(b1(i)-s)/A1(i,i);%回代求解，解向量中已知元素+1.用于回代求x(i-1)i=i-1;
end
fprintf('经过消元回代，此方程组的解为:\n')
disp(x);%其实这样子写会计算时多判断好多次，不够简洁，

列主元消去比高斯消去增加了16-29行，就是用来判断大小和排序的环节，看着这一排的for和if，真滴让人头大呀！
代码中用s来实现累加器的功能，用for来实现变量递增的循环，用while实现变量递减的循环；
我是直接定义了矩阵A，结果向量b，也可以把我直接给出的矩阵和结果向量注释掉，把下方我写的自己输入A,b的语句放出来。
写程序的过程中谁也不可能说永远能一步到位，一下子就写出正确的来，所以出现错误后，可以按下图方式打点，运行后按步进方式，一步一步的运行并观测数据，排查错误出现在那一步；

应该没什么要注意的了，因为我测试的数据比较少，不知道哪里还有没发现的bug，欢迎大家交流呀~~~

数值分析3-解线性方程组的高斯消去法、LU分解法及列主元消去法的matlab程序和调试方法相关推荐

紧凑存储的杜利特尔分解法Doolittle(LU分解法)_解线性方程组的直接解法
紧凑存储的杜利特尔分解法Doolittle(LU分解法)_解线性方程组的直接解法标签:计算方法实验 /* 紧凑存储的杜利特尔分解法Doolittle:如果初始矩阵不要求保留的话,可以紧凑存储.因为每 ...
matlab lu解线性方程,MATLAB报告用LU分解法求解线性方程组.doc
MATLAB报告用LU分解法求解线性方程组 <MATLAB>期末大报告报告内容:用LU分解法求解线性方程组学院(系): 专业: 班级: 学号: 学生姓名: 2014 年 6 月 ...
解线性方程组的直接方法:LU分解法及其C语言算法实现
在上一篇博客里面,笔者介绍了解线性方程组的列主元Guass消元法,这篇将介绍LU分解法及其算法实现. 什么是LU分解? 对于一个线性方程组Ax=b,其中A是非奇异系数矩阵,b是线性方程组右端项,在列主 ...
Matlab | Lab4——用LU 分解法、 Jacobi 迭代、 Gauss-Seidel 迭代解线性病态方程组（系数矩阵为Hilbert矩阵）
1.要求考虑线性方程组Hx=b,其中H为n阶Hilbert矩阵,即通过先给定解(例如取x的各个分量为1),再计算出右端向量b的办法给出一个精确解已知的问题． (1)分别编写Doolittle LU ...
计算方法：列主元消去法，LU分解法，雅可比迭代法，高斯塞德尔迭代法解线性方程（C++）
Matrix.h包括矩阵类Matrix的定义,Matrix.cpp包括该类成员函数的实现,LinearEqu.h包括线性方程类LinearEqu的定义,继承自Matrix类,其中solve()方法为列 ...
Doolittle分解法（LU分解法）的Python实现
在解一般的非奇异矩阵线性方程组的时候,或者在迭代改善算法中,需要使用LU分解法. 对于一个一般的非奇异矩阵A=(a11, a12,-,a1n,a21,-ann),可分解为一个下三角矩阵L和一个上三角矩 ...
用直接分解法求方程组的C语言程序,c语言编程求解线性方程组论文
计算机编程求解线性方程组第一章绪论在自然科学.工程技术.经济和医学各领域中产生的许多实际问题都可以通过数学语言描述为数学问题,也就是说,由实际问题建立数学模型,然后应用各种数学方法和技巧来求解 ...
Python02 雅克比迭代法 Gauss-Seidel迭代法列选主元法 LU分解法(附代码)
1. 实验结果 (1)在定义的矩阵类中设置需要求解的方程为: (2)在 test.py 中选择雅克比迭代法求解: 输入:最大容许迭代次数和精度要求: 输出:根据谱半径判断方法是否收敛,收敛时得到满足精 ...
LU分解法 | matlab
% LU分解法 % M为输入的增广矩阵 % precision为输入的精度要求,如不输入或输入有误,则默认为10位if nargin == 2trydigits(precision);catchdis ...

数值分析3-解线性方程组的高斯消去法、LU分解法及列主元消去法的matlab程序和调试方法

数值分析3-解线性方程组的高斯消去法、LU分解法及列主元消去法的matlab程序和调试方法相关推荐

最新文章

热门文章