数值分析2-解非线性方程的不动点迭代法，斯蒂芬森(steffensen)迭代法，牛顿法以及简化牛顿法的matlab程序

对于形如f(x)=0的单变量非线性方程，可以构造不同的迭代函数进行迭代求根，以f(x)=x³-x-1为例，我们可以简单的通过等式变形构造出x=x³-1和x=(x+1)^(1/3)两种等价形式；也可以自己增加x的表达式构建x=(1/2)×(x³+x-1），所以同一个非线性方程可以构建无穷多的等价形式，只要最终可以化简成f(x)=0的形式就行；
不动点迭代是最基础的迭代法，其代码如下：
代码块1：

%不动点迭代法1
clc;
clear;
syms x;%定义变量x
f(x)=input('请输入函数表达式,变量用x表示：f(x)=');
g(x)=input('请输入等价形式：x=');  %f(x)=0等价为x=g(x)
x1=input('请选择一个初始值:x1=');
eps=input('请输入停止精度要求:eps=');%|b-xk|<=eps/2
k=1;
x=[x1];  %定义数组x,存储每次迭代产生的xk
k=k+1;%迭代次数,注意k=2时对应第1次迭代
x2=g(x1);%第1次迭代，x2对应数组x中的x（2）
x=[x x2];   %把新产生的x2按顺序装进数组x
T=[1,x1,f(x1);2,x2,f(x2)];   %T初始值为2行3列，随着迭代增加行数
while (abs(x(k)-x(k-1))>eps/2)   %判断是达到求根的精度要求，未达到则执行循环（小于才是达到要求）
    if (k>1000)||(abs(f(x(k)))>1e15)    %1000次还没有达到精度要求的话，报错
        error('迭代出错!或许是迭代公式的问题')
    else
        k=k+1;%迭代次数加一
        xk=vpa(g(x(k-1)),30);%产生新的xk，第一次产生x3；注意：可以用vpa函数，控制每次迭代产生xk的有效数字位数，提升迭代速度
        x=[x xk]; %新产生的xk依次装入数组x中，注意：x(1)=x1,x(2)=x2,x(3)=x3
        T=[T;k,xk,f(xk)];%把此次迭代产生的数据装入存储矩阵T中
    end
end
fprintf('   k            xk              f(xk)\n');%输出表头
disp(vpa(T,10));%控制输出的有效数字为8位，根据精度可以自行设置
fprintf('经过%d次迭代，函数方程根的近似解为：x=%.8f\n',k,T(k,2))%注意迭代次数为k-1次
fprintf('f(x(%d))=%.8f\n',k,T(k,3))

写程序的时候要注意数组和矩阵的下标不能为0，所以迭代计数器k是从1开始的，最后输出时要注意减1；这里循环的判断是由数组x来实现的，我们可以对代码进行精简，去掉数组x，如代码块2所示：
代码块2：不动点迭代法-精简版

%不动点迭代法2
clc;%清屏
clear;%清工作区
syms x;%定义变量x
f(x)=input('请输入函数表达式,变量用x表示：f(x)=');
g(x)=input('请输入等价形式：x=');  %f(x)=0等价为x=g(x)
x0=input('请选择一个初始值:x0=');
eps=input('请输入停止精度要求:eps=');%|b-xk|<=eps/2
k=0;
T=[0,x0,f(x0)];
k=k+1;%k=1,第k次迭代
xk=g(x0);%第1次迭代
T=[T;k,xk,f(xk)];
while (abs(g(xk)-xk)>eps/2)||(abs(f(xk))>1)  %判断是达到求根的精度要求，未达到则执行循环（小于才是达到要求）
    if (k>1000)||(abs(f(xk))>1e15)    %1000次还没有达到精度要求的话，报错
        error('迭代出错!或许是迭代公式的问题')
    else
        k=k+1;%迭代次数加一
        xk=vpa(g(xk),30);%产生新的xk，第一次产生x3；注意：可以用vpa函数，控制每次迭代产生xk的有效数字位数，提升迭代速度
        T=[T;k,xk,f(xk)];%把此次迭代产生的数据装入存储矩阵T中
    end
end
T=[T;k+1,g(xk),f(g(xk))];%再加一行是因为while循环判断用的时k+1次的数据
fprintf('   k            xk              f(xk)\n');%输出表头
disp(vpa(T,10));%控制输出的有效数字为10位，根据精度可以自行设置
fprintf('经过%d次迭代，函数方程根的近似解为：x=%.8f\n',k+1,T(k+2,2))%注意迭代次数为k次
fprintf('f(x(%d))=%.8f\n',k+1,T(k+2,3))%当k做下标时，要加1

这两个思路是有一些不同的哈，代码中需要特别注意的点我在注释里标的已经很清楚了，就不再多说了
函数表达式的输入也可以通过内联函数inline实现，如代码块3所示：
代码块3：不动点迭代法-内联函数版

%不动点迭代法3
clc;
clear;
format short;
s=input('请输入函数表达式：f=','s');
f=inline(s);
m=input('请输入等价形式：x=','s');
g=inline(m); %g为方程的等价形式，f（x）=0等价为x=g（x）
x0=input('请选择一个初始值:x0=');
eps=input('请输入停止精度要求:eps=');    %|b-xk|<=eps/2
k=1;
x=[x0];                                 %定义数组x
k=k+1;
x1=g(x0);
x=[x x1];                               %把新产生的x1装进数组x中
T=[1,x0,f(x0);2,x1,f(x1)];   %T初始值为2行3列，随着迭代增加行数
while abs(x(k)-x(k-1))>eps/2
    if abs(f(x(k)))>1e15||(k>5000)%
        error('迭代出错!或许是迭代公式的问题')
    else
        k=k+1;
        xk=g(x(k-1));
        x=[x xk]; %x(1)=x0,x(2)=x1,x(3)=x3,
        T=[T;k,xk,f(xk)];
    end
end
fprintf('    k               xk               f(xk)\n');
disp(vpa(T,10));
fprintf('经过%d次迭代，函数方程根的近似解为：x=%.8f\n',k-1,T(k,2))

还需要提醒一下就是迭代很有可能不收敛，因为迭代不收敛时方程f(x)的值可能越来越大，最起码不会逼近0，所以可以通过迭代次数k的值以及方程f(x)的来判断是否出现迭代出错；
斯蒂芬森迭代、牛顿法以及简化牛顿法相较于不动点迭代而言，代码的实现思路时一样的，把迭代函数修改一下就行，其代码块分别如下：

代码块4：斯蒂芬森迭代-内联函数版

%斯蒂芬森迭代法1
clc;
clear;
format long;
s=input('请输入函数表达式：f=','s');
f=inline(s);
m=input('请输入等价形式：x=','s');
g=inline(m);                            %g为方程的等价形式，f（x）=0等价为x=g（x）
x1=input('请选择一个初始值:x1=');
eps=input('请输入停止精度要求:eps=');    %|b-xk|<=eps/2
k=1;
y1=g(x1);     %y1=g(x0);
z1=g(y1);     %z1=g(y1);
x=[x1];                                 %定义数组x
y=[y1];                                 %定义数组y
z=[z1];                                 %定义数组z
k=k+1;        %k=2
x2=x1-(y1-x1)^2/(z1-2*y1+x1);
y2=g(x2);     %;
z2=g(y2);     %;
x=[x x2];                               %把新产生的x2装进数组x中
y=[y y2];                               %把新产生的y2装进数组y中
z=[z z2];                               %把新产生的z2装进数组x中
fprintf('    k                   xk                  yk                  zk                     f(xk)\n');
T=[1,x1,y1,z1,f(x1);2,x2,y2,z2,f(x2)];   %T初始值为2行5列，随着迭代增加行数
while abs(x(k)-x(k-1))>eps/2
    if (abs(f(x(k)))>1e15)||(k>50)
        error('迭代出错!或许是迭代初值不对的问题')
    else
        k=k+1;
        xk=x(k-1)-(y(k-1)-x(k-1))^2/(z(k-1)-2*y(k-1)+x(k-1));
        x=[x xk]; %
        yk=g(x(k));
        y=[y yk];
        zk=g(yk);
        z=[z zk];
        T=[T;k,xk,yk,zk,f(xk)];
    end
end
disp(T);
fprintf('经过%d次迭代，函数方程根的近似解为：x=%.8f\n',k-1,T(k,2))

内联函数版都是我复习前为了交作业写的，内联函数也不太好使，还是看下面这一版比较好：
代码块5：斯蒂芬森迭代法

%steffensen迭代法
clc;%清屏
clear;%清工作区
syms x;%定义变量x
f(x)=input('请输入函数表达式,变量用x表示：f(x)=');
g(x)=input('请输入等价形式：x=');  %f(x)=0等价为x=g(x)
x1=input('请选择一个初始值:x1=');
eps=input('请输入停止精度要求:eps=');    %|b-xk|<=eps/2
k=1;
y1=g(x1);     %y1=g(x0);
z1=g(y1);     %z1=g(y1);
x=[x1];       %定义数组x
y=[y1];       %定义数组y
z=[z1];       %定义数组z
k=k+1;        %k=2
x2=x1-(y1-x1)^2/(z1-2*y1+x1);
y2=g(x2);     %;
z2=g(y2);     %;
x=[x x2];                               %把新产生的x2装进数组x中
y=[y y2];                               %把新产生的y2装进数组y中
z=[z z2];                               %把新产生的z2装进数组x中
T=[1,x1,y1,z1,f(x1);2,x2,y2,z2,f(x2)];   %T初始值为2行5列，随着迭代增加行数
while abs(x(k)-x(k-1))>eps/2
    if (abs(f(x(k)))>1e15)||(k>500)
        error('迭代出错!或许是迭代初值不对的问题')
    else
        k=k+1;
        xk=vpa(x(k-1)-(y(k-1)-x(k-1))^2/(z(k-1)-2*y(k-1)+x(k-1)),30);
        x=[x xk]; %
        yk=g(x(k));
        y=[y yk];
        zk=g(yk);
        z=[z zk];
        T=[T;k,xk,yk,zk,f(xk)];
    end
end
fprintf('  k             xk           yk           zk                   f(xk)\n');
disp(vpa(T,10));%控制输出的有效数字为10位，根据精度可以自行设置
fprintf('经过%d次迭代，函数方程根的近似解为：x=%.8f\n',k-1,T(k,2))

注意这一行：

 xk=vpa(x(k-1)-(y(k-1)-x(k-1))^2/(z(k-1)-2*y(k-1)+x(k-1)),30);

要控制一下精度哈，不然直接用分数计算，计算量太大了，算不动，太慢了
也可以考虑精简一下代码把数组x，y，z去掉，但是并不如这样直观，可以自己动手试着改一下，我懒得改了~

牛顿法代码：
代码块6：牛顿法

clc;clear;
format long;
syms x;
f(x)=input('请输入函数表达式,变量用x表示：f(x)=');
x1=input('请选择一个初始值:x1=');
eps=input('请输入停止精度要求:eps='); 
df(x)=diff(f,x);%f(x)的导数
g(x)=f(x)/df(x);
h(x)=x-g(x);%迭代函数
k=1;
xk=x1;
k=k+1;
xk=h(xk); 
T=[1,x1,f(x1);k,xk,f(xk)];   %T初始值为2行3列，随着迭代增加行数
while abs(h(xk)-xk)>eps/2||f(xk)>eps/2
    if (df(xk)==0)||(k>5000)        %或
        error('迭代出错!或许是迭代初值的问题')
    else
        k=k+1;
        xk=vpa(h(xk),30); 
        T=[T;k,xk,f(xk)];
    end
end
fprintf('     k           xk              f(xk)\n');
disp(vpa(T,10))
fprintf('经过%d次迭代，函数方程根的近似解为：x=%.8f\n',k-1,T(k,2))%注意迭代次数为k-1
fprintf('f(x(%d))=%.8f\n',k-1,T(k,3))

懒了懒了，不想写了，再贴一个代码就去打游戏去；
简化牛顿法代码：

clc;clear;
format long;
syms x;
f(x)=input('请输入函数表达式,变量用x表示：f(x)=');
x1=input('请选择一个初始值:x1=');
eps=input('请输入停止精度要求:eps='); 
df(x)=diff(f,x);%f(x)的导数
g(x)=f(x)/df(x);
h(x)=x-g(x);%迭代函数
k=1;
xk=x1;
k=k+1;
xk=h(xk); 
T=[1,x1,f(x1);k,xk,f(xk)];   %T初始值为2行3列，随着迭代增加行数
while abs(h(xk)-xk)>eps/2||f(xk)>eps/2
    if (df(xk)==0)||(k>5000)        %或
        error('迭代出错!或许是迭代初值的问题')
    else
        k=k+1;
        xk=vpa(h(xk),30); 
        T=[T;k,xk,f(xk)];
    end
end
fprintf('     k           xk              f(xk)\n');
disp(vpa(T,10))
fprintf('经过%d次迭代，函数方程根的近似解为：x=%.8f\n',k-1,T(k,2))%注意迭代次数为k-1次
fprintf('f(x(%d))=%.8f\n',k-1,T(k,3))

溜了溜了~~~