数值计算实验报告

问题一:

  1. 问题描述
    问题1

  2. 算法分析

    • 高斯消去法
      高斯消去法的过程与我们平常进行人工解方程组的过程相似,即通过行倍加,倍乘变换,行交换等将矩阵化简为一个三角矩阵再进行求解。化简过程通过每次迭代计算出乘数因子进行化简。以第一次迭代为例:
      第一次迭代的乘数因子为

      mi1=a(1)i1/a(1)11 m i 1 = a i 1 ( 1 ) / a 11 ( 1 )

      再将 -m_{ij} 乘第一个方程,加到第i个方程。
      这就是我们常说的行倍加变换。
      需要注意的是,每次迭代的每行首位都不能为0,因此常常需要交换行

    • 列主元消去法
      列主元消去法和高斯消去法特别相似,只是增加了一些交换行的步骤。每次选取列主元中最大的一个,判断它是否为零,若非零,与第一行(相对于子矩阵)进行交换。

  3. 数值实验

    1. 高斯消去法
      编写函数guasselim进行高斯消元,并计算其运行时间。迭代过程通过循环完成

          for k =1:n-1;
              %compute the kth column of M
              m(k+1:n) = A(k+1:n,k)/A(k,k);
              %compute
              %An=Mn*An-1
              %bn=Mn*bn-1;
              for i=k+1:n
                  A(i,k+1:n) = A(i,k+1:n)-m(i)*A(k,k+1:n);
              end;
              b(k+1:n) = b(k+1:n)-b(k)*m(k+1:n);
          end
    2. 列主元消去法
      列主元消去法的特点在于对最大列主元的选择以及行交换,相当于对高斯消元法的一种改进。由于matlab提供了十分方便的寻找最大元素的函数,以及十分简便的交换语法,因此执行上述操作非常方便

            %find the max element of the kth column
            [max_a,index] = max(abs(A(1:n,k))); %取第k列绝对值最大数的值和索引
            if max_a == 0
                mdet = 0;
                return;
            end;
            A([k,index],:) = A([index,k],:);
            b([k,index],:) = b([index,k],:);  %交换k行和index行
            mdet = -mdet;
  4. 结果测试分析
    编写了测试函数,对满足正态分布的矩阵和向量A,b分别采用高斯消元法和列主元消去法进行求解,矩阵大小n = [10,50,100,200].
    结果如下:
    结果
    多次实验发现(这里不一一列出),列主元消去法的时间通常要比高斯消元法的时间要多,原因可能在于,列主元需要时时进行行交换,导致效率有点下降。不过,与高斯消元法相比较,列主元消去法的结果要比高斯消元法的结果更好。因为列主元消去法每次都会将最大的第一位元素所在的行移到第一行,于是计算乘数因子的时候,不会出现分母过小导致的误差变大的问题


