Matlab非线性方程数值解法（2）

实验目的

　　用Matlab实现非线性方程的Aitken加速法、牛顿法、简化牛顿法和牛顿下山法。

实验要求

1. 　　给出Aitken加速法、牛顿法、简化牛顿法和牛顿下山法算法。
2. 　　 用Matlab实现Aitken加速法、牛顿法、简化牛顿法和牛顿下山法。

实验内容

 实验步骤

（1）Aitken加速算法，

　　MATLAB编程，

%Aitken算法
%迭代函数fun，初始值x，精度ep，最大迭代次数Nmax。
%满足精度（或已达到最大迭代次数）的近似根root，迭代次数k。
function [root,k]=Aitken(fun,x,ep,Nmax)
%默认迭代次数和默认精度
if nargin<4
    Nmax=500;
end
if nargin<3
    ep=1e-5;
end
%初始化
k=0;
r(1,1)=x;
%迭代
while abs(x-fun(x))>ep
    xk=x;
    x=fun(x);
    x1=x;
    x=fun(x);
    x2=x;
    x=x2-(x2-x1)^2/(x2-2*x1+xk);
    k=k+1;
    r(k+1,1)=x;
end
%root=x;
root=r(:,1);
if k==Nmax
    warning('已迭代上限次数');
end
end

　　运行示例（为了减小篇幅，初始值的选取由下面的讨论得到），

　　

　　当初始值为10时，节选截图,

　　

 　　可见该算法迭代了3610次，这个结果的好坏需要与其他算法对比，使用上一个实验所编写的bisect3.m（改进的二分法）：

　　

 　　可见迭代20次，已经接近所要求的近似根了。

　　为了了解该算法的特点，本报告对算法进行修改：缩小区间，再使用Aitken算法迭代。分析：易知上述到Aitken算法的时间复杂度与精度有关（O(g(e))，e的变化受迭代函数影响），初始值越接近所要求的近似根，运算的步骤越少。所以对初始值做处理是有必要的，这里如果使用20次2分法预报初值时，精度可以达到0.5^20，约1e-6。

%Aitken1算法
%迭代函数fun，初始值x，精度ep，最大迭代次数Nmax。
%满足精度（或已达到最大迭代次数）的近似根root，迭代次数k。
function [root,k]=Aitken1(fun,n,ep,Nmax)
%默认迭代次数和默认精度
if nargin<4
    Nmax=500;
end
if nargin<3
    ep=1e-5;
end
[x,k]=bisect3(fun,0,n,1e-4);
if abs(fun(x))<ep
    root=x;
    k=0;
    return;
end
%初始化
r(1,1)=x;
%迭代
while abs(x-fun(x))>ep
    xk=x;
    x=fun(x);
    x1=x;
    x=fun(x);
    x2=x;
    x=x2-(x2-x1)^2/(x2-2*x1+xk);
    k=k+1;
    r(k+1,1)=x;
end
%root=x;
root=r(:,1);
if k>=Nmax
    warning('已迭代上限次数');
end
end

%二分法改良
%在一开始给定的区间中寻找更小的有根区间
%输入：f(x)=0的f(x)，[a,b]的a,b，精度ep
%输出：近似根root，迭代次数k
%得到的根是优化区间里的最大根
function [root,k]=bisect3(fun,a,b,ep)
if nargin>3
    elseif nargin<4
        ep=1e-5;%默认精度
    else
        error('输入参数不足');%输入参数必须包括f(x)和[a,b]
end
%定义划分区间的分数
divQJ=1000;
%等分区间
tX=linspace(a,b,divQJ);
%计算函数值
tY=fun(tX);
%找到函数值的正负变化的位置
locM=find(tY<0);
locP=find(tY>0);
%定义新区间
if tY(1)<0
a=tX(locM(end));
b=tX(locP(1));
else
a=tX(locP(end));
b=tX(locM(1));
end
if fun(a)*fun(b)>0%输入的区间要求
    root=[fun(a),fun(b)];
    k=0;
    return;
