advancenewton.m

来自「改进的牛顿法求解： 方程的构造方法:给出[0,1]区间上的随机数(服从均匀分布」· M 代码 · 共 56 行

syms ii 
for ii=1:20;
clear x s e Y X w c p t;
r=vpa(rand(1,1),16)
syms x s e Y X w c p t;
a3=10;a2=1;a1=9;a0=-a3*r^3-10*(a2+1)*r^2-a1*r;
f1=a3*x^3+10*(a2+1)*x^2+a1*x+a0;
f2=diff(f1,'x');
f3=diff(f2,'x');
newton=x-(f1*f2)/(f2^2-f1*f3);
s(1)=0.5;
for n=1:20;
    s(n)=vpa(s(n),16)
    e(n)=vpa(e(n),16);
    s(n+1)=subs(newton,x,s(n));
    a=vpa((abs(s(n+1)-s(n))),16);
    e(n)=vpa((s(n)-r),16);
    e(n+1)=vpa((s(n+1)-r),16);
 
    if (subs(abs(e(n+1)))==0)
      break
    end;
    
    X(n)=log(abs(e(n)));
    Y(n)=log(abs(e(n+1)));
    
    w(n)=n^2;
    v00(n)=w(n);
    v10(n)=w(n)*X(n);
    v11(n)=w(n)*(X(n)^2);
    v02(n)=w(n)*Y(n);
    v12(n)=w(n)*X(n)*Y(n);
   
    if (subs(a)<5*10^(-16))
        
        break
    end;
     
    
end
    
    d00=sum(v00);
    d10=sum(v10);
    d11=sum(v11);
    d02=sum(v02);
    d12=sum(v12);
    q1=d00*c+d10*p-d02;
    q2=d10*c+d11*p-d12;
    S=solve(q1,q2,c,p);
    S.c;
    S.p;
    t=S.p;
T(ii)=t

end
sum(T)/20