專題一MATLAB求解方程

關于專題一MATLAB求解方程第1頁，講稿共41頁，2023年5月2日，星期三基本概念稀疏矩陣(sparsematrix)：矩陣中非零元素的個數遠遠小于矩陣元素的總數），并且非零元素的分布沒有規(guī)律。與之相區(qū)別的是，如果非零元素的分布存在規(guī)律（如上三角矩陣、下三角矩陣、對稱矩陣），則稱該矩陣為特殊矩陣。稠密矩陣:非0元素占所有元素比例較大的矩陣。若n階矩陣A的行列式不為零，即|A|≠0，則稱A為非奇異矩陣，否則稱A為奇異矩陣。 把矩陣A的行換成相應的列，得到的新矩陣稱為A的轉置矩陣，記作A‘。AA‘=E（E為單位矩陣，）或A’A=E，則n階實矩陣A稱為正交矩陣。第2頁，講稿共41頁，2023年5月2日，星期三基本概念求矩陣A的秩rank(A)求矩陣A的跡trace(A)求矩陣A的行列式det(A)求矩陣V的1范數 norm(V,1)求矩陣V的2范數 norm(V)或norm(V,2)求矩陣V的∞范數 norm(V,inf)第3頁，講稿共41頁，2023年5月2日，星期三魔方矩陣魔方矩陣是每行、每列及兩條對角線上的元素和都相等的矩陣。對于n階魔方陣，其元素由1,2,3,…,n2共n2

個整數組成.magic(n)：生成一個n階魔方陣,各行各列及兩條對角線和為(1+2+3+...+n2)/n第4頁，講稿共41頁，2023年5月2日，星期三范得蒙矩陣范得蒙(Vandermonde)矩陣最后一列全為1，倒數第二列為一個指定的向量，其他各列是其后列與倒數第二列的點乘積。

vander(V)生成以向量V為基礎向量的范得蒙矩陣。例>>A=vander([1;2;3;5])A=11118421279311252551第5頁，講稿共41頁，2023年5月2日，星期三希爾伯特矩陣Hilbert矩陣的每個元素hij=1/(i+j-1)hilb(n)生成n階希爾伯特矩陣

invhilb(n)求n階的希爾伯特矩陣的逆矩陣注意1：高階Hilbert矩陣一般為病態(tài)矩陣，所以直接求逆可能出現錯誤結論，故應該采用invhilb(n)注意2：由于Hilbert矩陣本身接近奇異，所以建議處理該矩陣時建議盡量采用符號運算工具箱，若采用數值解時應該驗證結果的正確性第6頁，講稿共41頁，2023年5月2日，星期三托普利茲矩陣(toeplitz)toeplitz矩陣除第一行第一列外，其他每個元素都與左上角的元素相同。toeplitz(x,y)生成一個以x為第一列，y為第一行的托普利茲矩陣。這里x,y均為向量，兩者不必等長。toeplitz(x)用向量x生成一個對稱的托普利茲矩陣。例>>T=toeplitz([1:5],[1:7])T=12345672123456321234543212345432123第7頁，講稿共41頁，2023年5月2日，星期三帕斯卡矩陣二次項(x+y)n展開后的系數隨n的增大組成一個三角形表，稱為楊輝三角形。由楊輝三角形表組成的矩陣稱為帕斯卡(Pascal)矩陣。

