数值计算讨论组

欢迎来到数值计算讨论组。 本论坛旨在增强师生交流的时效性,主要讨论数值计算相关问题,包括初等数值分析、矩阵计算、非线性方程(组)数值解法、微分方程数值解法等方面的内容。


    [旧]选主元消元法求矩阵的逆-程序示例

    分享

    Admin
    Admin

    帖子数 : 46
    注册日期 : 10-05-13

    [旧]选主元消元法求矩阵的逆-程序示例

    帖子 由 Admin 于 周四 五月 13, 2010 10:13 pm

    XuQiang

    今天上机课上,信计2班的王静同学给出了《数值分析引论》P300 习题8的一个Matlab程序,见附件liezhuyuanGJ.m。请同学们参考
    一下,并提一些意见。

    另外,我们提一个关于此程序的问题。程序中用选主元消元法求解矩阵A的逆矩阵时,是使用每一个对角线A(k,k)元素所在的行去消所有行中的A(i,k),其中1<= i <=n。请问这种消元法有什么优点和缺点?
    % liezhuyuanGJ.m
    % 王静 07级信计二班 0237
    % P300 习题8
    %A为非奇异阵时,用列主元高斯-约当法求A的逆
    %[A,B]用初等行变换将A变成E,对应的B变成A逆
    clear; clc;
    A = [2,1,-3,-1;3,1,0,7;-1,2,4,-2;1,0,-1,5];
    n = length(A);
    B = eye(n,n);
    C = [A,B];
    %d用于C的换行
    d = zeros(1,2*n);
    for k=1:n
    max = 0;
    maxrow = k;
    %求从A(k,k)到A(n,k)中绝对值最大的元素
    for i = k:n
    %dangqian=|A(i,k)|
    dangqian = A(i,k);
    if A(i,k) < 0
    dangqian = -1 * A(i,k);
    end
    %如果|A(i,k)|>max,将max改成该值,并记录该行
    if dangqian > max
    max = dangqian;
    maxrow = i;
    end
    end
    %从A(k,k)到A(n,k)元素都为0,A为奇异阵,退出
    if max == 0
    break;
    end
    %如果A(k,k)不是从A(k,k)到A(n,k)中绝对值最大的元素,将矩阵C的maxrow与k行互换
    if maxrow ~= k
    d = C(k,:);
    C(k,:) = C(maxrow,:);
    C(maxrow,:) = d;
    end
    %将C的C(k,k)应用行变换变为1
    ckk = C(k,k);
    for j = 1:2*n
    C(k,j) = C(k,j)/ckk;
    end
    %将C第k列除C(k,k)应用行变换全变为0
    for i = 1:n
    if i == k
    continue;
    end
    mik = C(i,k)/C(k,k);
    %将第i行减去mik倍的第k行
    for j = 1:2*n
    C(i,j) = C(i,j) - mik * C(k,j);
    end
    end
    end
    %invA为A的逆
    invA = C(1:n,5:8)

    ////////////////////////
    以下内容为回复内容 /////
    ////////////////////////



    friend726

    斗胆在这说一下自己的想法吧。
    首先,王静写的matlab程序,非常棒。巧妙的标记变量设置,程序层次清楚,注释恰当。非常欣赏。
    二。提一下建议。
    (1)最好将此程序设为函数,可使此程序的重用性提高。
    例如:function [inv] inverse(input);
    inv = zeros(size(input));
    inv = C(1:n,n+1:2*n);
    (2)下面语句中:我建议将“1”改为“k”。我认为会减少一些运算量。
    %将C的C(k,k)应用行变换变为1
    ckk = C(k,k);
    for j = 1:2*n
    C(k,j) = C(k,j)/ckk;
    end
    %将C第k列除C(k,k)应用行变换全变为0
    (3)程序中最后一句,我认为改为 C(1:n,n+1:2*n)比较好。
    其实我的方法,就是将原程序,转为函数,提高重用性。
    三。至于xuqiang老师的问题。我认为举个例子最好说明问题。
    例如:
    (0.000000000000000000333333333333333,2;10,4)
    10除于很小的数,误差是很大的。(或者很小的数做分母误差就会很大了。)
    优点:选主元方法,就会很好的避免A(k,k)很小的情况,使数值精确和稳定性更好。
    至于缺点,就是运算量更多一些。(多了n(n-1)/2步比较和至多n-1次行交换。我算的是这个。)

    其实我在此说的都是多余的,xuqiang老师的问题,数值老师有讲。而我的对于程序的改进,只要作者想变为函数,自然就会解决。
    所以,My name is nobody.



    Nobody

    补充一下,之所以高等代数没有选主元,而数值有选,主要是高代是理想的处理,而计算机确无法存储趋于无穷的小数或者无穷大的数。


    XuQiang
    选主元和变程序为函数的事情nobody说的很清楚了。

    关于(2)下面语句中:我建议将“1”改为“k”。我认为会减少一些运算量。
    %将C的C(k,k)应用行变换变为1
    ckk = C(k,k);
    for j = 1:2*n
    C(k,j) = C(k,j)/ckk;
    end
    %将C第k列除C(k,k)应用行变换全变为0

    我没有看明白是什么意思,这个消元求逆矩阵的过程是需要行变换为单位阵的,不知Nobody的理解是否有失偏颇。

    最后,U r not nobody, u will be somebody someday.


    friend726


    至于老师提出的问题,我想举个例子吧。
    例如:
    A = [1, 2, 3; 1, 4, 2 ;1, 3, 5;] (只是举个例子,是不是非奇异矩阵暂且不管吧,只是说一下过程)
    约化第一列会将A约化为 A = [1, 2, 3; 0, 2, -1; 0, 1, 3;]
    约化第二列时,要将第二行单位化,但是A(2,1)是 0 呀!,根本不需要单位化的。而程序中j是从1开始的。如果改为k,就会避免0的单位化,减少一些运算。我想虽然是0,但是Computer还是要运算一次的吧。
    谢谢老师的回复。有不当之处还望不吝指教。

    XuQiang

    为了节省计算量,是这样,没错。

      目前的日期/时间是周二 七月 17, 2018 7:48 am