不动点及牛顿法迭代工程实现

　　这篇文章给出不动点及牛顿法迭代是完整实现步骤，没有考虑运行效率及内存损耗，仅供参考。

一、举例计算

求下列方程的实根

(1) . x^2 - 3x + 2 - exp(x) = 0

(2). x^3 + 2x^2 +10x - 20 = 0

　　设计一种不动点及牛顿法迭代法，使迭代序列收敛，且计算到x(k)  - x(k-1)< 10-8为止，结果输出小数点 12 位及以上，通过迭代次数及精度等来比较方法的优劣。

二、不动点及牛顿法迭代实现

2.1 、实现环境

　　MATLAB 2018b

2.2 、实现原理

　　不动点迭代法的实现:通过合理构造j(x) 函数，满足不动点迭代局部收敛的条件：即在真值 x* 的邻域附近使得|j'(x)| < 1，在程序中首先需要大概计算下构造函数的导数，看是否满足条件，以其为后面迭代函数，而在不动点迭代函数MATLAB实现中，编写了一个函数myALG(x0, func, x_real, tol, MaxIter) ，输入参数有初值、构造的j(x) 函数、经过vpasolve 计算得出的真实值、误差容忍值及最大迭代次数参数。将这些参数传输至函数中后，函数通过定义中间变量来存储迭代值，将初值送入至构造的迭代函数后的输出赋值给中间变量，然后依据设置的误差容忍程度来进行循环。

　　而牛顿法迭代的实现需要重新构造j(x) 函数：j(x) = x - f(x)/f'(x)  ，需要对所要计算的方程进行求导等进行构造，而在 MATLAB的实现也是和不动点迭代法类似，都是函数通过定义中间变量来存储迭代值，将初值送入至构造的迭代函数后的输出赋值给中间变量，然后依据设置的误差容忍程度来进行循环。两种方法的具体实现见附录代码。

2.3 、数据测试

方程 1 的实现： x^2 - 3x + 2 - exp(x) = 0

不动点迭代构造j(x) 函数：j(x) = x^2 / 3 + 2 / 3 - exp(x) / 3

　　牛顿法迭代构造j(x) 函数：j(x) = x - (x^2 - 3x + 2 - exp(x) ) /(2x - 3 - exp(x))初始化myALG(2.0,@ myFunc, real_value,1e-8,20) ，之后通过 MATLAB 进行测试，经过vpasolve函数计算方程1，得到方程1解为x* =0.25753028543986076045536730493724，之后进行不动点及牛顿法迭代测试。

不动点迭代得到如下结果:

 

 

牛顿法迭代得到如下结果:

 

 

 

方程 2 的实现： x3 + 2x2 +10x - 20 = 0

 

不动点迭代构造j(x) 函数：j(x) = x + (x^3 + 2x^2 +10x - 20) × (-1/ 20)

牛顿法迭代构造j(x) 函数：j(x) = x - (x^3 + 2x^2 +10x - 20) /(3x^2 + 4x +10)

初始化myALG(1.3@ myFunc, real _ value,1e - 8,15) ，之后通过 MATLAB 进行测试， 经过 vpasolve 函数计算方程 2，得到方程 2有几个不同解， 以解x* =1.3688081078213726352274143300213 为例来进行测试。

不动点迭代得到如下结果:

 

 

 

牛顿法迭代得到如下结果:

 

 

2.4 、数据分析

　　通过构造了合适的不动点及牛顿法迭代，实现了方程 1 与方程 2 的求解，并使之与 xk - xk -1 < 10-8 及结果输出小数点 12 位及以上等要求。经过测试可以看出在一定条件下牛顿法收敛的速度要比不动点迭代法的收敛速度要快，甚至相差一定数量级，这与数值分析课本中对牛顿法及不动点迭代法收敛速度相对应。

附录

不动点迭代代码：

math_test_BDD.m

% syms x
% real_value=vpasolve(x^2-3*x+2-exp(x) == 0, x)   %不动点迭代方程式1
% cal_value=myALG(2.0,@myFunc,real_value,1e-8,15);
%  
% if abs(cal_value - real_value)< 1e-8
%     fprintf('abs(cal_value - real_value)*10^8=%f\n',abs(cal_value - real_value)*1e8);
% end
% abs(cal_value - real_value)


syms x
real_value=vpasolve(x^3 + 2*x^2 + 10*x - 20 == 0, x)   %不动点迭代方程式2
real_value_r = real_value(1)
cal_value=myALG_BDD(1.3,@myFunc_BDD,real_value_r,1e-8,15);
 
if abs(cal_value - real_value_r)< 1e-8
    fprintf('abs(cal_value - real_value_r)*10^8=%f\n',abs(cal_value - real_value_r)*1e8);
end
abs(cal_value - real_value_r)

