山東大學數值計算實驗二(matlab)
阿新 • • 發佈:2019-01-25
實驗題目:
1、高斯消元法 Computer Problems P100 2.2,(a)(b)
(1)用MATLAB語言提供的方法解方程組(x=A\b)
(2)寫出直接LU分解函式(參考例2.13),給出A的直接LU分解結果
(3)寫前代、回代函式。結合LU分解,前代、回代解(a)(b)方程組。
(4)直接使用軟體環境提供的函式LU分解函式、解方程組方法比較
說明:(2)、(3),這幾部分不得直接使用軟體環境提供的函式完成,可以參考課本上的演算法流程實現
實驗程式碼:
%實驗二第一問 解方程組(x=A\b)
A = [2,4,-2;4,9,-3;-2,-1,7];
b = [2;8;10];
c = [4;8;-6];
x = A\b;
y = A\c;
fprintf('結果為:');
x
y
mylu函式:
function [ L,U,p ] = mylu( A )
%UNTITLED2 此處顯示有關此函式的摘要
% 此處顯示詳細說明
% LU 三角分解
%[L,U,p] = mylu(A)
%生成單位下三角矩陣L,上三角矩陣U,以及排列向量p
%滿足:L*U = A(p,:)
%本函式可以解決主元位置上有0的情況
%獲得A的行數和列數
[n,n] = size(A);
%定義排列向量p
p = (1:n)';
%for迴圈,用 k 記錄消元步數
for k = 1:n-1
% 為了排除第k列全為0以及主元位置為0
%先找出第k列的對角元及以下元素的絕對值的最大值
[r,m] = max(abs(A(k:n,k))); % r為最大值,m 為相對行數
m = m+k-1; %此時 m 為最大值在A中對應的 行數
%如果第k列全為0 ,則跳過消元過程
if (A([m,k]) ~= 0)
%避免主元為0,將絕對值最大值交換到主元所在行,現在的主元是對角元
if( m ~= k)
A([k m],:) = A([m k],:); %矩陣A兩行進行交換
p([k m]) = p([m k]) ; %排列向量p 兩行進行交換
end
% 計算第k列元素除以主元后得到的數
i = k+1:n;
A(i,k) = A(i,k)/A(k,k);
%更新矩陣剩餘部分
j = k+1:n;
A(i,j) = A(i,j) - A(i,k)*A(k,j);
end
end
%分離結果
L = tril(A,-1) + eye(n,n); %tril()函式提取矩陣的下三角部分,eye()獲得單位矩陣
U = triu(A); %tril()函式提取矩陣的上三角部分
end
利用mylu函式程式碼如下:
%實驗二的(2)寫出直接LU分解函式,給出A的直接LU分解結果
A = [2,4,-2;4,9,-3;-2,-1,7];
[L,U,p] = my_lu(A);
fprintf('矩陣A的 L*U 分解為:');
L
U
fprintf(' L*U 形成矩陣 LU 為:');
LU = L*U
fprintf(' 排列向量 p 為:');
p'
fprintf('矩陣 A 為:');
A
(3)兩個函式,前代函式以及後代函式:
function x = myforward(L,x)
%my_forward 前向消元,前代
%對於下三角矩陣L, x = my_forward(L,b)給出 L*x = b 的解x
[n,n] = size(L);
%消元過程
for k = 1:n
j = 1:k-1 ;
x(k)= x(k) - L(k,j)*x(j);
x(k)= x(k)/L(k,k);
end
end
function x = myback(U,x)
% my_back 後向消元
%對於上三角矩陣U, x = my_back(U,y) 給出 U*x = y 的解x
[n,n] = size(U);
%消元過程
for k = n:-1:1
j = k+1:n ;
x(k) = (x(k) - U(k,j)*x(j))/U(k,k);
end
end
使用自編函式求解程式碼:
%第三問 寫前代、回代函式。結合LU分解,前代、回代解(a)(b)方程組。
A = [2,4,-2;4,9,-3;-2,-1,7];
b = [2;8;10];
c = [4;8;-6];
[L,U,p] = mylu(A);
%重新排列 b ,向前消元,L*U*x1 = b,其中 L*y1 = b
y1 = myforward(L,b(p));
%反向回代
x1 = myback(U,y1);
fprintf('利用前代、回代函式,結合LU分解,得到的答案 x 為:')
x = x1
%重新排列 c ,向前消元,L*U*x2 = c,其中 L*y2 = c
y2 = myforward(L,c(p));
%反向回代
x2 = myback(U,y2);
fprintf('利用前代、回代函式,結合LU分解,得到的答案 y 為:')
y = x2
%實驗二的第四問 直接使用軟體環境提供的函式LU分解函式
[L,U,P] = lu(A);
fprintf('得到的 L*U 分解為:');
L
U