end
k=1;
while abs(b-a)/2>ep%精度要求
    mid=(a+b)/2;%中点
    if fun(a)*fun(mid)<0
        b=mid;
    elseif fun(a)*fun(mid)>0
        a=mid;
    else
        a=mid;b=mid;
    end
    k=k+1;
end
root=(a+b)/2;
end

　　运行示例，

　　

 　　0的部分表示调用bisect3.m（精度要求1e-4），在bisect3.m迭代的次数为7，在Aitken部分迭代的次数为5。总共迭代12次远比3610次小。

　　对【题目】的解答，为　　 MATLAB编程，

%牛顿法
%输入：迭代函数fun,迭代函数的导函数dfun，初始值x，精度ep，最大迭代值Nmax。
%输出：近似根root，迭代次数k，迭代成功与否index。
function [root,k]=Newton2(fun,dfun,x0,ep,Nmax)
%默认迭代次数和默认精度
if nargin<5
    Nmax=500;
end
if nargin<4
    ep=1e-5;
end
%初始化
k=0;
x=x0;x0=x+2*ep;
%迭代
while abs(x0-x)>ep
    k=k+1;
    x0=x;
    x=x0-fun(x0)/dfun(x0);
end
root=x;
if k>=Nmax
    warning('已迭代上限次数');
end
end

　　运行示例，

　　

 　　所以解得【题目】，（与上述的改进二分法求得的结果相近）。

 　　MATLAB程序，

%牛顿平行弦法
%输入：迭代函数fname，迭代函数的一阶导函数dfname，初始值x0，精度ep，最大迭代值Nmax。
%输出：近似根root，迭代次数k，迭代成功与否index。
function [root,k,index]=newton4(fname,dfname,x0,ep,Nmax)
%默认迭代次数和默认精度
if nargin<5
    Nmax=500;
end
if nargin<4
    ep=1e-5;
end
%初始化
index=0;k=1;c=1/dfname(x0);x=x0;
while k<Nmax
    x1=x-c*fname(x);
    if abs(x1-x)<ep
        break;
    end
    x=x1;
    k=k+1;
end
if 0<c*dfname(x) && c*dfname(x)<2
    root=x;
    index=1;
else
    root=x;
    index=0;
    disp('不满足收敛条件');
end
if k==Nmax
    warning('已迭代上限次数');
end
end

　　运行示例，

　　

 　　所以解得【题目】，（效果与牛顿法相近）。

　　MATLAB编程，

%牛顿下山法
%输入：迭代函数fname，迭代函数的一阶导函数dfname，初始值x0，精度ep，最大迭代值Nmax。
%输出：近似根root，迭代次数k，%此处不展示函数迭代值单调变化f。
function [root,k]=newton3(fname,dfname,x0,ep,Nmax)
%默认迭代次数和默认精度
if nargin<5
    Nmax=500;
end
if nargin<4
    ep=1e-5;
end
%初始化
k=0;x=x0-fname(x0)/dfname(x0);
%f=abs(fname(x0));
%迭代
while abs(x0-x)>ep && k<Nmax
    k=k+1;
    r=1;
    while abs(fname(x))>abs(fname(x0))
        r=r/2;
        x=x0-r*fname(x0)/dfname(x0);
        if abs(fname(x))<abs(fname(x0))
            break;
        end
    end
    x0=x;
    x=x0-fname(x0)/dfname(x0);
    %f=abs(fname(x0));
end
root=x;
if k==Nmax
    warning('已迭代上限次数');
end
end

　　运行示例，

　　

　　所以，解得【题目】，（效果与牛顿法相差不大）。

 　　总结解得的近似根，

1. 　　,代入原式，1.4355526.
2. 3.4.,代入原式，-1.009655。总体来说，后面3种算法都比第1种算法更高效。 

 小结

　　在解第1题时，本报告对初始值的选取进行了一些讨论。使用二分法锁定耙子固然简单，但是二分法对区间也是有要求，起码是有根区间（在改进二分法中有根区间的限制可以放宽一点），例如，上题我的修改是（0到n）作为瞄准的区间，如果该题的根在负轴上呢？总而言之，初始值的优选区间也是一个值得讨论的问题。第2个值得讨论的问题是，不同算法的混用能否获得更好的结果？