%當(dāng)f(x)= 10*x.^3+80*x.^2+6*x+a0時(shí)[newton1.m]function sa = newton1(x0,epsilon,max)% NEWTON calculate the newton method approximation% x0 is the initial value % epsilon is the tolerance for x0 % max is the maximum number of iterationssyms x;syms pstar;syms q;syms a0;q=1000;pstar=rand(1,1000); % p*step=zeros(1,1000); %迭代步數(shù)order=zeros(3,1000); %收斂階數(shù)g=inline('-10*x.^3-80*x.^2-6*x');df=inline('30*x^2+160*x+6');U=zeros(1,max);U(1)=x0;for i=1:1000    a0=feval(g,pstar(i));        f=inline('10*x.^3+80*x.^2+6*x+a0 ');    for j=1:max        x1=x0-feval(f,a0,x0)/feval(df,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 qv=zeros(1,3);v=sum((order)'); averagesteps=sum(step)/qaverageorder1=v(1)/qaverageorder2=v(2)/qaverageorder3=v(3)/q