pascal(n)生成一個n階帕斯卡矩陣。第8頁，講稿共41頁，2023年5月2日，星期三矩陣分解矩陣分解通過將復雜矩陣表示成形式簡單或具有良好數學性質（統(tǒng)稱為簡單矩陣）的組合，以便于理論分析或數值計算。通常矩陣分解將復雜矩陣分解為幾個簡單矩陣的乘積。求解線性方程組不可避免地要用到矩陣分解的概念。MATLAB中，線性方程組的求解主要用到三種基本的矩陣分解，即對稱正定矩陣的cholesky分解、一般方程的gaussian消去法和矩陣的正交分解。這三種分解由函數chol、lu和qr完成。第9頁，講稿共41頁，2023年5月2日，星期三矩陣的逆求方陣A的逆可用inv(A)說明1:對于Hilbert求逆時，不建議用inv,可用invhilb直接產生逆矩陣說明2：符號工具箱中也對符號矩陣定義了inv()函數,即使對更高階的奇異矩陣也可以精確的求解出逆矩陣說明3：對于奇異矩陣用數值方法的inv()函數，會給出錯誤的結果，但符號工具箱中的inv()會給出正確的結論第10頁，講稿共41頁，2023年5月2日，星期三例A=[16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1];det(A)B=inv(A)A=[16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1];A=sym(A)inv(A)第11頁，講稿共41頁，2023年5月2日，星期三矩陣的偽逆pinv(A)若A不是一個方陣，或A是一個非滿秩的方陣時，矩陣A沒有逆矩陣，但可以找到一個與A的轉置矩陣A’同型的矩陣B，使得：

A·B·A=A,B·A·B=B此時稱矩陣B為矩陣A的偽逆。例求矩陣A的偽逆

A=[0,0,0;0,1,0;0,0,1];pinv(A)第12頁，講稿共41頁，2023年5月2日，星期三矩陣分解矩陣分解通過將復雜矩陣表示成形式簡單或具有良好數學性質（統(tǒng)稱為簡單矩陣）的組合，以便于理論分析或數值計算。通常矩陣分解將復雜矩陣分解為幾個簡單矩陣的乘積。求解線性方程組不可避免地要用到矩陣分解的概念。MATLAB中，線性方程組的求解主要用到三種基本的矩陣分解，即對稱正定矩陣的cholesky分解、一般方程的gaussian消去法和矩陣的正交分解。這三種分解由函數chol、lu和qr完成。第13頁，講稿共41頁，2023年5月2日，星期三矩陣分解之LU分解矩陣的三角分解又稱LU分解，它的目的是將一個矩陣分解成一個下三角矩陣L和一個上三角矩陣U的乘積，即A=LU。Matlab使用函數lu實現LU分解，其格式為：

[L,U]=lu(A)其中U為上三角陣，L為下三角陣或其變換形式，滿足LU=A。

[L,U,P]=lu(A)U為上三角陣，L為下三角陣，P為單位矩陣的行變換矩陣，滿足LU=PA。例:A=[123;456;789];[L,U]=lu(A)[L,U,P]=lu(A)第14頁，講稿共41頁，2023年5月2日，星期三矩陣分解之Cholesky分解對于正定矩陣A，則存在一個實的非奇異上三角陣R，滿足R‘*R=A，稱為Cholesky分解。Matlab使用函數chol實現Cholesky分解，其格式為：

R=chol(A)若A非正定，則產生錯誤信息。

[R,p]=chol(A)不產生任何錯誤信息，若A為正定陣，則p=0，R與上相同；若A非正定，則p為正整數，R是有序的上三角陣。例：A=pascal(4);[R,p]=chol(A)第15頁，講稿共41頁，2023年5月2日，星期三矩陣分解之QR分解將矩陣A分解成一個正交矩陣Q與一個上三角矩陣R的乘積A=QR，稱為QR分解。Matlab使用函數qr實現QR分解，其格式為：

[Q,R]=qr(A) [Q,R,E]=qr(A)求得正交矩陣Q和上三角陣R，E為單位矩陣的變換形式，R的對角線元素按大小降序排列，滿足AE=QR。

