• Matlab非线性方程数值解法


    实验目的

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

    实验要求

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

    2. Matlab实现二分法

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

    实验内容

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

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

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

    实验步骤

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

     

    MATLAB程序,

     1 %二分法
     2 %输入:f(x)=0的f(x),[a,b]的a,b,精度ep
     3 %输出:近似根root,迭代次数k
     4 function [root,k]=bisect(fun,a,b,ep)
     5 if nargin>3
     6     elseif nargin<4
     7         ep=1e-5;%默认精度
     8     else
     9         error('输入参数不足');%输入参数必须包括f(x)和[a,b]
    10 end
    11 if fun(a)*fun(b)>0%输入的区间要求
    12     root=[fun(a),fun(b)];
    13     k=0;
    14     return;
    15 end
    16 k=1;
    17 while abs(b-a)/2>ep%精度要求
    18     mid=(a+b)/2;%中点
    19     if fun(a)*fun(mid)<0
    20         b=mid;
    21     elseif fun(a)*fun(mid)>0
    22         a=mid;
    23     else
    24         a=mid;b=mid;
    25     end
    26     k=k+1;
    27 end
    28 root=(a+b)/2;
    29 end
    二分法1

    运行示例(并未对输出格式做控制,由于精度要求,事后有必要控制输出的精度):

     优化代码,减小迭代次数(在迭代前,先搜寻更适合的有根区间)

     1 %二分法改良
     2 %在一开始给定的区间中寻找更小的有根区间
     3 %输入:f(x)=0的f(x),[a,b]的a,b,精度ep
     4 %输出:近似根root,迭代次数k
     5 %得到的根是优化区间里的最大根
     6 function [root,k]=bisect3(fun,a,b,ep)
     7 if nargin>3
     8     elseif nargin<4
     9         ep=1e-5;%默认精度
    10     else
    11         error('输入参数不足');%输入参数必须包括f(x)和[a,b]
    12 end
    13 %定义划分区间的分数
    14 divQJ=1000;
    15 %等分区间
    16 tX=linspace(a,b,divQJ);
    17 %计算函数值
    18 tY=fun(tX);
    19 %找到函数值的正负变化的位置
    20 locM=find(tY<0);
    21 locP=find(tY>0);
    22 %定义新区间
    23 if tY(1)<0
    24 a=tX(locM(end));
    25 b=tX(locP(1));
    26 else
    27 a=tX(locP(end));
    28 b=tX(locM(1));
    29 end
    30 if fun(a)*fun(b)>0%输入的区间要求
    31     root=[fun(a),fun(b)];
    32     k=0;
    33     return;
    34 end
    35 k=1;
    36 while abs(b-a)/2>ep%精度要求
    37     mid=(a+b)/2;%中点
    38     if fun(a)*fun(mid)<0
    39         b=mid;
    40     elseif fun(a)*fun(mid)>0
    41         a=mid;
    42     else
    43         a=mid;b=mid;
    44     end
    45     k=k+1;
    46 end
    47 root=(a+b)/2;
    48 end
    二分法2

    运行示例(同样没有控制输出)

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

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

      1 %二分法改良
      2 %在一开始给定的区间中寻找更小的有根区间
      3 %输入:f(x)=0的f(x),[a,b]的a,b,精度ep
      4 %输出:近似根root,迭代次数k
      5 %例子:
      6 %fun=inline('x.^3-x.^2-4*x+4');a=-3;b=3;ep=1e-5;
      7 %fun=inline('x.^2-1');a=-2;b=2;ep=1e-5;
      8 %fun=inline('-x.^2+1');a=-2;b=2;ep=1e-5;
      9 %fun=inline('x.^2-3*x+2');a=0;b=3;ep=1e-5;
     10 function [root,k,x,y]=bisect4(fun,a,b,ep)
     11 if nargin>3
     12     elseif nargin<4
     13         ep=1e-5;%默认精度
     14     else
     15         error('输入参数不足');%输入参数必须包括f(x)和[a,b]
     16 end
     17 %定义划分区间的分数
     18 divQJ=100;
     19 %等分区间
     20 tX=linspace(a,b,divQJ);
     21 %计算函数值
     22 tY=fun(tX);
     23 %找到函数值的正负变化的位置
     24 locM=find(tY<0);
     25 locP=find(tY>0);
     26 %得到的符号序列,找其中的间断点
     27 z=0;
     28 k=1;
     29 for i=1:length(locM)-1%扫描-符号序列
     30     %滚筒算法
     31     if(i==1 && locM(1)~=1)
     32         z=z+1;
     33         x(z)=locM(i);
     34     end
     35     tmp=locM(i+1)-locM(i);
     36     if(tmp~=1)
     37         z=z+1;
     38         x(z)=locM(i);
     39         z=z+1;
     40         x(z)=locM(i+1);
     41     end
     42 end
     43 
     44 l=0;
     45 for i=1:length(locP)-1%扫描+符号序列
     46     %滚筒算法
     47     if(i==1 && locP(1)~=1)
     48         l=l+1;
     49         y(l)=locP(i);
     50     end
     51     tmp=locP(i+1)-locP(i);
     52     if(tmp~=1)
     53         l=l+1;
     54         y(l)=locP(i);
     55         l=l+1;
     56         y(l)=locP(i+1);
     57     end
     58 end
     59 if(z==l-1)
     60     z=z+1;
     61     x(z)=locM(end);
     62 elseif(z==l+1)
     63     l=l+1;
     64     y(l)=locP(i);
     65 end
     66 
     67 if(z==0 && l==0)
     68     if(locM(end)<locP(1))
     69         x(1)=locM(end);
     70         y(1)=locP(1);
     71     else
     72         x(1)=locM(1);
     73         y(1)=locP(end);
     74     end
     75     z=1;
     76 elseif(z==0)
     77     x(1)=locM(1);
     78     x(2)=locM(end);
     79     z=2;
     80 elseif(l==0)
     81     y(1)=locP(1);
     82     y(2)=locP(end);
     83     z=2;
     84 end
     85 for i=1:z
     86 a=tX(x(i));
     87 b=tX(y(i));
     88 if fun(a)*fun(b)>0%输入的区间要求
     89     root=[fun(a),fun(b)];
     90     k=0;
     91     return;
     92 end
     93 if(a>b)
     94     tmp=b;
     95     b=a;
     96     a=tmp;
     97 end
     98 while abs(b-a)/2>ep%精度要求
     99     mid=(a+b)/2;%中点
    100     if fun(a)*fun(mid)<0
    101         b=mid;
    102     elseif fun(a)*fun(mid)>0
    103         a=mid;
    104     else
    105         a=mid;b=mid;
    106     end
    107     k=k+1;
    108 end
    109 root(i)=(a+b)/2;
    110 end
    111 end
    二分法3

    运行示例(仍然没有控制输出,笔者忽然觉得自己中了月亮的毒!!!)

     

     

     

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

    (2)不动点迭代算法与MATLAB程序。

     

     1 %不动点迭代法
     2 %输入:fun函数,初始值x0,精度ep,最大迭代次数Nmax
     3 %输出:近似根root,迭代次数k
     4 function [root, k]=iterate(fun,x0,ep,Nmax)
     5 if nargin<4
     6     Nmax=500;
     7 end
     8 if nargin<3
     9     ep=1e-5;
    10 end
    11 x=x0;
    12 x0=x+2*ep;
    13 k=0;
    14 while abs(x0-x)>ep && k<Nmax
    15     x0=x;
    16     x=fun(x0);
    17     k=k+1;
    18 end
    19 root=x;
    20 if k==Nmax
    21     warning('已达到迭代次数上限');
    22 end
    23 end
    不动点迭代法

    运行示例(对迭代方程的简单推导得到:):

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

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

     小结

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

    (不可否认,第三个二分法程序仍然有bug!!!待解决——有空再说。)

  • 相关阅读:
    Go---第七章:接口(小知识点笔记)
    Go---第六章:方法(小知识点笔记)
    Go---第五章:函数(小知识点笔记)
    解决paramiko获取远程脚本延时返回数据的问题
    python字典合并
    关于iperf的使用
    python安装MySQLdb:出错Microsoft Visual C++ 9.0 is required
    v2r
    Win10 HotCorner热角小程序
    去掉显卡桌面右键菜单
  • 原文地址:https://www.cnblogs.com/jianle23/p/12846283.html
Copyright © 2020-2023  润新知