myALG_BDD.m

function [y] = myALG_BDD(x0,func,x_real,tol,MaxIter)

    xn = x0;
    fprintf('Iter  0: %16.14f\n',x0);
    xnp1 = func(xn);
    fprintf('Iter  1: %16.14f\n',xnp1);
    %criterion = abs(xnp1-xn);
    criterion = abs(xnp1-x_real);
    xn = xnp1;
    Iter = 1;
    %while(criterion>tol)
    while(abs(criterion)>tol)
        xnp1 = func(xn);
        %criterion = abs(xnp1-xn);
        criterion = abs(xnp1-x_real);
        xn = xnp1;
        Iter = Iter + 1;
        fprintf('Iter %2.0d: %16.14f\n',Iter,xnp1);
        if Iter>=MaxIter
            break;
        end
    end
    y = xnp1;
end

myFunc_BDD.m

function [y] = myFunc_BDD(x)
   
    %y = x^2-3*x+2-exp(x);
    %y = x^2/3+2/3-exp(x)/3;    %不动点迭代方程式1
    
    %y = x^3 + 2*x^2 + 10*x - 20    y = -x^3/10 - x^2/5 + 2
    y = x + (x^3 + 2*x^2 + 10*x - 20)*(-1/20);   %不动点迭代方程式2
end

牛顿法迭代代码：

math_test_ND.m

% syms x
% real_value=vpasolve(x^2-3*x+2-exp(x) == 0, x)   %不动点迭代方程式1
% cal_value=myALG(2.0,@myFunc,real_value,1e-8,15);
%  
% if abs(cal_value - real_value)< 1e-8
%     fprintf('abs(cal_value - real_value)*10^8=%f\n',abs(cal_value - real_value)*1e8);
% end
% abs(cal_value - real_value)


% syms x
% real_value=vpasolve(x^3 + 2*x^2 + 10*x - 20 == 0, x)   %不动点迭代方程式2
% real_value_r = real_value(1)
% cal_value=myALG(1.3,@myFunc,real_value_r,1e-8,15);
%  
% if abs(cal_value - real_value_r)< 1e-8
%     fprintf('abs(cal_value - real_value_r)*10^8=%f\n',abs(cal_value - real_value_r)*1e8);
% end
% abs(cal_value - real_value_r)



% syms x
% real_value=vpasolve(x^2-3*x+2-exp(x) == 0, x)   %牛顿法迭代1
% real_value_r = real_value
% cal_value=myALG(2,@myFunc,real_value_r,1e-8,15);
%  
% if abs(cal_value - real_value_r)< 1e-8
%     fprintf('abs(cal_value - real_value_r)*10^8=%f\n',abs(cal_value - real_value_r)*1e8);
% end
% abs(cal_value - real_value_r)


syms x
real_value=vpasolve(x^3 + 2*x^2 + 10*x - 20 == 0, x)   %牛顿法迭代2
real_value_r = real_value(1)
cal_value=myALG_ND(1.3,@myFunc_ND,real_value_r,1e-8,15);
 
if abs(cal_value - real_value_r)< 1e-8
    fprintf('abs(cal_value - real_value_r)*10^8=%f\n',abs(cal_value - real_value_r)*1e8);
end
abs(cal_value - real_value_r)

myALG_ND.m

function [y] = myALG_ND(x0,func,x_real,tol,MaxIter)

    xn = x0;
    fprintf('Iter  0: %16.14f\n',x0);
    xnp1 = func(xn);
    fprintf('Iter  1: %16.14f\n',xnp1);
    %criterion = abs(xnp1-xn);
    criterion = abs(xnp1-x_real);
    xn = xnp1;
    Iter = 1;
    %while(criterion>tol)
    while(abs(criterion)>tol)
        xnp1 = func(xn);
        %criterion = abs(xnp1-xn);
        criterion = abs(xnp1-x_real);
        xn = xnp1;
        Iter = Iter + 1;
        fprintf('Iter %2.0d: %16.14f\n',Iter,xnp1);
        if Iter>=MaxIter
            break;
        end
    end
    y = xnp1;
end

myFunc_ND.m

function [y] = myFunc_ND(x)
   


    %y = x^2-3*x+2-exp(x);
    %y = x^2/3+2/3-exp(x)/3;    %不动点迭代方程式1
    
    %y = x^3 + 2*x^2 + 10*x - 20    y = -x^3/10 - x^2/5 + 2
    %y = x + (x^3 + 2*x^2 + 10*x - 20)*(-1/20);   %不动点迭代方程2

    %y = x-(x^2-3*x+2-exp(x))/(2*x-3-exp(x)); %牛顿法迭代1
    
    y = x - (x^3 + 2*x^2 + 10*x - 20)/(3*x^2+4*x+10);  %牛顿法迭代2
    
    
end