程序为:
a=0;
b=pi/2;
a1=a;
a3=b;
a2=a;
a4=b;
while a3-a1>0.01
f1=-atan(2*sin(2*a1)/(5-cos(2*a1)));
f3=-atan(2*sin(2*a3)/(5-cos(2*a3)));
a2=(a1+a3)/2;
f2=-atan(2*sin(2*a2)/(5-cos(2*a2)));
c1=(f3-f1)/(a3-a1);
c2=((f2-f1)/(a2-a1)-c1)/(a2-a3);
if c2~=0
a4=0.5*(a1+a3-c1/c2);
f4=-atan(2*sin(2*a4)/(5-cos(2*a4)));
if (4-a1)*(a3-a4)>0
if a2
if f2
a3=a4;
else
a1=a2;
end
else
if f2
a1=a4;
else
a3=a2;
end
end
end
end
end
if f2
x=a2;
else
x=a4;
end
f=-atan(2*sin(2*x)/(5-cos(2*x)));
x=x*180/pi %最终结果
f=f*180/pi