Matlab非线性方程数值解法

实验目的

用Matlab实现非线性方程的二分法、不动点迭代法

实验要求

1. 给出二分法算法和不动点迭代算法

2. 用Matlab实现二分法

3. 用Matlab实现不动点迭代法

实验内容

（1）在区间[0,1]上用二分法和不动点迭代法求的根到小数点后六位。

（2）二分法的基本思想：逐步二分区间[a,b]，通过判断两端点函数值的符号，进一步缩小有限区间，将有根区间的长度缩小到充分小，从而，求得满足精度要求的根的近似值。

（3）不动点迭代法基本思想：已知一个近似根，构造一个递推关系（迭代格式），使用这个迭代格式反复校正根的近似值，计算出方程的一个根的近似值序列，使之逐步精确法，直到满足精度要求（该序列收敛于方程的根）。

实验步骤

（1）二分法算法与MATLAB程序（二分法的依据是根的存在性定理，更深地说是介值定理）。

 

 

 

MATLAB程序，

%二分法
%输入：f(x)=0的f(x)，[a,b]的a,b，精度ep
%输出：近似根root，迭代次数k
function [root,k]=bisect(fun,a,b,ep)
if nargin>3
    elseif nargin<4
        ep=1e-5;%默认精度
    else
        error('输入参数不足');%输入参数必须包括f(x)和[a,b]
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

运行示例（并未对输出格式做控制，由于精度要求，事后有必要控制输出的精度）：

 优化代码，减小迭代次数（在迭代前，先搜寻更适合的有根区间）

%二分法改良
%在一开始给定的区间中寻找更小的有根区间
%输入：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

运行示例（同样没有控制输出）

明显地，迭代次数减小许多。

继续优化代码，以获得区间里所有的根（关键是对有根区间的判断与记忆）

%二分法改良
%在一开始给定的区间中寻找更小的有根区间
%输入：f(x)=0的f(x)，[a,b]的a,b，精度ep
%输出：近似根root，迭代次数k
%例子：
%fun=inline('x.^3-x.^2-4*x+4');a=-3;b=3;ep=1e-5;
%fun=inline('x.^2-1');a=-2;b=2;ep=1e-5;
%fun=inline('-x.^2+1');a=-2;b=2;ep=1e-5;
%fun=inline('x.^2-3*x+2');a=0;b=3;ep=1e-5;
function [root,k,x,y]=bisect4(fun,a,b,ep)
if nargin>3
    elseif nargin<4
        ep=1e-5;%默认精度
    else
        error('输入参数不足');%输入参数必须包括f(x)和[a,b]
end
%定义划分区间的分数
divQJ=100;
%等分区间
tX=linspace(a,b,divQJ);
%计算函数值
tY=fun(tX);
%找到函数值的正负变化的位置
locM=find(tY<0);
locP=find(tY>0);
%得到的符号序列，找其中的间断点
z=0;
k=1;
for i=1:length(locM)-1%扫描-符号序列
    %滚筒算法
    if(i==1 && locM(1)~=1)
        z=z+1;
        x(z)=locM(i);
    end
    tmp=locM(i+1)-locM(i);
    if(tmp~=1)
        z=z+1;
        x(z)=locM(i);
        z=z+1;
        x(z)=locM(i+1);
    end
end

l=0;
for i=1:length(locP)-1%扫描+符号序列
    %滚筒算法
    if(i==1 && locP(1)~=1)
        l=l+1;
        y(l)=locP(i);
    end
    tmp=locP(i+1)-locP(i);
    if(tmp~=1)
        l=l+1;
        y(l)=locP(i);
        l=l+1;
        y(l)=locP(i+1);
    end
end
if(z==l-1)
    z=z+1;
    x(z)=locM(end);
elseif(z==l+1)
    l=l+1;
    y(l)=locP(i);
end

if(z==0 && l==0)
    if(locM(end)<locP(1))
        x(1)=locM(end);
        y(1)=locP(1);
    else
        x(1)=locM(1);
        y(1)=locP(end);
    end
    z=1;
elseif(z==0)
    x(1)=locM(1);
    x(2)=locM(end);
    z=2;
elseif(l==0)
    y(1)=locP(1);
    y(2)=locP(end);
    z=2;
end
for i=1:z
a=tX(x(i));
b=tX(y(i));
if fun(a)*fun(b)>0%输入的区间要求
    root=[fun(a),fun(b)];
    k=0;
    return;
end
if(a>b)
    tmp=b;
    b=a;
    a=tmp;
end
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(i)=(a+b)/2;
end
end

运行示例（仍然没有控制输出，笔者忽然觉得自己中了月亮的毒！！！）

 

 

 

 

 使用二分法解得在[0,1]上的根为0.090526。（终于在最后回答问题时做了<_>     控制输出的方法：在调用函数前，在控制台输入format long;得到一long型的数据，根据精度读取。）

（2）不动点迭代算法与MATLAB程序。

 

%不动点迭代法
%输入：fun函数,初始值x0，精度ep，最大迭代次数Nmax
%输出：近似根root，迭代次数k
function [root, k]=iterate(fun,x0,ep,Nmax)
if nargin<4
    Nmax=500;
end
if nargin<3
    ep=1e-5;
end
x=x0;
x0=x+2*ep;
k=0;
while abs(x0-x)>ep && k<Nmax
    x0=x;
    x=fun(x0);
    k=k+1;
end
root=x;
if k==Nmax
    warning('已达到迭代次数上限');
end
end

运行示例（对迭代方程的简单推导得到：）：

 就没优化的二分法与不动点迭代法比较，明显地迭代法的迭代次数更少。

使用不动点迭代法，求得近似根0.090525。

 小结

就上述两个程序而言，二分法的编程更具挑战。在优化二分法的程序过程中，发现二分法的关键就在于对所给定区间的讨论划分。得到一个能求解给定区间里的所有零点的程序的关键是（分类讨论的思想），对符号变化的记录与使用，换言之，解决这个问题的关键就在于分支结构的使用。

（不可否认，第三个二分法程序仍然有bug！！！待解决——有空再说。）