hypot.m

来自「抛物线法求解 方程的构造方法:给出[0,1]区间上的随机数(服从均匀分布)作为」· M 代码 · 共 51 行

syms ii 
for ii=1:20;
clear x y s e Y X w c p t;
r=vpa(rand(1,1),16)
syms x y 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=a3*y^3+10*(a2+1)*y^2+a1*y+a0;
hypoth=x-f1*(x-y)/(f1-f2);
s(1)=0,s(2)=1 ;
for n=2:20;
    s(n)=vpa(s(n),16)
    s(n+1)=subs(hypoth,{x,y},{s(n),s(n-1)});
    a=vpa((abs(s(n+1)-s(n))),16)
    e(n-1)=vpa((s(n)-r),16);
    e(n)=vpa((s(n+1)-r),16);
 
    if (subs(abs(e(n)))<10^(-16))
      break
    end;
    
    X(n-1)=log(abs(e(n-1)));
    Y(n-1)=log(abs(e(n)));
    
    w(n-1)=(n-1)^2;
    v00(n-1)=w(n-1);
    v10(n-1)=w(n-1)*X(n-1);
    v11(n-1)=w(n-1)*(X(n-1)^2);
    v02(n-1)=w(n-1)*Y(n-1);
    v12(n-1)=w(n-1)*X(n-1)*Y(n-1);
   
    if (subs(a)<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