700字范文,内容丰富有趣,生活中的好帮手!
700字范文 > MATLAB数值分析学习笔记:线性代数方程组的求解和高斯消元法

MATLAB数值分析学习笔记:线性代数方程组的求解和高斯消元法

时间:2022-08-10 02:26:02

相关推荐

MATLAB数值分析学习笔记:线性代数方程组的求解和高斯消元法

工程和科学计算的许多基本方程都是建立在守恒定律的基础之上的,比如质量守恒等,在数学上,可以建立起形如[A]{x}={b}的平衡方程。其中{x}表示各个分量在平衡时的取值,它们表示系统的状态响应;右端向量{b}由无关系统性态的常数组成通常表示为外部激励。矩阵A则表示为由系统各部分相互作用或耦合关系的参数组成的系数矩阵。在工程上则意味着[相互作用][响应]=[激励]

对于单个方程,可以采用前面介绍的一些求根法加以求解,然而事实上还有一些关系式是彼此相互耦合的,比如复杂电路的基尔霍夫定律。这就需要将这些关系式表示为一个线性代数方程组。下面就此问题介绍MATLAB求解线性代数方程组的一些方法,重点介绍高斯消元法

目录

一、取逆和“左除”

取逆

左除

二、克莱姆法则

代码实现:

问题求解:

三、高斯消元法

原理:

代码实现:

问题求解:

四、三对角方程组的求解

代码实现:

问题求解:

一、取逆和“左除”

对于形如 A*x=b的线性代数方程组,MATLAB提供了两种直接求解的方法:取逆 和“左除”

取逆

MATLAB的矩阵取逆函数为inv(A)

故x=inv(A)*b

左除

x=A\b

二、克莱姆法则

学习过线性代数的同学应该对克莱姆法则并不陌生,这里不再赘述,只是讲一下在MATLAB中的实现:

问题:求解

代码实现:

A=[0.3 0.52 1;0.5 1 1.9;0.1 0.3 0.5];b=[-0.01 0.67 -0.44]';D=det(A);A(:,1)=b;x(1)=det(A)/D;A=[0.3 0.52 1;0.5 1 1.9;0.1 0.3 0.5];A(:,2)=b;x(2)=det(A)/D;A=[0.3 0.52 1;0.5 1 1.9;0.1 0.3 0.5];A(:,3)=b;x(3)=det(A)/D;x

问题求解:

>> Clemx =-14.9000 -29.5000 19.8000

三、高斯消元法

对于小型方程,使用克莱姆法则还是很方便的,但工程上往往面对的是庞大的方程组问题,这个时候仅仅用克莱姆法则是万万不够的。

高斯消元法尽管是求解联立方程组最古老的方法之一,但直到今天依然是实际应用中最为重要的方法之一。

原理:

一般方程组(A*x=b):

步骤:向前消元-->向后回代

向前消元:第一阶段是将A矩阵消为一个上三角矩阵。方法是将第一个方程左右同时乘以

然后用第二个方程去减它这样就消掉了.

类似的,反复这样进行,就可以将原方程组变换为新的“上三角方程组”。

向后回代:第二个阶段是利用“上三角方程组”解出x。

首先解出然后将其带入倒数第二个方程,求解,以此类推。

这样就成功的解出了x。

然而,不知道大家发没发现,如果在执行向前消元时遇到=0,或者很小很小,这样就遇到了一个问题——“除数不可以为0”,这就需要每一次执行的时候选出最大的数作为“主元”,将主元作为新的“除数”。

称之前的高斯法叫朴素的高斯消元法,后者为选主元的高斯消元法。

代码实现:

function x = GaussPivot(A,b)%%建立高斯消元法的实现代码%求解方程A*x=b%步骤:%(1)向前消元;%(2)向后回代。%%输入:%A=系数矩阵%b=右向量%%输出:%x=方程的解向量[m,n]=size(A);if m~=n,error("系数矩阵须为方阵");endnb=n+1;Aug=[A b];%向前消元for k=1:n-1%选主元(max函数可以传回最大值和最大值索引)[big,i] = max(abs(Aug(k:n,k)));ipr = i+k-1;if ipr~=kAug([k,ipr],:) = Aug([ipr,k],:);%交换:将绝对值大的数作主元endfor i=k+1:nfactor = Aug(i,k)/Aug(k,k);Aug(i,k:nb) = Aug(i,k:nb)-factor*Aug(k,k:nb);endend%向后回代x = zeros(n,1);%创建解向量x(n)=Aug(n,nb)/Aug(n,n);%先求出最后一个值for i=n-1:-1:1x(i)=(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);endend

问题求解:

问题:

%%高斯消元法A=[0.3 0.52 1;0.5 1 1.9;0.1 0.3 0.5];b=[-0.01 0.67 -0.44]';x = GaussPivot(A,b)

结果:

>> GaussPivot_testx =-14.9000-29.500019.8000

四、三对角方程组的求解

工程上经常遇到形如如下的方程组:

求解步骤依然是向前消元-->向后回代,只不由于A的稀疏性,再使用高斯消元法就太消耗时间了,书中专门为其设计了代码

代码实现:

function x = Tridiag(e,f,g,r)%求解三对角方程组% 由于三对角方程的系数矩阵稀疏,运算量与n成正比,而不是高斯消元的n^3%%输入:%e =下对角线%f = 主对角线%g =上对角线%r = 右向量%%输出:%x=方程的解向量n=length(f);for k=2:nfactor=e(k)/f(k-1);f(k)=f(k)-factor*g(k-1);r(k)=r(k)-factor*r(k-1);endx(n)=r(n)/f(n);for k= n-1:-1:1x(k)=(r(k)-g(k)*x(k+1))/f(k);end

问题求解:

就上述问题进行求解

e=[0,-1,-1,-1];f=[2.04,2.04,2.04,2.04];g=[-1,-1,-1,0];r=[40.8,0.8,0.8,200.8];x = Tridiag(e,f,g,r)

>> Tridiag_testx =65.9698 93.7785 124.5382 159.4795

声明:文章来源于笔者学习【美】Steven C. CHapra所著,林赐译 《工程于科学数值方法的MATLAB实现》(第4版)的笔记,如有谬误或想深入了解,请翻阅原书。

本内容不代表本网观点和政治立场,如有侵犯你的权益请联系我们处理。
网友评论
网友评论仅供其表达个人看法,并不表明网站立场。