[Q,R]=qr(A,0)例：A=[123;456;789;101112];[Q,R]=qr(A)第16頁，講稿共41頁，2023年5月2日，星期三代數方程的求解一般的代數方程包括：線性方程非線性方程超越方程超越方程一般沒有解析解，而只有數值解或近似解，只有特殊的超越方程才可以求出解析解來。matlab是獲得數值解的一個最強大的工具。常用的命令有fsolve,fzero等，但超越方程的解很難有精確的表達式，因此在matlab中常用eval（）函數得到近似數值解，再用vpa（）函數控制精度。第17頁，講稿共41頁，2023年5月2日，星期三求解線性方程→第18頁，講稿共41頁，2023年5月2日，星期三求解線性方程對于Ax=b，A為m*n矩陣。1、m=n（恰定方程組）2、m<n（不定方程）3、m>n（超定方程），多用于曲線擬合。第19頁，講稿共41頁，2023年5月2日，星期三線性方程的組的求解若系數矩陣的秩r=n(n為方程組中未知變量的個數)，則有唯一解；若系數矩陣的秩r<n，則可能有無窮解；求線性方程組的唯一解或特解（直接法）求線性方程組的通解線性方程組的無窮解=對應齊次方程組的通解+非齊次方程組的一個特解齊次線性方程組：常數項全部為零的線性方程組第20頁，講稿共41頁，2023年5月2日，星期三直接法求線性方程組求逆矩陣的思想:x=inv(A)*b或x=A^(-1)*b左除法（原理上是運用高斯消元法求解,但Matlab在實際執(zhí)行過程中是通過LU分解法進行的）：x=A\b符號矩陣法（此法最接近精確值，但運算速度慢）x=sym(A)\sym(b)化為行最簡形：C=[A,b]R=rref(C)則R的最后一列元素就是所求之解。如果無解，那么返回的是一個單位矩陣。第21頁，講稿共41頁，2023年5月2日，星期三例1：A=[-0.040.040.12;0.56-1.560.32;-0.241.24-0.28];b=[3;1;0];x11=inv(A)*b%法1：求逆思想x12=A^(-1)*b%法1：求逆思想x3=A\b%法2：左除法x4=sym(A)\sym(b)%法3：符號矩陣法C=[A,b];rref(C)%法4：化為行最簡形第22頁，講稿共41頁，2023年5月2日，星期三例:2：求下列方程組的解。A=[5600015600015600015600015];B=[10001]';R_A=rank(A) %求秩X=A\B %求解第23頁，講稿共41頁，2023年5月2日，星期三例3：求下列方程組的解。A=[11-3-1;3-1-34;15-9-8];B=[140]'; R_A=rank(A) %求秩X=A\B %由于系數矩陣不滿秩,該解法可能存在誤差X=[00-0.53330.6000]‘%一個特解近似值A=[11-3-1;3-1-34;15-9-8];B=[140]';C=[A,B];%構成增廣矩陣R=rref(C)若用rref求解，則比較精確：第24頁，講稿共41頁，2023年5月2日，星期三LU分解求方程組的解LU分解又稱Gauss消去分解，可把任意方陣分解為下三角矩陣和上三角矩陣的乘積。即

A=LU，L為下三角陣，U為上三角陣。則：

A*X=bL*U*X=b所以

X=U\(L\b)這樣可以大大提高運算速度。若采用[l，u，p]=lu(A)，則Ax=b的解為x=u(l\(p*b))第25頁，講稿共41頁，2023年5月2日，星期三利用LU分解求方程組的解例4:求下列方程組的一個特解。

解：A=[42-1;3-12;1130];B=[2108]';[L,U]=lu(A)X=U\(L\B)第26頁，講稿共41頁，2023年5月2日，星期三

例5：利用LU分解法求解方程組A=[2426;49615;26918;6151840]%先錄入系數矩陣b=[9232247]’[L,U]=lu(A)%對A矩陣做LU分解y=L\b%求解方程組Ly=b

x=U\y%求解方程組Ux=y得到方程組的最終解故方程組的最終解為:

x=(0.5000,2.0000,3.0000,-1.0000)’第27頁，講稿共41頁，2023年5月2日，星期三或A=[2426;49615;26918;6151840]%先錄入系數矩陣b=[9232247]’[L,U,P]=lu(A)%對A矩陣做LU分解y=L\P*b%求解方程組Ly=Pb

