Gauss列主元消去法函数

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了Gauss列主元消去法函数相关的知识,希望对你有一定的参考价值。

 

 

 

%列主元消去法解方程组Ax=b,实现PA=LU

function [x,detA] =gauss(A,b)

n=length(b);[p,q]=size(A);

if p~=q||p~=n      fprintf(‘方阵的维数不同,请重新输!); %检错

end 

%为提高运行速度,给L,U,x,c,d1赋初值

L=zeros(n,n);

U=zeros(n,n);

x=zeros(n,1);

c=zeros(1,n);

d1=0;

%按列选主元,并进行行交换,记录行信息

for i=1:n-1     

    max=abs(A(i,i));

    m=i;

    for j=i+1:n         

        if max<abs(A(j,i))

            max=abs(A(j,i));

            m=j;        

        end

    end     

    if (m~=i)         

        for k=i:n             

            c(k)=A(i,k);                        

            A(i,k)=A(m,k);             

            A(m,k)=c(k);                     

        end

        d1=b(i);

        b(i)=b(m);

        b(m)=d1;

    end     

  %进行消元计算      

  for k=i+1:m                          

      for j=i+1:n              

          A(k,j)=A(i,j)-A(i,j)*A(i,i);        

      end

     b(k)=b(k)-b(i)*A(k,i)/A(i,i);

     A(k,i)=0;

  end

  %回代求解  

  x(n)=b(n)/A(n,n);

  for i=n-1:-1:1     

      sum=0.0;    

      for j=i+1:n          

          sum=sum+A(i,j)*x(j);    

      end

      x(i)=(b(i)-sum)/A(i,i);

  end

  %计算行列式的值

  detA=1;

  for k=1:n

      detA=detA*A(k,k);

  end

  %输出PA=LU中的L,U的信息

  for i=1:n     

      for j=1:n          

          if i<j                       

              U(i,j)=A(i,j);         

          elseif i==j            

              L(i,j)=1;            

              U(i,j)=A(k,j);         

          else

              L(i,j)=A(i,j);                    

          end

      end

  end

end

 

本题方程组求解(脚本文件):

A=[-0.002 2 2

    1 0.78125 0

    3.996 5.5625 4];

b=[0.4;1.3816;7.4178];

%[L,U,x,detA]=myLU(A,b)

[x,detA]=gauss(A,b)

disp(‘列主元高斯消去法得到的解为:)

 

六.运行结果:

 x =

    1.8167

    0.0527

   -0.0337

detA =

   798.0666

 

以上是关于Gauss列主元消去法函数的主要内容,如果未能解决你的问题,请参考以下文章

高斯消去法与矩阵三角分解法(LU分解)

《数值分析》-- 高斯消去法与矩阵三角分解法(LU分解)

线性方程组的分解法——列主元消去法

高斯消去追赶法 matlab

Matlab数值微分

用列主元消去法分别解方程组Ax=b,用MATLAB程序实现(最有效版)