文章目錄
- 線性方程組的解法
- 認識一些基本的矩陣函數
- MATLAB 實現
- 機電工程學院教學函數構造
- 1.高斯消元法
- 2.列主元消去法
- 3. L U LU LU分解法
線性方程組的解法
看到以下線性方程組的一般形式:設有以下的 n n n階線性方程組:
A x = b \mathbf{Ax}=\mathbf{b} Ax=b
求解線性方程組的方法可以分為兩類:直接法
和迭代法
。
- 直接法是指假設計算中不產生舍入誤差,結果有限次的運算可以得到方程組的精確解的方法,主要用于解低階稠密矩陣。
- 迭代法是一種通過構造迭代序列逐步逼近方程組精確解的方法。它將求解方程組的問題轉化為一個迭代格式,從一個初始近似解出發,按照一定的迭代公式反復計算,得到一系列近似解,當迭代次數足夠多時,這些近似解逐漸收斂到方程組的精確解。迭代法主要用于解高階稀疏矩陣方程組,因為對于高階稀疏矩陣,直接法可能會面臨計算量過大、存儲需求過高的問題,而迭代法可以利用矩陣的稀疏性,減少計算量和存儲空間。
認識一些基本的矩陣函數
函數 | 功能 |
---|---|
rank ( A ) \texttt{rank}(\mathbf{A}) rank(A) | 求 A \mathbf{A} A的秩,即 A \mathbf{A} A中線性無關的行數和列數 |
det ( A ) \texttt{det}(\mathbf{A}) det(A) | 求 A \mathbf{A} A的行列式 |
inv ( A ) \texttt{inv}(\mathbf{A}) inv(A) | 求 A \mathbf{A} A的逆矩陣,若 A \mathbf{A} A近似奇異,會拋出錯誤 |
pinv ( A ) \texttt{pinv}(\mathbf{A}) pinv(A) | 求 A \mathbf{A} A的偽逆 |
trace ( A ) \texttt{trace}(\mathbf{A}) trace(A) | 求 A \mathbf{A} A的跡,即對角線元素和 |
MATLAB 實現
在 M A T L A B MATLAB MATLAB中,使用運算符\
直接求解線性系統,該運算符功能強大,具有智能性。
x=A\b %求解線性系統 Ax=b
X=A\B %求解系統:AX=B
- 直接解法:
問題
{ x 1 + 3 x 2 ? 3 x 3 ? x 4 = 1 3 x 1 ? 6 x 2 ? 3 x 3 + 4 x 4 = 4 x 1 + 5 x 2 ? 9 x 3 ? 8 x 4 = 0 \begin{cases} x_1 + 3x_2 - 3x_3 - x_4&= 1 \\ 3x_1 - 6x_2 - 3x_3 + 4x_4 &= 4\\ x_1 + 5x_2 - 9x_3 - 8x_4 &= 0 \end{cases} ? ? ??x1?+3x2??3x3??x4?3x1??6x2??3x3?+4x4?x1?+5x2??9x3??8x4??=1=4=0?
A=[1,3,-3,-1;3,-6,-3,4;1,5,-9,-8];
B=[1,4,0]';
X=A\B
- 逆矩陣:
注意這種方法首先要看一下 A \mathbf{A} A是不是方陣。
x=A^-1*b
x=inv(A)*b
問題
{ x 1 + 2 x 2 = ? 1 3 x 1 + 4 x 2 = ? 1 \begin{cases} x_1 + 2x_2&= -1 \\ 3x_1+4x_2 &=-1 \end{cases} {x1?+2x2?3x1?+4x2??=?1=?1?
A=[1,2;3,4];b=[-1;-1];x=A^-1*b
3. L U LU LU分解
[L U]=lu(A)
X=U\(L\B)
問題
{ 4 x 1 + 2 x 2 ? x 3 = 2 3 x 1 ? x 2 + 2 x 3 = 10 11 x 1 + 3 x 2 + x 3 = 8 \left\{ \begin{array}{c} 4x_1+2x_2-x_3=2\\ 3x_1-x_2+2x_3=10\\ 11x_1+3x_2+x_3=8\\ \end{array} \right. ? ? ??4x1?+2x2??x3?=23x1??x2?+2x3?=1011x1?+3x2?+x3?=8?
>> A=[4 2 -1;3 -1 2; 11 3 1];
>> B=[2 10 8]';
>> D=det(A)D =-10.0000>> [L U]=lu(A)L =0.3636 -0.5000 1.00000.2727 1.0000 01.0000 0 0U =11.0000 3.0000 1.00000 -1.8182 1.72730 0 -0.5000>> X=U\(L\B)X =4.0000-10.0000-6.0000
機電工程學院教學函數構造
1.高斯消元法
代碼模板
function x = pureGaussianElimination(A, b)% 獲取矩陣 A 的行數n = size(A, 1);% 構建增廣矩陣 [A, b]augmentedMatrix = [A, b];% 前向消元過程for k = 1:n - 1for i = k + 1:n% 計算消元因子factor = augmentedMatrix(i, k) / augmentedMatrix(k, k);% 消去第 i 行第 k 列的元素augmentedMatrix(i, k:end) = augmentedMatrix(i, k:end) - factor * augmentedMatrix(k, k:end);endend% 回代求解x = zeros(n, 1);x(n) = augmentedMatrix(n, end) / augmentedMatrix(n, n);for i = n - 1:-1:1x(i) = (augmentedMatrix(i, end) - augmentedMatrix(i, i + 1:n) * x(i + 1:n)) / augmentedMatrix(i, i);end
end
2.列主元消去法
代碼模板
function x = gaussianElimination(A, b)% 獲取矩陣 A 的行數和列數[n, m] = size(A);% 構建增廣矩陣 [A, b]augmentedMatrix = [A, b];% 前向消元過程for k = 1:n-1% 選主元[~, pivotIndex] = max(abs(augmentedMatrix(k:n, k)));pivotIndex = pivotIndex + k - 1;% 交換行if pivotIndex ~= ktemp = augmentedMatrix(k, :);augmentedMatrix(k, :) = augmentedMatrix(pivotIndex, :);augmentedMatrix(pivotIndex, :) = temp;end% 消元for i = k+1:nfactor = augmentedMatrix(i, k) / augmentedMatrix(k, k);augmentedMatrix(i, k:end) = augmentedMatrix(i, k:end) - factor * augmentedMatrix(k, k:end);endend% 回代求解x = zeros(n, 1);x(n) = augmentedMatrix(n, end) / augmentedMatrix(n, n);for i = n-1:-1:1x(i) = (augmentedMatrix(i, end) - augmentedMatrix(i, i+1:n) * x(i+1:n)) / augmentedMatrix(i, i);end
end
3. L U LU LU分解法
代碼模板
function x = luDecomposition(A, b)% 獲取矩陣 A 的行數n = size(A, 1);% 初始化 L 為單位矩陣,U 為 AL = eye(n);U = A;% LU 分解過程for k = 1:n - 1for i = k + 1:n% 計算 L 矩陣的元素L(i, k) = U(i, k) / U(k, k);% 更新 U 矩陣U(i, k:end) = U(i, k:end) - L(i, k) * U(k, k:end);endend% 求解 Ly = by = zeros(n, 1);for i = 1:ny(i) = (b(i) - L(i, 1:i - 1) * y(1:i - 1)) / L(i, i);end% 求解 Ux = yfunction x = GaussianElimination(A, b)% 獲取矩陣 A 的行數n = size(A, 1);% 構建增廣矩陣 [A, b]augmentedMatrix = [A, b];% 前向消元過程for k = 1:n - 1for i = k + 1:n% 計算消元因子factor = augmentedMatrix(i, k) / augmentedMatrix(k, k);% 消去第 i 行第 k 列的元素augmentedMatrix(i, k:end) = augmentedMatrix(i, k:end) - factor * augmentedMatrix(k, k:end);endend% 回代求解x = zeros(n, 1);x(n) = augmentedMatrix(n, end) / augmentedMatrix(n, n);for i = n - 1:-1:1x(i) = (augmentedMatrix(i, end) - augmentedMatrix(i, i + 1:n) * x(i + 1:n)) / augmentedMatrix(i, i);end
endx = zeros(n, 1);for i = n:-1:1x(i) = (y(i) - U(i, i + 1:n) * x(i + 1:n)) / U(i, i);end
end