计算方法单点截弦matlab实现

function Dan(fun,a,b,e)
% 使用单点截弦法求方程 fun=0 在区间[x1,x2]处的根 e为误差限，默认为 eps
% 对于任意的 x 属于 [a,b] ， diff(fun,x) 与 diff(fun,x,2) 连续且不变号
% 选择x0为不动点，xn为动点， 要求 fun(x_0)*diff(fun,x_0,2) > 0
% 其中迭代公式为 y=a-(x-a)/(fun(x)-fun(a))*fun(a)

% 作者：野渡无人

% 最后修改日期：2008.3.31

% Dan('x^3-2*x-5',2,3,10^(-5))
% n= 0 x0= 3.00000000000000000
% n= 1 x1= 2.00000000000000000
% n= 2 x2= 2.05882352941176450
% n= 3 x3= 2.08126365984502290
% n= 4 x4= 2.08963921009084740
% n= 5 x5= 2.09273957431800550
% n= 6 x6= 2.09388370846184820
% n= 7 x7= 2.09430545112525970
% n= 8 x8= 2.09446084576648770
% n= 9 x9= 2.09451809338575720
% n= 10 x10= 2.09453918229401690

if nargin==3
e=eps;
end

if a>=b
error('a>=b！请重新输入a,b值！');
end

if nargin<3
error('请正确输入至少三个参数！');
end

syms x y
y=fun; %定义符号函数，便于判断单点截弦法时隔根区间的要求以及初值的选定
y_diff=diff(y,x,1) ; % 函数的一阶导
y_diff2=diff(y,x,2) ;% 函数的二阶导

% 判断单点截弦法时隔根区间的要求，即：对于x属于[a,b],f(x)的一阶导以及二阶导应连续且不变号
for i=a:(b-a)/100:b
if sign(subs(y_diff,x,i))~=sign(subs(y_diff2,x,i))
error('对于x属于[a,b],f(x)的一阶导以及二阶导应连续且不变号！故所输入的参数有误！');
end
end

if sign(subs(y,x,a))==sign(subs(y_diff,x,a))
x0=a ; % a为初值x0
x1=b ; % b为x1
elseif sign(subs(y,x,b))==sign(subs(y_diff,x,b))
x0=b;
x1=a;
else
error('隔根区间端点不能作初值！');
end

format long;
n=0;
fprintf('n=%3d x%d= %12.17f /n',n,n,x0);
n=n+1;
fprintf('n=%3d x%d= %12.17f /n',n,n,x1);
xn=x0-(x1-x0)/(subs(y,x1)-subs(y,x0))*subs(y,x0);
while abs(xn-x1) > e
n=n+1;
fprintf('n=%3d x%d= %12.17f /n',n,n,xn);
x1=xn;
xn=x0-(x1-x0)/(subs(y,x1)-subs(y,x0))*subs(y,x0);
end