问题二:

  1. 问题描述:
    2

  2. 算法分析
    说到迭代求解线性方程组,首先要讲的是一个基本方法。此后的雅各比迭代,高斯塞德尔迭代,超松弛迭代,都是基于这个方法进行操作的,
    假设有线性方程组

    Ax=b A x = b

    将A分裂为
    A=MN A = M − N

    则原线性方程组可以转化为求解
    Mx=Nx+b M x = N x + b


    x=M1Nx+M1b x = M − 1 N x + M − 1 b

    根据上式,我们就能够进行一阶定常迭代

    接下来,根据不同的算法来做具体分析

    1. 雅各比迭代法
      雅各比迭代法将A 化为 D - L - U, D是对角矩阵, L 和 U分别是下三角矩阵和上三角矩阵
      取M = D, N = L + U, 满足 A = M - N
      利用

      x=M1Nx+M1b x = M − 1 N x + M − 1 b
      进行迭代计算

    2. 高斯塞德尔迭代法
      令M = D - L, N = U, 同上进行迭代运算

    3. 超松驰迭代法

      M=1w(DwL) M = 1 w ( D − w L )
      , w称为松弛因子。代入迭代式中就能进行迭代。需要注意的是,松弛银子的范围设为(0,2)。当松弛因子值为1时,SOR方法即为高斯塞戴尔迭代法

    4. 共轭梯度法
      CG方法是一种十分常用的变分方法。所谓变分法,就是为了最终寻求极值函数,对应于求一个二次函数的极值。
      CG方法用于求解大型稀疏对称矩阵非常有效。

  3. 数值实验

    1. 雅各比迭代法

      for j=1:N
          x(j) = (b(j)-A(j,[1:j-1,j+1:N]) * P([1:j-1, j+1:N]))/A(j,j);
      end
    2. 高斯塞德尔迭代法

          for i=1:n
      
              sigma=0;
      
              for j=1:i-1
                      sigma=sigma+A(i,j)*x(j);
              end
      
              for j=i+1:n
                      sigma=sigma+A(i,j)*x_old(j);
              end
      
              x(i)=(1/A(i,i))*(b(i)-sigma);
          end
    3. 超松驰迭代法

      D = diag(diag(A));
      L = -tril(A,-1);
      U = -triu(A,1);
      B = (D - L * w)\((1 - w) * D + w * U);
      f = w * inv((D - L * w)) * b;
      x = B * x0 + f;
      n = 1;
    4. 共轭梯度法

          for k = 1:numel(b);
             r = r - t*z;
             err = norm(r);
             if( norm(r) < tol )
                  return;
             end
             B = (r'*z)/s;
             y = -r + B*y;
             z = A*y;
             s = y'*z;
             t = (r'*y)/s;
             x = x + t*y;
          end
  4. 结果测试分析
    这里写图片描述
    这里写图片描述
    由结果来看,随着n的增加,Jacobi 迭代法的运行效率持续下降,甚至从收敛变为发散。高斯塞戴尔次之,而共轭梯度法表现最好。对比上下两图,区别在于我将y轴范围扩大。可以看出,雅各比迭代法在迭代早期就开始发散。 而观察共轭梯度法,当n较小时,共轭梯度法只需要几步迭代就能接近最优解。而随着n增加,仍然保持着较好的效果。原因在于我们将A设为对称正定矩阵,而共轭梯度法十分适合这类矩阵

    观察以下时间运行比较
    这里写图片描述
    可见,共轭梯度法仍然是最优秀的算法

ps.如何生成一个特征值在0到1之间均匀分布的对称正定矩阵?

function [A] = generateSPDmatrix(n)
%A = rand(n,n); 
%A = 0.5*(A+A'); 
%A = A + n*eye(n); 
d = rand(1,n);
D = diag(d); %d作为对角元素
P = rand(n);  
[Q,~] = qr(P); %将P进行QR分解获得正交矩阵
A = Q*D*Q';  %对对角矩阵D进行相似变换获得A

end

问题三

  1. 问题描述
    这里写图片描述

    给定的网页显示如下信息
    这里写图片描述

  2. 算法分析

    1. 基本概念
      pagerank算法也是通过迭代来对多个节点进行排名。假设P为结果向量,向量元素为P的价值。那么有

      Pn+1=APn P n + 1 = A P n

      A=dS+1aNeeT A = d S + 1 − a N e e T

      • 矩阵S记录了不同节点之间的关系,从一个例子中来说明:
        这里写图片描述
        引用自pagerank从原理到实现
      • a表示权重,一般为0.85
      • (1-a) / N是为没有出度的节点设置的,假设他们对其他所有节点都有贡献(模拟了用户进入某个不链接到其他网页的网页上时,会通过输入网址随机访问其他网页)
    2. 实现过程
      显然,基本概念中的描述是对矩阵进行操作,可是本次实验要求给出的节点数为75879个,边数为508837。当你尝试着构造S矩阵时你会发现,matlab报错,原因提示为75879 * 75879 矩阵需要占据50G的内存,超出最大限制。这种方法行不通。后来与他人交流,发现能够使用稀疏矩阵进行操作,不过下面要讲的这种方法将更加简便。
      从上述的S矩阵例子可以看出,矩阵的每一列凡是非零的元素,值都相等。原因在于,sij表示的是j对于i的出链,这里假设为出链之和是1,出链的概率相同。因此,假如j有n个出链,那么每个出链的贡献即为1/n
      因此,我们可以得出我们的实现过程。

      1. 读入数据文件,得到一个75879 * 2 的矩阵
      2. 遍历数据矩阵的第一列,计算每种元素的个数并存入一个75879长度的数组。所获得数组即为每种元素的出度数量。根据出度数量就可以获得每种元素每条出链的贡献值
      3. 遍历数据矩阵的第二列,对于不同的数字j,我们可以找到对应的第一列数字i,这个数字即为i对j的入度。根据i去上一步的数组中寻找对应的贡献值,与旧结果数组的对应位置相乘后加到新的结果数组的对应元素中。
      4. 最后一步,即为结果数组中的每一个值添加一个 (1-a)/ N * x 。x为对应旧结果数组的每个值。

      可以发现,我们的实现方法完全模拟矩阵与向量相乘,但是这是在不构造矩阵的情况下完成。因为发现(1-a)/ N 对于任何元素(不管有没有出度),因此留到了最后直接与所有旧数组元素相乘之后加到新数组的每个元素。这个实现过程很简单,调用sum函数计算旧数组中所有元素之和,再乘以(1-a) / N

    3. 实现过程中的问题。

      1. 寻找每个元素入度过程中出现的问题
        这是未修改前的代码:

        %    for i = 1:MAXINDEX  %对结果数组的每一个进行操作
        %        index = i - 1;
        %        tempResult = 0;
                for j = 1:EDGE
                    index = soc_Epinions1(j,2);
                    target = soc_Epinions1(j,1);  %target是从0开始
                    tempA = a * (1/ODnum(target + 1));
                    tempJ = tempA * oldResult(target + 1);
                    Result(index+1) = Result(index+1) + tempJ;
                end
        %        tempResult = tempResult + ((1-a)/N) * sum(oldResult);
        %        Result(i) = tempResult;
                for j = 1:MAXINDEX
                    if Result(j) ~= 0
                        Result(j) = Result(j) + ((1-a)/N) * sum(oldResult);
                    end
                end
        
                change = change + abs(sum(Result - oldResult));
        %    end

        从代码可以看出,为了寻求不同元素的入度,我在外层(变量i所在层)进行循环,每次选择特定的元素i,然后遍历数据矩阵去寻找i的入度。数据矩阵有50万条数据,节点有7万个,因此一次迭代下来,我需要做350万次计算
        毫无疑问,结果是跑不出来的
        因此我将代码做出以下修改:

                for j = 1:EDGE
                    index = soc_Epinions1(j,2);
                    target = soc_Epinions1(j,1);  %target是从0开始
                    tempA = a * (1/ODnum(target + 1));
                    tempJ = tempA * oldResult(target + 1);
                    Result(index+1) = Result(index+1) + tempJ;
                end
                for j = 1:MAXINDEX
                    if Result(j) ~= 0
                        Result(j) = Result(j) + ((1-a)/N) * sum(oldResult);
                    end
                end

        思想就在于,不再选择特定的元素,而是在扫描数据矩阵的第二列元素的时候,扫到什么元素就算什么元素,将计算结果加到Result对应元素里面。因此,一次扫描就能完成所有元素的值的更新。最后,再将 (1-a) / N 加到Result的每个元素里。这里有需要一次循环。因此总共迭代大概为57万次。

      2. 数据矩阵的下标是从0开始的,而matlab数组的下标从1开始。
        这个问题,只需要用进行1的增加或减少即可。在代码中,我用index表示数据下标(从0开始),用i表示数组下标(从1开始)

      3. 数据节点数为75879个,而查询数据文件发现,最大下标为75887。这是,我开始担心,节点标号并不是连续的。毫无疑问,我必须申请一个长度大于75888的数组来存储,但因此也必然会出现,数组中的某些元素,在实际数据文件中并不存在,因此迭代过程中必须让这些元素保持为0(就如同不存在这个网页,就不能对他进行操作)。因此,在初始化和迭代过程中,我小心的避免了这些节点

        for i = 1:EDGE  %初始化
            index1 = soc_Epinions1(i,1); %index 从零开始
            index2 = soc_Epinions1(i,2);
            Result(index1 + 1) = 1/N;
            Result(index2 + 1) = 1/N;
            ODnum(index1 + 1) = ODnum(index1 + 1) + 1;
        end
        
                for j = 1:MAXINDEX  %迭代
                    if Result(j) ~= 0
                        Result(j) = Result(j) + ((1-a)/N) * sum(oldResult);
                    end
                end
    4. 结果分析
      这里写图片描述
      由下至上,排名逐渐降低,第一列为元素标号,第二列为对应价值。从图中看出,第一名为18号元素
      精度设置为 1e-4
      迭代次数为1次

课题一: 线性方程组的迭代法 一、实验内容 1、设线性方程组 = x = ( 1, -1, 0, 1, 2, 0, 3, 1, -1, 2 ) 2、设对称正定阵系数阵线方程组 = x = ( 1, -1, 0, 2, 1, -1, 0, 2 ) 3、三对角形线性方程组 = x = ( 2, 1, -3, 0, 1, -2, 3, 0, 1, -1 ) 试分别选用Jacobi 迭代法,Gauss-Seidol迭代法和SOR方法计算其解。 二、实验要求 1、体会迭代法求解线性方程组,并能与消去法做以比较; 2、分别对不同精度要求,如 由迭代次数体会该迭代法的收敛快慢; 3、对方程组2,3使用SOR方法时,选取松弛因子 =0.8,0.9,1,1.1,1.2等,试看对算法收敛性的影响,并能找出你所选用的松弛因子的最佳者; 4、给出各种算法的设计程序和计算结果。 三、目的和意义 1、通过上机计算体会迭代法求解线性方程组的特点,并能和消去法比较; 2、运用所学的迭代法算法,解决各类线性方程组,编出算法程序; 3、体会上机计算时,终止步骤 (予给的迭代次数),对迭代法敛散性的意义; 4、体会初始解 x ,松弛因子的选取,对计算结果的影响。 课题二:数值积分 一、实验内容 选用复合梯形公式,复合Simpson公式,Romberg算法,计算 (1) I = (2) I = (3) I = (4) I = 二、实验要求 1、 编制数值积分算法的程序; 2、 分别用两种算法计算同一个积分,并比较其结果; 3、 分别取不同步长 ,试比较计算结果(如n = 10, 20等); 4、 给定精度要求 ,试用变步长算法,确定最佳步长。 三、目的和意义 1、 深刻认识数值积分法的意义; 2、 明确数值积分精度与步长的关系; 3、 根据定积分的计算方法,可以考虑二重积分的计算问题。 四、流程图设计
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值