improvednewton2.m

给定一个三次方程,改进的Newton法求解,并研究各种算法的收敛阶数.

%当f(x)= (10*x.^3+80*x.^2+6*x+a0)*(x-p*)时
[improvednewton2.m]
function sa = improvednewton2(x0,epsilon,max)
% IMPROVEDNEWTON calculate the improvednewton method approximation
% x0 is the initial value 
% epsilon is the tolerance for x0 
% max is the maximum number of iterations
syms x;
syms pstar;
q=1000;
pstar=rand(1,q);% p*
step=zeros(1,q); %迭代步数
order=zeros(3,q); %收敛阶数
g=inline('-10*x^3-80*x^2-6*x');
df=inline('1-(10*x^3+80*x^2+6*x+a0)*(x-pstar)*((60*x+160)*(x-pstar)+2*(30*x.^2+160*x+6))/((30*x.^2+160*x+6)*(x-pstar)+(10*x.^3+80*x.^2+6*x+a0))^2');
U=zeros(1,max);
U(1)=x0;
for i=1:q
    a0=feval(g,pstar(i));    
    f=inline('(10*x.^3+80*x.^2+6*x+a0)*(x-pstar)/((30*x.^2+160*x+6)*(x-pstar)+(10*x.^3+80*x.^2+6*x+a0))');
    for j=1:max
        x1=x0-feval(f,a0,pstar(i),x0)/feval(df,a0,pstar(i),x0);
        U(j+1)=x1;
        err=abs(x1-x0);
        x0=x1;
       if err<epsilon
           break;
       end
    end
    for k=1:j
        Us(k)=log(abs(U(k)-pstar(i)));
    end
    X=Us(1:j-1);
    Y=Us(2:j); 
    if  j==max    %不收敛
        step(i)=0;
        order(:,i)=0;
        q=q-1; 
    else if j==2 %收敛
            X=Us(1:j);
            Y=Us(2:j+1);            
    else          %收敛
        step(i)=j; 
        xmean=mean(X);
        ymean=mean(Y);
        Gs=X-xmean;
        G=zeros(1,j-1);
        for h=1:3
             for d=1:j-1
                 G(d)=Gs(d)*(d^(h-1));
             end
            sumx2=G*(X-xmean)';
            sumxy=(Y-ymean)*G';
            order(h,i)=sumxy/sumx2; 
        end
        end
    end    
end 
q
v=zeros(1,3);
v=sum((order)'); 
averagesteps=sum(step)/q
averageorder1=v(1)/q
averageorder2=v(2)/q
averageorder3=v(3)/q