syms x;
f1=(@(x)x^2-3*x+2-exp(x));
f2=@(x)x^3+2*x^2+10*x-20;
df1=matlabFunction(diff(f1,x));
df2=matlabFunction(diff(f2,x));
DieDai1=@(x)(x^2+2-exp(x))/3;
DieDai2=@(x)20/(x^2+2*x+10);

%

BuDongDieDai(DieDai1,1.5);
BuDongDieDai(DieDai2,1.5);

%

[x_star,inedx,it] = steffensen(DieDai1,1.5,1e-8);
[x_star,inedx,it] = steffensen(DieDai2,1.5,1e-8);

%

r=newton(f1,df1,1.5);
r=newton(f2,df2,1.5);

%

steffensen(DieDai2,50,1e-8)
newton(f2,df2,50);

%
% 扩展实验:比较不同算法的迭代次数优劣
%

x0=-50:1:50;
n=length(x0);
K_BuDongDieDai=zeros(n,1);
K_steffensen=zeros(n,1);
K_newton=zeros(n,1);
for i=1:n
    K_BuDongDieDai(i)=BuDongDieDai(DieDai1,x0(i));
    K_steffensen(i)=steffensen(DieDai1,x0(i),1e-8);
    K_newton(i)=newton(f1,df1,x0(i));
end
plot(x0,K_BuDongDieDai,x0,K_steffensen,x0,K_newton)
xlabel('初值')
ylabel('迭代次数,上限500次')
legend('不动点法','斯特芬森加速迭代','牛顿迭代')

%

x0=-200:1:200;
n=length(x0);
K_BuDongDieDai=zeros(n,1);
K_steffensen=zeros(n,1);
K_newton=zeros(n,1);
for i=1:n
    K_BuDongDieDai(i)=BuDongDieDai(DieDai2,x0(i));
    K_steffensen(i)=steffensen(DieDai2,x0(i),1e-8);
    K_newton(i)=newton(f2,df2,x0(i));
end

%

plot(x0,K_BuDongDieDai,x0,K_steffensen,x0,K_newton)
xlabel('初值')
ylabel('迭代次数,上限500次')
legend('不动点法','斯特芬森加速迭代','牛顿迭代')

%

DieDai3=@(x)(exp(x)-2)/(x-3)
x0=-200:1:200;
n=length(x0);
K_BuDongDieDai=zeros(n,1);
K_steffensen=zeros(n,1);
K_newton=zeros(n,1);
for i=1:n
    K_BuDongDieDai(i)=BuDongDieDai(DieDai3,x0(i));
    K_steffensen(i)=steffensen(DieDai3,x0(i),1e-8);
    K_newton(i)=newton(f1,df1,x0(i));
end
plot(x0,K_BuDongDieDai,x0,K_steffensen,x0,K_newton)
xlabel('初值')
ylabel('迭代次数,上限500次')
legend('不动点法','斯特芬森加速迭代','牛顿迭代')

%

x0=-50:0.5:0;
n=length(x0);
K_BuDongDieDai=zeros(n,1);
K_steffensen=zeros(n,1);
K_newton=zeros(n,1);
for i=1:n
    K_BuDongDieDai(i)=BuDongDieDai(DieDai3,x0(i));
    K_steffensen(i)=steffensen(DieDai3,x0(i),1e-8);
    K_newton(i)=newton(f1,df1,x0(i));
end
plot(x0,K_BuDongDieDai,x0,K_steffensen,x0,K_newton)
xlabel('初值')
ylabel('迭代次数,上限500次')
legend('不动点法','斯特芬森加速迭代','牛顿迭代')

%

BuDongDieDai(DieDai1,50)
steffensen(DieDai1,100,1e-8)
newton(f1,df1,100)

%

steffensen(DieDai3,-200,1e-8)

%

newton(f1,df1,-200)

%
%

function k=BuDongDieDai(DieDaiFun,x0)
r=inf;
k=0;%迭代次数
while abs(r-x0)>1e-8
     r=x0;
     y=DieDaiFun(x0);
     x0=y;
     %fprintf('\n次数k=%d;方程的根x=%.11f;f(x)=%.8f;误差为%f\n',k,r,y,abs(r-x0));%保留11位小数
     k=k+1;
     if abs(y)>=1e29
         k=500;
         break;%防止数值过大
     end
end 
end

function k=newton(f,df,x0)
r=inf;
k=0;%迭代次数
while abs(r-x0)>1e-8
     r=x0;
     y=f(x0);
     x0=x0-f(x0)/df(x0);
     %fprintf('\n次数k=%d;方程的根x=%.11f;f(x)=%.8f;误差为%f\n',k,r,y,abs(r-x0));%保留11位小数
     k=k+1;
     if abs(y)>=1e29
         k=500;
         break;%防止数值过大
     end
end 
end


function k=steffensen(phi,x,ep)
%斯特芬森加速迭代方法
% x为初始点
% x_star为当迭代成功时，输出方程的根
%   当迭代失败时，输出最后的迭代值;
% index为指标变量，当index=1时，表明迭代成功
% it为迭代次数
index=0;k=1;
%fprintf('\n次数k=0;方程的根x=%.11f;f(x)=%.11f\n',x,phi(x));
while k< 500
    x1=x;
    y=feval(phi,x);
    z=feval(phi,y);
    x=x-(y-x)^2/(z-2*y+x);

    %fprintf('\n次数k=%d;方程的根x=%.11f;f(x)=%.8f;误差为%f\n',k,x,y,x-x1);%保留11位小数
    if max(abs(y),abs(z))>=1e9
         k=500;
    break;%防止数值过大
    end
    if abs(x-x1)<ep
        index=1;break;
    end
    k=k+1;
end
x_star=x;it=k;
end