x=U\y%求解方程組Ux=y得到方程組的最終解第28頁，講稿共41頁，2023年5月2日，星期三Cholesky分解求方程組的解若A為對稱正定矩陣，則Cholesky分解可將矩陣A分解成上三角矩陣和其轉置的乘積，即：

A=R'*R其中R為上三角陣。方程A*X=b變成R‘*R*X=b所以 X=R\(R'\b)第29頁，講稿共41頁，2023年5月2日，星期三QR分解求方程組的解對于任何長方矩陣A，都可以進行QR分解，其中Q為正交矩陣，R為上三角矩陣的初等變換形式，即：

A=QR方程A*X=b變形成

QRX=b所以

X=R\(Q\b)若采用[Q，R，E]=qr(X)，則Ax=b的解為 x=E(R\(Q\b)) 第30頁，講稿共41頁，2023年5月2日，星期三解線性方程組的一般函數文件如下：function[x,y]=line_solution(A,b)[m,n]=size(A);y=[];

ifnorm(b)>0%b非零，非齊次方程組

ifrank(A)==rank([A,b])%方程組相容（有解）

ifrank(A)==n%有唯一解

x=A\b;

else%方程組有無窮多個解,基礎解系

disp('原方程組有有無窮個解，其齊次方程組的基礎解系為y，特解為x');y=null(A,'r');x=A\b;

end

第31頁，講稿共41頁，2023年5月2日，星期三

else%方程組不相容（無解），給出最小二乘解

disp('方程組的最小二乘法解是：');x=A\b;

endelse%齊次方程組

ifrank(A)==n%列滿秩

x=zero(n,1)%0解

else%非0解，無窮多個解

disp('方程組有無窮個解，基礎解系為x');x=null(A,'r');

end

endreturn第32頁，講稿共41頁，2023年5月2日，星期三如在MATLAB命令窗口，輸入命令

A=[2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2];b=[4,6,12,6]';[x,y]=line_solution(A,b)及：

A=[2,7,3,1;3,5,2,2;9,4,1,7];b=[6,4,2]';[x,y]=line_solution(A,b)分別顯示其求解結果。第33頁，講稿共41頁，2023年5月2日，星期三迭代法求解線性方程組迭代法適用于解大型稀疏方程組(萬階以上的方程組,系數矩陣中零元素占很大比例,而非零元按某種模式分布)背景:電路分析、邊值問題的數值解和數學物理方程常用迭代法：Jacobi迭代法Guass-Seidel迭代法SOR迭代法第34頁，講稿共41頁，2023年5月2日，星期三編寫實現Jacobo迭代法的函數jacobi.m如下:functions=jacobi(a,b,x0,err)%jacobi迭代法求解線性方程組，a為系數矩陣，b為%ax=b右邊的%矩陣b,x0為迭代初值，err為迭代誤差ifnargin==3 %nargin判斷輸入變量個數的函數

err=1.0e-6;elseifnargin<3errorreturnendD=diag(diag(a)); %構造對角陣DD1=inv(D); %求對角陣D的逆矩陣L=tril(a,-1); %構造嚴格下三角陣U=triu(a,1); %構造嚴格上三角陣B=-D1*(L+U); %求出迭代矩陣第35頁，講稿共41頁，2023年5月2日，星期三f=D1*b;s=B*x0+f;n=1;%n為迭代次數whilenorm(s-x0)>=errn=n+1x0=s;s=B*x0+f;s'endn第36頁，講稿共41頁，2023年5月2日，星期三A=[9-1-1;-110-1;-1-115];b=[7;8;13];x0=[0;0;0];err=0.00005;s=jacobi(A,b,x0,err)結果：n

1

0.77780.80000.866720.96300.96440.971930.99290.99350.995240.99870.99880.999150.99980.99980.9998