1. 程式人生 > >採用GAUSS列主消元法求解線性方程組(MATLAB)

採用GAUSS列主消元法求解線性方程組(MATLAB)

%%求解任意線性方程組的解
clc;
clear all;
format long e
disp('線性方程組求解,請輸入引數');
n=input('維數n=');
A=input('矩陣A=');
b=input('右端項b=');
eps=input('控制精度eps=');


b=b';   %%變為列向量
A=[A b];    %%矩陣增廣


for k=1:n-1
    B=A(k:n,k); %%先將第k列可能作為主元的元素取出方至矩陣B
    P=max(abs(B));  %%選主元P
    if(P<eps)   %%控制小主元
        disp('無解');
        break;
    else
        u=find((abs(B))==P); %%計算主元所在行相對與k行的位置
            if(u~=1)
                A([k,u],:)=A([u,k],:);  %%換行
            end
            m=A(k+1:n,k)/A(k,k);    %%求出各行行乘數並放至矩陣m
            for i=1:length(m)
            A(k+i,k:n+1)=A(k+i,k:n+1)-m(i)*A(k,k:n+1);  %%消元按行進行
            end
    end
end
if A(n,n)==0
    disp('無解'); %%若矩陣A不滿秩,則無解
else
    x(n)=A(n,n+1)/A(n,n);   %%由最後一行首先求出方程組的第一個解x(n)
    for i=n-1:-1:1  %%計算第i個解x(i)
        for j=1:1:n-i  %%利用回代思想
            A(i,n+1)=A(i,n+1)-A(i,i+j)*x(i+j);  %%減去已知部分
        end
        x(i)=A(i,n+1)/A(i,i);
    end
end
disp('方程組的解');
x=x'   %%輸出方程組的解