数值分析A-实验作业1 - 3.1-Toy模板网

这篇具有很好参考价值的文章主要介绍了数值分析A-实验作业1 - 3.1。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

数值分析A-实验作业1- 3.1

注：为了便于校验结果，本文为实验报告的补充

	可能用到的MATLAB函数
zeros(m,n)	生成m行，n列的零矩阵
ones(m,n)	生成m行，n列的元素全为1的矩阵
eye(n)	生成n阶单位矩阵
rand(m,n)	生成m行，n列(0,1)上均匀分布的随机矩阵
diag(x)	返回由向量x的元素构成的对角矩阵
tril(A)	提取矩阵A的下三角部分生成下三角矩阵
tril(A)	提取矩阵A的上三角部分生成上三角矩阵
rank(A)	返回矩阵A的秩
det(A)	返回方阵A的行列式
inv(A)	返回可逆方阵A的逆矩阵
[V,D]=eig(A)	返回方阵A的特征值和特征向量
norm(A,p)	矩阵或向量A的p范数
cond(A,p)	矩阵A的条件数
[L,U,P]=lu(A)	选取列主元LU分解
R=chol(X)	平方根分解
Hi=hilb(n)	生成n阶Hilbert矩阵

实验 3.1（主元的选取与算法的稳定性）

问题提出：

Gauss 消去法是我们在线性代数中已经熟悉的，但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保 Gauss 消去法作为数值算法的稳定性呢？ Grauss 消去法从理论算法到数值算法，其关键是主元的选择，主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

实验内容：

考虑线性方程组
$A\in R^{n \times n}，b \in R^n$
编制一个能自动选取主元，又能手动选取主元的求解线性代数方程组的 Gauss 消去过程.

实验要求1.

取矩阵 $A=\begin{bmatrix} 6& 1 & & & \\ 8 & 6& 1& & \\ & \ddots& \ddots &\ddots & \\ & & 8& 6& 1\\ & & & 8&6 \end{bmatrix}$ , $b=\begin{bmatrix} 7\\ 15\\ \vdots\\ 15\\ 14 \end{bmatrix}$ ，则方程有解 $x^*=(1,1,\cdots,1)^T$ ，取 $n = 10$ 计算矩阵的条件数.让程序自动选取主元，结果如何？

解：

（1）计算条件数

取 $n = 10$ ，计算矩阵A的条件数，

%-----  code 计算条件数 -----

% --- 输入矩阵A和b --- 
n = input(' 请输入n: ');
A = zeros(n,n);
b = zeros(n,1);

for i =1:n
    A(i,i) = 6;
end
for i=1:n-1
    A(i,i+1) = 1;
    A(i+1,i) = 8;
end

b(1,1) = 7;
b(2:n-1,1) = 15;
b(n,1) = 14;

% --- 计算矩阵A的3种条件数 --- 
disp('cond(A,1) ');
disp(cond(A,1));
disp('cond(A,2) ');
disp(cond(A,2));
disp('cond(A,inf) ');
disp(cond(A,inf));

得到：

$cond(A)_1=2557.5\\ cond(A)_2=1727.55602491356\\ cond(A)_\infty=2557.5$
$A$ 的条件数很大，说明矩阵 $A$ 为病态矩阵，方程组 $A x = b$ 为病态方程组，方程两端的微小扰动都可能对结果造成很大的误差。

（2）自动选取主元求解线性方程

输出说明：使用format long 输出（16位有效数字），计算精度评价指标：$\left |x-x^* \right |_{\infty} $

 顺序高斯消去法

%-----  code 顺序Gauss消去 -----
%-----  code 顺序Gauss消去 -----
%-----  code 顺序Gauss消去 -----



clc %清除命令窗口的内容
clear all; %清除工作空间的所有变量，函数，和MEX文件
format long ; %设置输出显示格式为16位有效数字

%---------------------  数值分析 实验3.1 主元的选取与算法的稳定性 ---------------------
    %A  系数矩阵
    %b  右端项
    %n  阶数
%------------------------  1. -------------------------


% --- --- --- --- --- ---输入矩阵A和b --- --- --- --- --- ---

n = input('Input n = ');
A = zeros(n,n);
b = zeros(n,1);

for i =1:n
    A(i,i) = 6;
