山东大学数值计算实验二（matlab）

网址：

数值计算实验二代码以及实验报告地址

实验题目：

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