para.m

抛物线法求解 方程的构造方法:给出[0,1]区间上的随机数(服从均匀分布)作为方程的根p*. 设你的班级数为a3,学号的后两位数分别为a2与a1,从而得到你的三次方程 例如:你的31班的1

syms ii 
for ii=1:20;
clear clear x y s e Y X w c p t;
r=rand(1,1)
syms x y z 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;


s(1)=0,s(2)=0.5,s(3)=1;

for n=3:20;
    s(n)=vpa(s(n),16)
    
    x=s(n);
    y=s(n-1);
    z=s(n-2);
    f1=a3*x^3+10*(a2+1)*x^2+a1*x+a0;
    f2=a3*y^3+10*(a2+1)*y^2+a1*y+a0;
    f3=a3*z^3+10*(a2+1)*z^2+a1*z+a0;
    f4=(f2-f1)/(y-x)+(((f3-f1)/(z-x)-(f2-f1)/(y-x))/(z-x))*(x-y);
    if (subs(f4)>=0)
        parab=x-2*f1/(f4+sqrt(f4*f4-4*f1*((f3-f1)/(z-x)-(f2-f1)/(y-x))/(z-x)));
    else
        parab=x-2*f1/(f4-sqrt(f4*f4-4*f1*((f3-f1)/(z-x)-(f2-f1)/(y-x))/(z-x)));
    end
    s(n+1)=subs(parab,{x,y,z},{s(n),s(n-1),s(n-2)});
    a=vpa((abs(s(n+1)-s(n))),16);
    e(n-2)=vpa((s(n-2)-r),16);
    e(n-1)=vpa((s(n-1)-r),16);
    
    if (subs(abs(e(n-1)))<10^(-16))
      break
    end;
    
    X(n-2)=log(abs(e(n-2)));
    Y(n-2)=log(abs(e(n-1)));
    
    w(n-2)=(n-2)^2;
    v00(n-2)=w(n-2);
    v10(n-2)=w(n-2)*X(n-2);
    v11(n-2)=w(n-2)*(X(n-2)^2);
    v02(n-2)=w(n-2)*Y(n-2);
    v12(n-2)=w(n-2)*X(n-2)*Y(n-2);
    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