end
for i=1:n-1
    A(i,i+1) = 1;
    A(i+1,i) = 8;
end

b(1,1) = 7;
b(2:n-1,1) = 15;
b(n,1) = 14;

% --- --- --- --- --- --- --- --- --- ------ --- --- --- ---


% --- --- --- --- --- --- 求解 --- --- --- --- --- --- --- 
disp('------- 精确解： ');
x = ones(n,1)

disp('------- 顺序高斯消元法解： ');
x0 = Gauss0(A,b)

disp('------- 误差 \left\| x-x^* \right\| _{\infty}  ： ');
norm(x-x0,inf)
% --- --- --- --- --- ---定义判断函数 --- --- --- --- --- ---
function r= iszero(a)
    if(abs(a-0)<1e-10)
        r=1;
        disp('被除数为0，方程组无解 或有无穷解')
    else
        r=0;
    end
end
% --- --- --- --- --- --- --- --- --- ------ --- --- --- ---


% --- --- --- --- --- ---定义顺序高斯消去函数 --- --- --- --- --- 

function x = Gauss0(A,b)
    %Gauss0为顺序高斯消去法
    %A  系数矩阵
    %b  右端项
    n = size(A,1);
    x = ones(n,1);
    l = ones(n,1);
    
    %消去过程
    for k = 1:n-1
        for i=k+1:n
            if(iszero(A(k,k)))%如果A(k,k)=0则停机
                return;
            else
                l(i) = A(i,k)/A(k,k);
                b(i) = b(i) - l(i)*b(k);
                for j =k:n
                    A(i,j) = A(i,j) - l(i)*A(k,j); 
                end
            end
        end
    end
    
    %回代过程
    if(iszero(A(n,n)))%如果(A(n,n)=0则停机
        return;
    else
       x(n) = b(n)/A(n,n);
    end

    for i=n-1:-1:1
        w=0;
        for j = i+1:n
            w = w + A(i,j)*x(j);
        end
        x(i) = (b(i)-w)/A(i,i);
    end

end

得到结果：

 列主元消去法

%-----  code 列主元Gauss消去 -----
%-----  code 列主元Gauss消去 -----
%-----  code 列主元Gauss消去 -----



clc %清除命令窗口的内容
clear all; %清除工作空间的所有变量，函数，和MEX文件
format long ; %设置输出显示格式为16位有效数字

%---------------------  数值分析 实验3.1 主元的选取与算法的稳定性 ---------------------
    %A  系数矩阵
    %b  右端项
    %n  阶数
%------------------------  1. -------------------------


% --- --- --- --- --- ---输入矩阵A和b --- --- --- --- --- ---

n = input('Input n = ');
A = zeros(n,n);
b = zeros(n,1);

for i =1:n
    A(i,i) = 6;
end
for i=1:n-1
    A(i,i+1) = 1;
    A(i+1,i) = 8;
end

b(1,1) = 7;
b(2:n-1,1) = 15;
b(n,1) = 14;

% --- --- --- --- --- --- --- --- --- ------ --- --- --- ---


% --- --- --- --- --- --- 求解 --- --- --- --- --- --- --- 
disp('------- 精确解为x ------- ');
x = ones(n,1)

disp('------- 列主元Gauss消元法解为x1 -------');
x1 = Gauss1(A,b)

disp('------- 误差为 \left\| x-x^* \right\| _{\infty} -------');
norm(x-x1,inf)
% --- --- --- --- --- ---定义判断函数 --- --- --- --- --- ---
function r= iszero(a)
    if(abs(a-0)<1e-10)
        r=1;
        disp('被除数为0，方程组无解 或有无穷解')
    else
        r=0;
    end
end
% --- --- --- --- --- --- --- --- --- ------ --- --- --- ---


% --- --- --- --- --- ---定义列主元Gauss消去函数 --- --- --- --- --- 
function x = Gauss1(A,b)  %Gauss1为列主元高斯消去法
    n = size(A,1);
    x = zeros(n,1);
    l = zeros(n,1);
        
    %消去过程
    for k = 1:n-1
        %找列主元
        MAX = max(abs(A(k:n,k)));%找到第k列，k到n行中绝对值最大的元素值
        index = find((abs(A(k:n,k)))==MAX);%得到第k到n行中最大元素的行标
        index = index+k-1;%转化为在矩阵A中的行标
        
        %列主元交换
        if abs(A(k,k)) ~= MAX %比index ~= k更好，因为可能有多个大小相同的最大值的情况
            index=index(1);
            A([k,index],:) = A([index,k],:);
            b([k,index],:) = b([index,k],:);
        end

        %消去过程
        for i=k+1:n
            if(iszero(A(k,k)))%如果A(k,k)=0则停机
                return;
            else
                l(i) = A(i,k)/A(k,k);
                b(i) = b(i) - l(i)*b(k);
                for j =k:n
                    A(i,j) = A(i,j) - l(i)*A(k,j); 
                end
            end
        end
    end

    %回代过程
    if(iszero(A(n,n)))%如果(A(n,n)=0则停机
        return;
    else
       x(n) = b(n)/A(n,n);
    end

    for i=n-1:-1:1
        w=0;
        for j = i+1:n
            w = w + A(i,j)*x(j);
        end
        x(i) = (b(i)-w)/A(i,i);
    end

end

得到结果：

 完全选主元消去法

%-----  code 完全选主元Gauss消去 -----
%-----  code 完全选主元Gauss消去 -----
%-----  code 完全选主元Gauss消去 -----




clc %清除命令窗口的内容
clear all; %清除工作空间的所有变量，函数，和MEX文件
format long ; %设置输出显示格式为16位有效数字

%---------------------  数值分析 实验3.1 主元的选取与算法的稳定性 ---------------------
    %A  系数矩阵
    %b  右端项
    %n  阶数
%------------------------  1. -------------------------


% --- --- --- --- --- ---输入矩阵A和b --- --- --- --- --- ---

n = input('Input n = ');
A = zeros(n,n);
b = zeros(n,1);

for i =1:n
    A(i,i) = 6;
end
for i=1:n-1
    A(i,i+1) = 1;
    A(i+1,i) = 8;
end

b(1,1) = 7;
b(2:n-1,1) = 15;
b(n,1) = 14;

% --- --- --- --- --- --- --- --- --- ------ --- --- --- ---


% --- --- --- --- --- --- 求解 --- --- --- --- --- --- --- 
disp('------- 精确解为x ------- ');
x = ones(n,1)

disp('------- 完全选主元Gauss消元法解为2 -------');
x2 = Gauss2(A,b)

disp('------- 误差为 \left\| x-x^* \right\| _{\infty} -------');
norm(x-x2,inf)
% --- --- --- --- --- ---定义判断函数 --- --- --- --- --- ---
function r= iszero(a)
    if(abs(a-0)<1e-10)
        r=1;
        disp('被除数为0，方程组无解 或有无穷解')
    else
        r=0;
    end
end
% --- --- --- --- --- --- --- --- --- ------ --- --- --- ---


% --- --- --- --- --- ---定义完全选主元Gauss消去函数 --- --- --- --- --- 
function x = Gauss2(A,b) %Gauss2为完全主元高斯消去法
    
    n = size(A,1);
    x = zeros(n,1);
    %l为乘数
    l = zeros(n,1);
    %e为初等变换矩阵，用于列交换之后将最后的解x的位置更正
    e = eye(n);
    %迭代n-1次
    for k = 1:n-1
        %寻找最大主元
        MAX1 = max(abs(A(k:n,k:n)));
        MAX = max(MAX1);
        [index1,index2] = find(abs(A(k:n,k:n))==MAX,1);
        index1 = index1+k-1;
        index2 = index2+k-1;
        %交换行、列
        if k ~= index1
            A([k,index1],:) = A([index1,k],:);
            b([k,index1],:) = b([index1,k],:);
        end
        if k ~= index2
            A(:,[k,index2]) = A(:,[index2,k]);
            e(:,[k,index2]) = e(:,[index2,k]);
        end
        for i=k+1:n
            l(i) = A(i,k)/A(k,k);
            b(i) = b(i) - l(i)*b(k);
            for j =k:n
                A(i,j) = A(i,j) - l(i)*A(k,j); 
            end
        end
    end

    x(n) = b(n)/A(n,n);
    for k=n-1:-1:1
        w=0;
        for j = k+1:n
            w = w + A(k,j)*x(j);
        end
        x(k) = (b(k)-w)/A(k,k);
    end
    x=e*x;

end

得到结果

实验要求2.

现选择程序中手动选取主元的功能，每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果.若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果.

（1）按模最小或按模尽可能小的元素作为主元

由完全选主元法就是按模最大的元素作为主元，若将取max换成取min，即按模最小或按模尽可能小的元素作为主元。因此对完全选主元函数myGauss2稍加修改，即可实现手动选取最小主元素。

同理，可以用程序模拟手动选取过程：迭代 $k$ 次时，在 $\begin{bmatrix} a_{kk}^{(k)} & \cdots &a_{kn}^{(k)}\\ \vdots & \ddots & \vdots \\ a_{nk}^{(k)} & \cdots &a_{nn}^{(k)} \end{bmatrix}$ 中按模最小选取主元。

%-----  code 手动按模最小或尽可能小Gauss消去 -----
%-----  code 手动按模最小或尽可能小Gauss消去 -----
%-----  code 手动按模最小或尽可能小Gauss消去 -----



clc %清除命令窗口的内容
clear all; %清除工作空间的所有变量，函数，和MEX文件
format long ; %设置输出显示格式为16位有效数字

%---------------------  数值分析 实验3.1 主元的选取与算法的稳定性 ---------------------
    %A  系数矩阵
    %b  右端项
    %n  阶数
%------------------------  实验要求（2） -------------------------


% --- --- --- --- --- ---输入矩阵A和b --- --- --- --- --- ---

n = input('Input n = ');
A = zeros(n,n);
b = zeros(n,1);

for i =1:n
    A(i,i) = 6;
end
for i=1:n-1
    A(i,i+1) = 1;
    A(i+1,i) = 8;
end

b(1,1) = 7;
b(2:n-1,1) = 15;
b(n,1) = 14;

% --- --- --- --- --- --- --- --- --- ------ --- --- --- ---


% --- --- --- --- --- --- 求解 --- --- --- --- --- --- --- 
disp('------- 精确解为x ------- ');
x = ones(n,1)

disp('------- 按模最小或尽可能小选取主元Gauss消元法解为x3 -------');
x3 = Gauss3(A,b)

disp('------- 误差为 \left\| x-x^* \right\| _{\infty} -------');
norm(x-x3,inf)
% --- --- --- --- --- ---定义判断函数 --- --- --- --- --- ---
function r= iszero(a)
    if(abs(a-0)<1e-10)
        r=1;
        disp('被除数为0，方程组无解 或有无穷解')
    else
        r=0;
    end
end
% --- --- --- --- --- --- --- --- --- ------ --- --- --- ---


% --- --- 定义 按模最小或尽可能小选取主元的Gauss消去 函数 --- --- --- 
 function x = Gauss3(A,b)
    format long;
    n = size(A,1);
    x = zeros(n,1);
    %l为乘数
    l = zeros(n,1);
    %e为初等变换矩阵，用于列交换之后将最后的解x的位置更正
    e = eye(n);
    %迭代n-1次
    for k = 1:n-1

        %按模最小或尽可能小选取主元
        MIN = Inf;
        for p=k:n
            for q=k:n
                if (abs(A(p,q))<MIN) && (abs(A(p,q))~=0)
                    MIN = abs(A(p,q));
                    index1 = p;
                    index2 = q;
                else
                    index1 = k;
                    index2 = k;
                end
            end
        end

        %交换行、列
        if k ~= index1
            A([k,index1],:) = A([index1,k],:);
            b([k,index1],:) = b([index1,k],:);
        end
        if k ~= index2
            A(:,[k,index2]) = A(:,[index2,k]);
            e(:,[k,index2]) = e(:,[index2,k]);
        end
        for i=k+1:n
            l(i) = A(i,k)/A(k,k);
            b(i) = b(i) - l(i)*b(k);
            for j =k:n
                A(i,j) = A(i,j) - l(i)*A(k,j); 
            end
        end
    end

    x(n) = b(n)/A(n,n);
    for k=n-1:-1:1
        w=0;
        for j = k+1:n
            w = w + A(k,j)*x(j);
        end
        x(k) = (b(k)-w)/A(k,k);
    end
    x=e*x;

end