精华内容
下载资源
问答
  • MATLAB数学实验——Jacobi迭代法&Gauss-Seidel迭代法 一、迭代算法的数学知识 二、Jacobi迭代法的MATLAB实现 三、Gauss迭代法的MATLAB实现 四、例子

    MATLAB数学实验——Jacobi迭代法&Gauss-Seidel迭代法

    一、迭代算法的数学知识

    线性方程组的数值求解方法,有经典的Jacobi和Gauss-Seidel迭代方法。
    二者通过迭代,从而求解方程。
    基本思路:
    ①将矩阵方程AX=b中的A分解为(U+L+D),即上三角矩阵、下三角矩阵和对角矩阵之和;
    ②将等式化为:Xk+1=BXk+d的格式,从而求得X。

    (1) Jacobi迭代法

    A=U+L+D
    AX=b
    →(U+L+D)X=b
    →DX=-(U+L)X+b
    →DX=-(A-D)X+b
    →X=(E+D-1A)X+D-1b=BJX+dJ

    Xk+1=BJXk+dJ
    BJ=E+D-1A;
    dJ=D-1b;

    (2) Gauss-Seidel迭代法

    A=U+L+D
    AX=b
    →(U+L+D)X=b
    →(L+D)X=-UX+b
    →X=-(L+D)-1UX+(L+D)-1b

    Xk+1=BGXk+dG
    BG=-(L+D)-1U;
    dG=(L+D)-1b;

    二、Jacobi迭代法的MATLAB实现

    调用格式:
    X=Jacobi_2(A,b,X0,Norm,Error,Max,p)
    Norm:范数的名称,Norm=1,2,inf;
    error:近似解的误差;
    Max:迭代的最大次数;
    p:是否需要画图(不输入则不画),p=1,0,不输入

    %用Jacobi迭代法解线性方程组Ax=b
    %Norm:范数的名称,Norm=1,2,inf;
    %error:近似解的误差;
    %Max:迭代的最大次数;
    function X=Jacobi_2(A,b,X0,Norm,Error,Max,p)
    if nargin==6
        p=0;
    end
    a=[];
    x=[];
    [N N]=size(A);
    X=X0;
    [L,D,U]=LUD(A);
    B=eye(N)-inv(D)*A;
    d=inv(D)*b;
    X1=A\b;
    Result=lt_con(B);
    if Result~=1
        error('迭代算法不收敛');
        return
    end
    disp '迭代算法收敛,继续计算...'
    for i=1:Max
        X=B*X+d;
        errX=norm(X-X1,Norm);
        %X0=X;
        a(i)=errX;
        x=i;
        if errX<Error
            diedaicishu=i;
            if p==1
                disp('迭代次数:')
                disp(diedaicishu)
                plot(1:x,a);
            end
            return
        end
    end
    if errX>=Error
        disp('请注意:Jacobi迭代次数已经超过最大迭代次数Max.')
    end
    if p==1
        plot(1:x,a);
    end
    
    end
    

    三、Gauss迭代法的MATLAB实现

    调用格式:
    X=G_S(A,b,X0,Norm,Error,Max,p)
    Norm:范数的名称,Norm=1,2,inf;
    error:近似解的误差;
    Max:迭代的最大次数;
    p:是否需要画图(不输入则不画),p=1,0,不输入

    %用Gauss-Seidel迭代法解线性方程组Ax=b
    %Norm:范数的名称,Norm=1,2,inf;
    %error:近似解的误差;
    %Max:迭代的最大次数;
    function X=G_S(A,b,X0,Norm,Error,Max,p)
    if nargin==6
        p=0;
    end
    a=[];
    x=[];
    [N N]=size(A);
    X=X0;
    [L,D,U]=LUD(A);
    B=-inv(D+L)*U;
    d=inv(D+L)*b;
    X1=A\b;
    Result=lt_con(B);
    if Result~=1
        error('迭代算法不收敛');
        return
    end
    disp '迭代算法收敛,继续计算...'
    % disp '...'
    % pause(0.1)
    % disp '..'
    % pause(0.1)
    % disp '.'
    % pause(1)
    for i=1:Max
        X=B*X+d;
        errX=norm(X-X1,Norm);
        %X0=X;
        a(i)=errX;
        x=i;
        if errX<Error
            diedaicishu=i;
            if p==1
                disp('迭代次数:')
                disp(diedaicishu)
                plot(1:x,a);
            end
            return
        end
    end
    if errX>=Error
        disp('请注意:Gauss-Seidel迭代次数已经超过最大迭代次数Max.')
    end
    if p==1
        plot(1:x,a);
    end
    end
    

    四、例子

    A=[2,-1,1; 3,6,2;3,3,7];
    b=[15,5,8]';
    X0=[0,0,0]';
    Jacobi_2(A,b,X0,inf,1e-3,1000,1);
    hold on
    G_S(A,b,X0,inf,1e-3,1000,1);
    

    例子结果:
    (1)命令行窗口:
    在这里插入图片描述

    (2)绘图窗口(蓝色为Jacobi迭代误差,红色为Gauss-Seidel迭代误差,横轴为迭代次数):
    在这里插入图片描述

    [20200212]更新:补充lt_con()子程序代码:

    %计算迭代的收敛性,作为迭代的子程序
    %
    function Result=lt_con(B)
    syms k;
    l=length(B);
    C=zeros(size(B));
    for i=1:l
        C(i)=limit(B(i)^k,k,inf);
    end
    Crit=C;
    %Crit=limit(B^k,k,inf);
    if Crit==0
        Result=1;
    else
        Result=0;
    end
    end
    

    [20210104]更新:补充LUD()子程序代码:

    function [L U D]=LUD(A)
    [n m]=size(A);
    L=zeros(size(A));
    U=zeros(size(A));
    D=zeros(size(A));
    if n<=0
        error('输入矩阵错误')
        return;
    end
    if n~=m
        error('输入矩阵错误')
        return;
    end
    for i=1:n-1
        L(i+1:end,i)=A(i+1:end,i);
        U(i,i)=A(i,i);
        D(i,i+1:end)=A(i,i+1:end);
    end
    U(n,n)=A(n,n);
    end
    

    确保有以上子程序,迭代法才能顺利运行。

    转载请注明出处

    展开全文
  • 为了找到任何非线性方程的根,仅在第一次迭代中获得的值不足以找到完美的根,必须执行迭代次数以获得清晰准确的根
  • 数学笔记9——牛顿迭代法

    万次阅读 多人点赞 2017-09-25 18:22:21
    牛顿迭代法(Newton's method)又称为牛顿-拉夫逊(拉弗森)方法(Newton-Raphson method),它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。 示例1:求解平方根  先来看如何用牛顿迭代法求解...

      牛顿迭代法(Newton's method)又称为牛顿-拉夫逊(拉弗森)方法(Newton-Raphson method),它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。

    示例1:求解平方根

      先来看如何用牛顿迭代法求解5的平方根。在计算器上的结果是2.236067…

      问题可以看作解方程x2=5,下面尝试用牛顿迭代法求解。

      首先令f(x)= x2 – 5 = 0,这是标准步骤,取得一个新函数,令该函数为0。这是一个抛物线:

      抛物线与x轴的交点x就是方程的解,它比2稍大一点。

      现在在x=2处对f(x)做切线:

    f(x)的切线

    切线与x轴的交点

      x0=2,y0= x02 – 5 = -1,设k是切线斜率:

      在x1处做f(x)的切线,重复上面步骤,

      这就是牛顿迭代法的公式。通过作图可以看出,每一次迭代,x都将更靠近最终解。

      f’(x)=2x,将公式代入目标方程f(x)=x2-5:

      已经相当接近计算器的结果。

    示例2:2cosx=3x

      解方程2cosx=3x

      由图像可知,方程存在唯一解。

      f(x)=2cosx-3x=0,f’(x)=-2sinx-3,x0=π/6≈0.52

    注意事项

      牛顿迭代法几乎可以求解所有方程,但它仍然有一些限制。

      通过前两个例子可以看到,在使用牛顿迭代法时,需要选取一个较为解接近真实解的x0作为迭代基数,x0如何选取呢?一句参考是:“f’不能太小,f’’不能太大,x0要在x附近”,这似乎要凭经验和感觉了,没有什么太好的办法;实际上,如果x0和x的差距过大,可能会得到一个没谱的解。

      设第n次迭代的误差是En=|x-xn|,那么需要满足En+1<En。如果选择和计算都正确,误差缩小的速度将非常快。

      以计算5的平方根为例,如果选择x0=-2,结果将偏向于-2.236067…;如果选择x0=0,则f’(0)=0,没法继续迭代,函数曲线如下图所示:

    选择了错误的x0

    代码示例:牛顿迭代法开平方

      设x2=a,则f(x)= x2-a,根据牛顿迭代法公式:

     1 const float EPS = 0.00001; 
     2 double sqrt(double x) { 
     3     if(x == 0) 
     4         return 0; 
     5     double result = x; 
     6     double lastValue; 
     7     do{ 
     8         lastValue = result; 
     9         result = result / 2.0f + x / 2.0f / result; 
    10     }while(abs(result - lastValue) > EPS);
    11     return (double)result;
    12  } 

       上面方法开平方会很快,但 https://www.2cto.com/kf/201206/137256.html 中提到了一个更快的方法。

      1999年12月,美国id Software公司发布了名为“雷神之锤III”的电子游戏。它是第一个支持软件加速的游戏,取得了极大成功。(由于影响力过大,文化部于2004年将它列入了非法游戏名单)雷神之锤III并不是id Software公司的第一次成功。早在1993年开始,这家公司就以“毁灭战士”系列游戏名闻天下。1995年,“毁灭战士”的安装数超过了当年微软的windows 95。据传比尔盖茨才曾经考虑买下id software。(id software公司后来被推出过“上古卷轴”系列的Bethesda公司买下)

      id Software所取得的成功很大程度上要归功于它的创始人约翰·卡马克。马克尔也是一个著名的程序员,他是id Software游戏引擎的主要负责人。 回到刚才提到的雷神之锤,马克尔是开源软件的积极推动者,他于2005年公布了雷神之锤III的源代码。至此人们得以通过研究这款游戏引擎的源文件来查看它成功的秘密。

      在其中一个名字为q_math.c的文件中发现了如下代码段:

     1 float Q_rsqrt( float number ) { 
     2     long i; float x2, y; const float threehalfs = 1.5F;
     3     x2 = number * 0.5F; 
     4     y = number; 
     5     i = * ( long * ) &y; // evil floating point bit level hacking 
     6     i = 0x5f3759df - ( i >> 1 ); // what the fuck? 
     7     y = * ( float * ) &i; 
     8     y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration 
     9     // y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
    10     #ifndef Q3_VM #
    11     ifdef __linux__ assert( !isnan(y) ); // bk010122 - FPE?
    12     #endif
    13     #endif return y; 
    14 }

      这段代码的作用就是求number的平方根,并且返回它的倒数。

      经过测试,它的效率比上述牛顿法程序要快几十倍。也比c++标准库的sqrt()函数要快好几倍。此段代码有一个奇怪的句子:

      i = 0x5f3759df - ( i >> 1 ); // what the fuck? 

      没错,一般的求平方根都是这么循环迭代算的但是卡马克(quake3作者)真正牛B的地方是他选择了一个神秘的常数0x5f3759df 来计算那个猜测值,就是我们加注释的那一行,那一行算出的值非常接近1/sqrt(n),这样我们只需要2次牛顿迭代就可以达到我们所需要的精度。好吧如果这个还不算NB,接着看:

      普渡大学的数学家Chris Lomont看了以后觉得有趣,决定要研究一下卡马克弄出来的这个猜测值有什么奥秘。Lomont也是个牛人,在精心研究之后从理论上也推导出一个最佳猜测值,和卡马克的数字非常接近, 0x5f37642f。卡马克真牛,他是外星人吗?

      传奇并没有在这里结束。Lomont计算出结果以后非常满意,于是拿自己计算出的起始值和卡马克的神秘数字做比赛,看看谁的数字能够更快更精确的求得平方根。结果是卡马克赢了... 谁也不知道卡马克是怎么找到这个数字的。

      最后Lomont怒了,采用暴力方法一个数字一个数字试过来,终于找到一个比卡马克数字要好上那么一丁点的数字,虽然实际上这两个数字所产生的结果非常近似,这个暴力得出的数字是0x5f375a86。

      Lomont为此写下一篇论文,"Fast Inverse Square Root"。

      论文下载地址:
      http://www.math.purdue.edu/~clomont/Math/Papers/2003/InvSqrt.pdf
      http://www.matrix67.com/data/InvSqrt.pdf

    向牛顿致敬

      牛顿是近代科学的先驱,智商290,在多个领域都有非凡的成就。

      他在1687年发表的论文《自然定律》里,对万有引力和三大运动定律进行了描述。这些描述奠定了此后三个世纪里物理世界的科学观点,并成为了现代工程学的基础。他通过论证开普勒行星运动定律与他的引力理论间的一致性,展示了地面物体与天体的运动都遵循着相同的自然定律;为太阳中心说提供了强有力的理论支持,并推动了科学革命。

      在力学上,牛顿阐明了动量和角动量守恒的原理,提出牛顿运动定律[1]  。在光学上,他发明了反射望远镜,并基于对三棱镜将白光发散成可见光谱的观察,发展出了颜色理论。他还系统地表述了冷却定律,并研究了音速。

      在数学上,牛顿与戈特弗里德·威廉·莱布尼茨分享了发展出微积分学的荣誉。他也证明了广义二项式定理,提出了“牛顿法”以趋近函数的零点,并为幂级数的研究做出了贡献。

      在经济学上,牛顿提出金本位制度。

      在天文成上,牛顿1672年创制了反射望远镜。他还用万有引力原理说明潮汐的各种现象,指出潮汐的大小不但同月球的位相有关,而且同太阳的方位有关。牛顿预言地球不是正球体。

      在哲学成上,牛顿的哲学思想基本属于自发的唯物主义,他承认时间、空间的客观存在。如同历史上一切伟大人物一样,牛顿虽然对人类作出了巨大的贡献,但他也不能不受时代的限制。例如,他把时间、空间看作是同运动着的物质相脱离的东西,提出了所谓绝对时间和绝对空间的概念;他对那些暂时无法解释的自然现象归结为上帝的安排,提出一切行星都是在某种外来的“第一推动力”作用下才开始运动的说法。《自然哲学的数学原理》牛顿最重要的著作,1687年出版。

      向伟大的牛顿致敬!

    总结

    1. 牛顿迭代法的公式:
    2. 注意事项,f’不能太小,f’’不能太大,x0要在x附近

        作者:我是8位的

      出处:http://www.cnblogs.com/bigmonkey

      本文以学习、研究和分享为主,如需转载,请联系本人,标明作者和出处,非商业用途! 

     

    展开全文
  • 第三节课程,介绍的是迭代法。 前两节笔记的文章: 程序员的数学笔记1–进制转换 程序员的数学笔记2–余数 03 迭代法 什么是迭代法 迭代法,简单来说,其实就是不断地用旧的变量值,递推计算新的变量值。 这里...

    第三节课程,介绍的是迭代法。

    前两节笔记的文章:


    03 迭代法

    什么是迭代法

    迭代法,简单来说,其实就是不断地用旧的变量值,递推计算新的变量值。

    这里采用一个故事来介绍什么是迭代法,这个故事是讲述一个国王要重赏一个做出巨大贡献的臣子,让臣子提出他想得到的赏赐,这个聪明的臣子说出了他想得到的赏赐–在棋盘上放满麦子,但要求是每个格子的麦子数量都是前一个格子的两倍。国王本以为这个赏赐可以轻而易举的满足,但真正开始放麦子后,发现即便是拿出全国的粮食也无法满足的臣子的这个赏赐。

    这里我们可以用f(n)表示当前各自的麦子数量,而前一个格子的麦子数量就是f(n-1),那么臣子的要求就可以这么表示:

    f(n) = f(n-1) * 2
    f(1) = 1
    

    这也就是迭代法了,而如果用编程来实现,其实就是实现一个循环运算的过程。

    用 Python 实现这个计算麦子的代码如下所示:

    def get_number_of_wheat(grid):
        '''
        \计算放到给定格子数量需要的麦子数量
        :param grid: 格子数
        :return:
        '''
        # f(1) = 1
        wheat_numbers = 1
    
        sums = wheat_numbers
        for i in range(2, grid+1):
            wheat_numbers *= 2
            sums += wheat_numbers
    
        print('when grid = %d, wheats numbers = %d' % (grid, sums))
    
        return sums
    

    简单的测试例子:

    if __name__ == '__main__':
        print('compute numbers of wheat!')
        numbers_grid = 63
        get_number_of_wheat(numbers_grid)
        print('finish')
    

    给定格子数量是 63 个,输出结果如下:

    compute numbers of wheat!
    when grid = 63, wheats numbers = 9223372036854775807
    finish
    

    所以这个天文数字是 19 位数–9223372036854775807,真的是非常的多!假设一袋 50 斤的麦子估计有 130 万粒麦子,那么这个计算结果是相当于 70949 亿袋 50 斤的麦子!

    迭代法的应用

    看完上述例子,相信应该对迭代法的基本概念比较了解了,而迭代法的基本步骤也很简单,分为三个步骤:

    • 确定用于迭代的变量。上述例子中,这个迭代变量就是f(n)f(n-1)
    • 建立迭代变量之间的递推关系。上述例子中,这个递归关系是f(n)=f(n-1)*2
    • 控制迭代的过程。这里需要确定迭代的初始条件和终止条件,上述例子,初始条件就是f(1)=1,而终止条件就是达到给定的格子数了。

    那么迭代法有什么应用呢?

    其实,它在数学和计算机领域都有很广泛的应用,如:

    • 求数值的精确或者近似解。典型的方法包括二分法(Bisection method)和牛顿迭代法(Newton’s method);
    • 在一定范围内查找目标值。典型方法包括二分查找,其实也是二分法在搜索方面的应用;
    • 机器学习算法中的迭代。比如 Kmeans 聚类算法(不断迭代来对数据进行聚类)、马尔科夫链(Markov chain)、梯度下降法(Gradient descent)等。迭代法在机器学习中有广泛的应用,其实是因为机器学习的过程,就是根据已知数据和一定的假设,求一个局部最优解。迭代法可以帮助学习算法逐步搜索,直到发现这种解。

    接下来会重点介绍求数值的解和查找匹配记录,这两个应用其实都是采用二分法来实现。

    求方程的精确或者近似解

    迭代法除了用于计算庞大的数字,还可以帮助我们进行无穷次地逼近,求得方程的精确或者近似解

    举个例子,我们要计算一个给定的正整数n(n>1)的平方根,并且不能采用编程语言自带的函数,应该如何计算呢?

    首先我们可以明确的是,对于给定的正整数n,它的平方根肯定是小于它,但大于1,也就是这个平方根的取值范围是 1 到 n ,在这个范围内求一个数值的平方等于n

    这里就可以通过采用刚刚说的二分法。每次查看区间内的中间值,检查它是否符合标准。

    比如我们要求 10 的平方根,寻找的区间就是[1,10],第一个中间值就是(1+10)/2=11/2=5.5,而 5.5 的平方等于 30.25,明显比 10 大,所以寻找区间变成 5.5 的左侧,也就是[1, 5.5],中间值就是 3.25,但 3.25 的平方是 10.5625,依然大于 10,寻找区间变为[1, 3.25],中间值变为 2.125, 2.125 的平方是 4.515625,小于 10,所以区间就是[2.125, 3.25],这样继续寻找和计算中间值的平方,直到发现某个数的平方正好是 10。

    具体步骤如下图:

    这里用代码实现,如下图所示:

    def get_square_root(n, threshold, max_try):
        '''
        计算大于 1 的正整数的平方根
        :param n: 给定正整数
        :param threshold: 误差的阈值
        :param max_try: 最大尝试次数
        :return:
        '''
        if n <= 1:
            return -1.0
        # interval boundary 区间的左右边界
        left = 1.0
        right = float(n)
        for idx in range(max_try):
            # 防止溢出
            middle = left + (right - left) / 2
            square = middle * middle
            # 误差
            delta = abs(square / n - 1)
            if delta <= threshold:
                return middle
            else:
                if square > n:
                    right = middle
                else:
                    left = middle
    
        return -2.0
    

    简单的测试例子:

    square_root = get_square_root(10, 0.000001, 10000)
    if square_root == -1.0:
        print('please input a number > 1')
    elif square_root == -2.0:
        print('cannot find the square root')
    else:
        print('square root==', square_root)
    

    输出结果是:

    square root== 3.1622767448425293
    

    这里代码中,设置了两个控制迭代结束的参数:

    1. threshold:误差的阈值,用于控制解的精度。理论上二分法可以通过无限次迭代求到精确解,但实际应用还需要考虑时间和计算资源,所以一般我们只需要一个近似解,而不需要完全精确的数据;
    2. max_try:控制迭代的次数。设置这个参数也是为了避免使用while True循环可能导致的死循环,当然理论上设置了threshold是可以避免死循环的,但这是一个良好的编程习惯,主动避免产生的可能性。
    查找匹配记录

    二分法通过迭代式逼近,不仅可以求得方程的近似解,还可以帮助查找匹配的记录

    这里老师给的例子是在自然语言处理中,处理同义词或者近义词的扩展问题。这时,你是会有一个词典,用于记录每个单词的同义词或者近义词。对于一个待查找单词,我们需要在字典找到这个单词,以及对应的所有同义词和近义词,然后进行拓展,例如对于单词–西红柿,它的同义词包括了番茄tomato

    词典如下表格所示:

    词条 同义词1 同义词2 同义词3
    西红柿 番茄 tomato

    当处理文章的时候,遇到“西红柿”这个单词,就在字典里查找,返回“番茄”和“tomato"等同义词或者近义词,并添加到文章作为同义词/近义词的拓展。

    这里要解决的问题就是如何在字典查询匹配单词的问题。一种做法就是哈希表。而如果不用哈希表的方法,还可以采用二分查找法。二分查找法进行字典查询的思路如下:

    1. 对整个字典先进行排序(假设是从小到大)。二分法的一个关键前提条件就是所查找区间必须是有序的,这样每次折半的时候,可以知道是往左还是右继续查找。
    2. 使用二分法逐步定位到被查找的单词。同样是每次都选择查找区间的中间值,判断是否和待查找单词一致,如果一致就返回;如果不一致,就进行判断大小,如果比待查找单词小,就需要往中间值右边区间查找;否则就在左边区间查找。
    3. 重复第二步操作,迭代式查找,直到找到单词,或者没有找到,就返回不存在。

    相比于利用二分法查找方程解,二分查找必须要求数据是有序的!

    用代码实现如下:

    def search_word(dictionary, word):
        '''
        查找匹配单词
        :param dictionary: 排序后的字典
        :param word:待查找单词
        :return:
        '''
        if dictionary is None:
            return False
        if len(dictionary) < 1:
            return False
    
        left = 0
        right = len(dictionary) - 1
        while left <= right:
            middle = int(left + (right - left) / 2)
            if dictionary[middle] == word:
                return True
            else:
                if dictionary[middle] > word:
                    right = middle - 1
                else:
                    left = middle + 1
    
        return False
    
    

    简单的测试代码:

    print('find word in dictionary')
    dict_list = ['i', 'am', 'coder']
    dict_list = sorted(dict_list)
    print('sorted dict:', dict_list)
    word_to_find = 'am'
    found = search_word(dict_list, word_to_find)
    if found:
        print('word "%s" found in dictionary--%s!' % (word_to_find, dict_list))
    else:
        print('cannot find the word "%s"' % word_to_find)
    

    输出结果:

    find word in dictionary
    sorted dict: ['am', 'coder', 'i']
    word "am" found in dictionary--['am', 'coder', 'i']!
    finish
    

    迭代法的介绍就到这里了!上述源代码地址:

    https://github.com/ccc013/CodesNotes/blob/master/Maths/lesson_iterations.py


    欢迎关注我的微信公众号–机器学习与计算机视觉,或者扫描下方的二维码,大家一起交流,学习和进步!


    往期精彩推荐

    数学学习笔记
    学习笔记
    Github项目 & 资源教程推荐
    展开全文
  • Jacobi迭代法与Gauss-Seidel迭代法

    万次阅读 多人点赞 2016-01-25 19:29:58
    今天刚好有朋友和我讨论泊松图像融合算法,我说我过去文章里给出的是最原始、最直观...但是过去我浅尝辄止了,也没深究,今天刚好再提到,小看了一下,似乎涉及高斯-塞德尔迭代法。这是一种常用的基于迭代的数值方法

    之前我在博客里介绍过牛顿-拉弗逊迭代法,对数据挖掘技术熟悉的同学应该还知道有梯度下降法(其实也是一种迭代算法)。今天刚好有朋友和我讨论泊松图像融合算法,我说我过去文章里给出的是最原始、最直观的实现算法。对于理解泊松融合的原理比较有帮助,但是效率可能并不理想。印象中,泊松融合是有一个以矩阵为基础的快速算法的。但是过去我浅尝辄止了,也没深究,今天刚好再提到,小看了一下,似乎涉及高斯-塞德尔迭代法。好吧,博主君暂且把知道的这部分内容做个介绍吧。


    一、雅各比迭代法


    考虑如下方程组:
    4xy+z=74x8y+z=212x+y+5z=15 4x-y+z=7\\4x-8y+z=-21\\-2x+y+5z=15
    上述方程可表示成如下形式:
    x=7+yz4y=21+4x+z8z=15+2xy5 x=\frac{7+y-z}{4} \\y=\frac{21+4x+z}{8}\\z=\frac{15+2x-y}{5}
    这样就提出了下列雅可比迭代过程:

    xk+1=7+ykzk4yk+1=21+4xk+zk8(3)zk+1=15+2xkyk5 \begin{aligned} x_{k+1} & =\frac{7+y_k-z_k}{4} \\ y_{k+1} & =\frac{21+4x_k+z_k}{8} & (3) \\ z_{k+1} & =\frac{15+2x_k-y_k}{5} \end{aligned}

    如果从P0=(x0,y0,z0)=(1,2,2)P_0=(x_0,y_0,z_0)=(1,2,2)开始,则上式中的迭代将收敛到解(2,4,3)(2,4,3)

    x0=1,y0=2x_0=1,y_0=2z0=2z_0=2代入上式中每个方程的右边,即可得到如下新值:
    x1=7+224=1.75y1=21+4+28=3.375z1=15+225=3.00 x_1=\frac{7+2-2}{4}=1.75\\y_1=\frac{21+4+2}{8}=3.375\\z_1=\frac{15+2-2}{5}=3.00

    新的点P1=(1.75,3.375,3.00)P_1=(1.75,3.375,3.00)P1P_1更接近(2,4,3)(2,4,3)。使用迭代过程(3)生成点的序列{P0P_0}将收敛到解(2,4,3)(2,4,3)

    这里写图片描述

    这个过程称为雅可比迭代,可用来求解某些类型的线性方程组。从上表中可以看出,经过19步选代,选代过程收敛到一个精度为9 位有效数字的近似值(2.00000000, 4.00000000, 3.00000000)。但有时雅可比迭代法是无效的。通过下面的例子可以看出,重新排列初始线性方程组后,应用雅可比迭代法可能会产生一个发散的点的序列。

    设重新排列的线性方程组如下:
    2x+y+5z=154x8y+z=214xy+z=7 -2x+y+5z=15\\4x-8y+z=-21\\4x-y+z=7

    这些方程可以表示为如下形式:
    x=15+y5z2y=21+4x+z8z=74x+y5 x=\frac{-15+y-5z}{2} \\y=\frac{21+4x+z}{8}\\z=\frac{7-4x+y}{5}

    这可以用如下雅可比迭代过程求解:
    \begin{aligned}
    x_{k+1} & =\frac{-15+y_k+5z_k}{2} \
    y_{k+1} & =\frac{21+4x_k+z_k}{8} \
    z_{k+1} & =7-4x_k+y_k
    \end{aligned}

    如果从P0=(x0,y0,z0)=(1,2,2)P_0=(x_0,y_0,z_0)=(1,2,2)开始,则上式中的迭代将对解(2,4,3)(2,4,3)发散。将x0=1x_0=1y0=2y_0=2z0=2z_0=2带入上式中每个方程的右边,即可得到新值x1x_1y1y_1z1z_1

    x1=15+2+102=1.5y1=21+4+28=3.375z1=74+2=5.00 x_1=\frac{-15+2+10}{2}=-1.5\\y_1=\frac{21+4+2}{8}=3.375\\z_1=7-4+2=5.00

    新的点P1=(1.5,3.375,5.00)P_1=(-1.5,3.375,5.00)P0P_0更远地偏离(2,4,3)(2,4,3)。使用上述迭代过程生成点的序列是发散的。

    这里写图片描述


    二、高斯-塞德尔迭代法


    有时候通过其他方面可以加快迭代的收敛速度。观察由雅可比迭代过程(3)产生的三个序列{xkx_k},{yky_k}和{zkz_k},它们分别收敛到2,4和3。由于xk+1x_k+1被认为是比xkx_k更好的xx之近似值,所以在计算yk+1y_k+1时用xk+1x_k+1来替换xkx_k是合理的。同理,可用xk+1x_k+1yk+1y_k+1计算zk+1z_k+1。下面的例子演示了对上述例子中给出的方程组使用上述方法的情况。

    设给定上述线性方程组并利用高斯-塞德尔(Gauss-Seidel)迭代过程求解:

    xk+1=7+ykzk4yk+1=21+4xk+1+zk8zk+1=15+2xk+1yk+15x_{k+1}=\frac{7+y_k-z_k}{4}\\ y_{k+1}=\frac{21+4x_{k+1}+z_k}{8}\\ z_{k+1}=\frac{15+2x_{k+1}-y_{k+1}}{5}

    如果从P0=(x0,y0,z0)=(1,2,2)P_0=(x_0,y_0,z_0)=(1,2,2)开始,用上式中的迭代可收敛到解(2,4,3)(2,4,3)
    y0=2y_0=2z0=2z_0=2代入上式第一个方程可得

    x1=7+224=1.75x_1=\frac{7+2-2}{4}=1.75

    x1=1.75x_1=1.75z0=2z_0=2代入第二个方程可得

    y1=21+4×1.75+28=3.75y_1=\frac{21+4\times1.75+2}{8}=3.75

    x1=1.75x_1=1.75y1=3.75y_1=3.75代入第三个方程可得

    z1=15+2×1.753.755=2.95z_1=\frac{15+2\times1.75-3.75}{5}=2.95

    新的点P1=(1.75,3.75,2.95)P_1=(1.75,3.75,2.95)P0P_0更接近解(2,4,3)(2,4,3),而且比之前例子中的值更好。用迭
    代(7)生成序列{PkP_k收敛到(2,4,3)(2,4,3)

    这里写图片描述

    正如前面讨论的,应用雅各比迭代法计算有时可能是发散的。所以有必要建立一些判定条件来判断雅可比迭代是否收敛。在给出这个条件之前,先来看看严格对角占优矩阵的定义。
    设有NN×\timesNN维矩阵AA,如果

    这里写图片描述
    其中,ii是行号,jj是列号,则称该矩阵是严格对角占优矩阵。显然,严格对角占优的意思就是指对角线上元素的绝对值不小于所在行其他元素的绝对值和。

    \begin{aligned}
    a_{11}x_1+a_{12}x_2+…+a_{1j}x_j+…+a_{1N}x_N& =b_1 \
    a_{21}x_1+a_{22}x_2+…+a_{2j}x_j+…+a_{2N}x_N& =b_2 \
    … … \
    a_{j1}x_1+a_{j2}x_2+…+a_{jj}x_j+…+a_{jN}x_N& =b_j \
    … … \
    a_{N1}x_1+a_{N2}x_2+…+a_{Nj}x_j+…+a_{NN}x_N& =b_N
    \end{aligned}

    设第k点为Pk=(x1(k),x2(k),...,xj(k),...,xN(k),)P_k=(x_1^{(k)},x_2^{(k)},...,x_j^{(k)},...,x_N^{(k)},),则下一点为Pk+1=(x1(k+1),x2(k+1),...,xj(k+1),...,xN(k+1),)P_{k+1}=(x_1^{(k+1)},x_2^{(k+1)},...,x_j^{(k+1)},...,x_N^{(k+1)},)。向量PkP_k的上标(k)(k)可用来标识属于这一点的坐标。迭代公式根据前面的值(x1(k),x2(k),...,xj(k),...,xN(k),)(x_1^{(k)},x_2^{(k)},...,x_j^{(k)},...,x_N^{(k)},),使用上述线性方程组中第jj行求解xj(k+1)x_j^{(k+1)}

    雅可比迭代:

    xj(k+1)=bjaj1x1(k)...ajj1xj1(k)ajj+1xj+1(k)...ajNxN(k)ajjx_j^{(k+1)}=\frac{b_j-a_{j1}x_1^{(k)}-...-a_{jj-1}x_{j-1}^{(k)}-a_{jj+1}x_{j+1}^{(k)}-...-a_{jN}x_N^{(k)}}{a_{jj}}

    其中j=1,2,...,Nj=1,2,...,N
    雅可比迭代使用所有旧坐标来生成所有新坐标,而高斯-塞德尔迭代尽可能使用新坐标得到更新的坐标。

    高斯-塞德尔迭代:

    xj(k+1)=bjaj1x1(k+1)...ajj1xj1(k+1)ajj+1xj+1(k)...ajNxN(k)ajjx_j^{(k+1)}=\frac{b_j-a_{j1}x_1^{(k+1)}-...-a_{jj-1}x_{j-1}^{(k+1)}-a_{jj+1}x_{j+1}^{(k)}-...-a_{jN}x_N^{(k)}}{a_{jj}}

    其中j=1,2,...,Nj=1,2,...,N

    下面的定理给出了雅可比迭代收敛的充分条件。
    (雅可比选代) 设矩阵AA具有严格对角优势,则AX=BAX=B有惟一解X=PX = P。利用前面给出的迭代式可产生一个向量序列{PkP_k},而且对于任意初始向量P0P_0,向量序列都将收敛到PP

    当矩阵AA具有严格对角优势时,可证明高斯-塞德尔迭代法也会收敛。在大多数情况下,高斯-塞德尔迭代法比雅可比迭代法收敛得更快,因此通常会利用高斯-塞德尔迭代法。但在某些情况下,雅可比迭代会收敛,而高斯-塞德尔迭代不会收敛。

    特别说明:以上内容主要取材自《数值方法(MATLAB版)(第四版)》马修斯等著,电子工业出版社2010年出版发行。

    展开全文
  • 用二分法和牛顿迭代法都可以解决。题目: Strange fuction Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others) Total Submission(s): 739 Accepted Submission(s): 5
  • 数学知识】非线性方程求解的二分法以及牛顿迭代法导入包二分法问题1编程作业牛顿迭代法编程作业 本博客不谈及理论推导,只提供代码实现,关于理论推导,大家可以查看其它博客文章。 导入包 import sys import math...
  • 1、简单迭代法: 参考代码: #include<stdio.h> #include<math.h> #define x0 3//初值 #define MAXREPT 1000//迭代次数 #define EPS 0.5E-6//精度 #define G(x) pow(sin(x)+12*x-1,1.0/3)//你的迭代函数...
  • [导读] 前面刚转了一篇文章提到了牛顿-拉夫逊(拉弗森)(Newton-Raphson method)方法,感觉这个数学方法很有必要相对深入写一篇文章来总结分享印证一下自己的理解。这是写本文的由来,如果发现文章中有错误之处,...
  • Jacobi迭代法和Gauss-Seidel迭代法.docx
  • 谈谈常见迭代优化方法

    万次阅读 2016-02-02 09:54:22
    本文主要介绍三种方法,分别是梯度下降,共轭梯度(Conjugate Gradient Method)和近似牛顿(Quasi-Newton)。具体在stanford-nlp中都有对应的实现,由于前两种方法都涉及到梯度的概念,我们首先从介绍梯度开始。 ...
  • 余数在生活中非常常见,日历,日期,其中包含一个定理 同余定理。同余定理:余数总是在一个固定的范围内,任意两个整数 a 和 b,如果它们除以正整数 m 得到的余数相等,我们就称 a 和 b 对于模 m 同余。同余定理主要...
  • 递归把计算交给计算机,归纳把计算交给⼈,前者是拿...数学归纳:证明基本情况(n=1)是否成立,假设n=k-1成立,证明n=k也是成立。 递归调用和数学归纳的逻辑是一样的! 数学归纳正确,递归调用逻辑也正确。 ...
  • 昨天菜鸟自己为自己放了一天假,所以没有更新,也没有看学习内容(QωQ),望大家理解!...文章目录程序员的数学基础课1、到底什么是迭代法?2、迭代法有什么具体应用?相信大家有点不解,为什么二分法是迭代法
  • 迭代法

    千次阅读 2019-01-09 14:36:50
      刚认识女朋友时候,那段时间很流行微信红包(好吧,一直都很流行),每天早上发...  因为那将是很恐怖的后果,这里面就包含着我们这次要讲的重点:迭代法(Iterative Method)。 迭代法定义   那么何为迭代...
  • 雅克比迭代法

    2017-12-06 18:25:28
    考虑线性方程组Ax = b时,一般当A为低阶稠密矩阵时,用主元消去法解此方程组是有效方法。但是,对于由工程技术中产生...雅克比迭代法就是众多迭代法中比较早且较简单的一种,其命名也是为纪念普鲁士著名数学家雅可比。
  • 开个坑 有空补
  • 算法--迭代法

    万次阅读 多人点赞 2018-10-24 18:52:32
    迭代法(Iteration)是一种不断用变量的旧值递推出新值的解决问题的方法。迭代算法是用计算机解决问题的一种基本方法,一般用于数值计算。累加、累乘都是迭代算法的基础应用。典型案例:牛顿迭代法”。 步骤: 确定...
  • 迭代法Ⅰ 象棋&米粒-迭代法Ⅱ 编程求麦粒的数量Ⅲ 迭代法的应用① 求方程的精确或者近似解② 查找匹配记录 Ⅰ 象棋&米粒-迭代法 这里引用一个大家从小就听过的小故事。 传说,印度的舍罕国王打算重赏国际...
  • 关注、星标嵌入式客栈,精彩及时送达[导读] 前面刚转了一篇文章提到了牛顿-拉夫逊(拉弗森)(Newton-Raphson method)方法,感觉这个数学方法很有必要相对深入写一篇文章来...
  • 牛顿迭代法

    2020-06-16 07:53:30
    牛顿迭代法 方法讲解:https://blog.csdn.net/ccnt_2012/article/details/81837154(马同学高等数学) 摘要 迭代公式: Xn+1=Xn−f(Xn)f′(Xn) X_{n+1}=X_n-\frac{f(X_n)}{f'(X_n)} Xn+1​=Xn​−f′(Xn​)f(Xn​)...
  • 文章目录一、点迭代法二、点迭代的步骤 一、点迭代法 MATLAB 作图显示: y1 = @(x) exp(x)/3; % fplot 绘制表达式函数 fplot(y1,[0,1],'b') hold on % 画出参考线,观察 y2= @(x) x; fplot(y2,[0,1],'r') 二、点...
  • 随后用迭代y=(y+x/y)/2 进行迭代循环求出x的平方根y。 二、 求x的立方根,先假定初始值为y 随后用跌打y=(2y+x/y^2)/2 进行迭代循环求出x的立方根y。   迭代次数越多,越接近真实值。误差=|y^2-x| ...
  • 西京学院数学软件实验任务书课程名称数学软件实验班级数0901学号0912020107姓名李亚强实验课题线性方程组的J-迭代,GS-迭代,SOR-迭代方法。实验目的熟悉线性方程组的J-迭代,GS-迭代,SOR-迭代方法。实验要求运用...
  • D. M....由于超松弛迭代法公式简单,编制程序容易,很多工程学、计算数学中都会应用超松弛迭代方法。使用超松弛迭代法的关键在于选取合适的松弛因子,如果松弛因子选取合适,则会大大缩短计算时间。
  • #include #include using namespace std; double f(double x); double f1(double x); double f2(double x); int main() { double x0, x1,e; e = 0.0001; //终止门限 x0 = 15;
  • 该方法称为牛顿迭代方法。2.条件f(x)函数是连续可导函数。f(x)在局部收敛,当时,局部收敛。注意:牛顿迭代法的局部收敛性,很依赖于初始值的取法。也就是说,初始值的选取,决定该区域的收敛性。3.思想其总思想还是...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 104,649
精华内容 41,859
关键字:

常见的数学迭代方法