实验题目:
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