- 积分
- 5
- 注册时间
- 2003-12-30
- 仿真币
-
- 最后登录
- 1970-1-1
|
zz from simulate of freecity
two versions
I programed the first and blinderbrain did the second
最后输出的x去掉了重复的直,x0是有解的初值和最后的解的对照
iter迭代次数
由于计算同时删除了一些不可能的节点
时效上会快些
如果能够较快判断出一些明显不收敛的点而不是运行maxit次迭代
可能效率会更好
以后有空再看看
function [x,x0,iter]=mynewton(f,x0,maxit,tol);
df=diff(f);
if nargin<4
tol=10^(-12);
if nargin<3
maxit=10000;
if nargin<2
x0=0;
end
end
end
x0=x0(:);
px=x0;
px0=[];
x=[];
for i=1:maxit
fx=subs(f,px);
dfx=subs(df,px);
x=[x;px(abs(fx)<tol)];
px0=[px0;x0(abs(fx)<tol)];
id=all([abs(dfx)>10*eps abs(fx)>tol],2);
px=px(id);
if isempty(px)
break;
end
x0=x0(id);
px=px-fx(id)./dfx(id);
end
x0=[px0 x];
x=sort(x);
x=x(abs(diff([-inf;x]))>tol);
iter=i;
下面是我前段时间帮人家编写的用newton法求解一维多项式根的函数,
优点是初始点可以向量化,试过超过10000个初始点的情况,
运算时间不超过一分钟,然后初始点还可以是复平面里面的点。
大家参考一下,如果有什么不对的地方多多指教。
由于仓促,所以没有写太多的解释性的东西
%-------------begin-------------------
function [x,x0,fval,exitflag]=newton(func,x0,tol,maxiter)
% func为函数的字符串格式,如果输入x0为向量则要求写成向量运算的格式
% x0为迭代的初始点,可以是向量的形式
% tol为误差容限
% maxiter为最大允许迭代次数
dfunc=diff(func);
x0=x0(:);
p00=x0;
p0=x0;
if nargin<4
maxiter=20000;
if nargin<3
tol=10^(-12);
if nargin<2
disp('Error:Not Enough Input Parameters!');
end
end
end
ns=length(x0);
Data=[];
for k=1:maxiter
if any(abs(subs(dfunc,p0))<10*eps)
idx1=find(abs(subs(dfunc,p0))<=10*eps);
Data=[Data;[p00(idx1),p0(idx1),subs(func,p0(idx1)),...
-5*ones(length(idx),1)]];
p0(idx1)=[];
p00(idx1)=[];
else p=p0-subs(func,p0)./subs(dfunc,p0);
end
if length(p00)==0
[x,x0,fval,exitflag]=lastperform(Data);
return;
end
if any(abs(p-p0)<tol)
idx2=find(abs(p-p0)<tol);
Data=[Data;[p00(idx2),p(idx2),subs(func,p(idx2)),...
1*ones(length(idx2),1)]];
p(idx2)=[];
p00(idx2)=[];
end
if length(p00)==0
[x,x0,fval,exitflag]=lastperform(Data);
return;
end
p0=p;
end
Data=[Data;[p00,p,subs(func,p),-1*ones(length(p),1)]];
[x,x0,fval,exitflag]=lastperform(Data);
% Private Function
function [x,x0,fval,exitflag]=lastperform(Data)
Data=sortrows(Data,1);
x=Data(:,2);fval=Data(:,3);exitflag=Data(:,4);
x0=Data(:,1);
%--------------end------------------
[ 本帖最后由 bainhome 于 2006-10-15 12:06 编辑 ] |
评分
-
1
查看全部评分
-
|