系列文章目錄
綜合實例應(yīng)用:方程組的求解
前言
無論工程應(yīng)用問題,還是數(shù)學(xué)計算問題,方程組都是解決問題轉(zhuǎn)化的重要途徑之一,將復(fù)雜問題轉(zhuǎn)化為簡單的方程組矩陣求解問題。
一、求解四元一次線性方程組
>> %創(chuàng)建方程組系數(shù)矩陣
>> A=[2 1 -5 1;1 -3 0 -6;0 2 -1 2;1 4 -7 6];
>> b=[8 9 -5 0]';
>> %判斷方程是否有解
>> %求方程組的秩
>> r=rank(4)
r =
1
>> B=[A,b];%創(chuàng)建增廣矩陣
>> s=rank(B)
s =
4
>> %r=s=n(未知數(shù))=4,則該齊次線性方程組有唯一解。
>> %利用矩陣的逆
>> x0=pinv(A)*b
x0 =
3.0000
-4.0000
-1.0000
1.0000
二、利用矩陣分解求解
利用矩陣分解來求解線性方程組,是工程計算中最常用的計算。
1.LU分解法
LU分解法是先將系數(shù)矩陣A進行LU分解,得到LU=PA,然后解Ly=Pb,最后再解Ux=y得到原方程組的解。
編寫利用LU分解法求解線性方程組Ax=b的自定義函數(shù)M文件,操作方法:
函數(shù)solvebyLU的程序,如下所示
function x=solvebyLU(A,b)
%該函數(shù)利用LU分解法求線性方程組Ax=b的解
flag=isexist(A,b); %調(diào)用自定義函數(shù)isexist()判斷方程組解的情況
if flag==0
disp('該方程組無解!');
x=[];
return;
else
r=rank(A);
[m,n]=size(A);
[L,U,P]=lu(A);
b=P*b;
%解Ly=b
y(1)=b(1);
if m>1
for i=2:m
y(i)=b(i)-L(i,1:i-1)*y(1:i-1)';
end
end
y=y';
%解Ux=y得原方程組得一個特解
x0(r)=y(r)/U(r,r);
if r>1
for i=r-1:-1:1
x0(i)=(y(i)-U(i,i+1:r)*x0(i+1:r)')/U(i,i);
end
end
x0=x0';
if flag==1 %若方程組有唯一解
x=x0;
return;
else %若方程組有無窮多解
format rat;
Z=null(A,'r'); %求出對應(yīng)齊次方程組的基礎(chǔ)解系
[mZ,nZ]=size(Z);
x0(r+1:n)=0;
for i=1:nZ
t=sym(char([107 48+i]));
k(i)=t; %取k=[k1,k2...,];
end
x=x0;
for i=1:nZ
x=x+k(i)*Z(:,i); %將方程組的通解表示為特解加對應(yīng)齊次通解形式
end
end
end
函數(shù)isexist()的程序,如下所示
function y=isexist(A,b)
%該函數(shù)用來判斷線性方程組Ax=b的解的存在性
%若方程組無解則返回0,若有唯一解則返回1,若有無窮多解則返回Inf。
[m,n]=size(A);
[mb,nb]=size(b);
if m~=mb
error('輸入有誤!');
return;
end
r=rank(A);
s=rank([A,b]);
if r==s &&r==n
y=1;
elseif r==s&&r<n
y=Inf;
else
y=0;
end
命令行代碼,如下所示
>> A=[2 1 -5 1;1 -3 0 -6;0 2 -1 2;1 4 -7 6];
>> b=[8 9 -5 0]';
>> x2=solvebyLU(A,b)
x2 =
3
-4
-1
1
2.QR分解法
利用QR分解法先將系數(shù)矩陣A進行QR分解A=QR,然后解Qy=b,最后解Rx=y得到原方程組的解
1.編寫求解線性方程組Ax=b的函數(shù)solvebyQR,代碼如下:文章來源:http://www.zghlxwxcb.cn/news/detail-802911.html
function x=solvebyQR(A,b)
%該函數(shù)利用QR分解法求線性方程組Ax=b的解
flag=isexist(A,b); %調(diào)用自定義函數(shù)isexist()
if flag==0
disp('方程組無解');
x=[];
return;
else
r=rank(A);
[m,n]=size(A);
[Q,R]=qr(A);
b=Q'*b;
%解Rx=b得原方程組得一個特解
x0(r)=b(r)/R(r,r);
if r>1
for i=r-1:-1:1
x0(i)=(b(i)-R(i,i+1:r)*x0(i+1:r)')/R(i,i);
end
end
x0=x0';
if flag==1 %若方程組有唯一解
x=x0;
return;
else %若方程組有無窮多解
format rat;
Z=null(A,'r'); %求出對應(yīng)齊次方程組得基礎(chǔ)解系
[mZ,nZ]=size(Z);
x0(r+1:n)=0;
for i=1:nZ
t=sym(char([107 48+i]));
k(i)=t; %取k=[k1,...,kr];
end
x=x0;
for i=1:nZ
x=x+k(i)*Z(:,i); %將方程組的通解表示為特解加對應(yīng)齊次通解形式
end
end
end
總結(jié)
綜合實例—方程組的求解,到這里就結(jié)束啦!感謝觀看,希望這篇文章對大家有幫助。文章來源地址http://www.zghlxwxcb.cn/news/detail-802911.html
到了這里,關(guān)于MATLAB:方程組的求解的文章就介紹完了。如果您還想了解更多內(nèi)容,請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章,希望大家以后多多支持TOY模板網(wǎng)!