% 用来求π--割圆术--周长/直径【大体完成,但估计是由于精度问题,造成后续的数值出错了,ppib(i)-ppi(i)值小于0】
clc;clear all;close all
format long
ra = 10; % 圆半径
bianN = 30; % 分割次数bianN,边数为rn,内接边长为rr,内接周长为l
rr = zeros(1,bianN); rr(1) = ra; % 边长
rn = zeros(1,bianN); rn(1) = 6; % 边数
l = zeros(1,bianN); l(1) =6*rr(1); % 周长
ppi= zeros(1,bianN); ppi(1) =l(1)/(2*ra);% 周长/直径
rrb= zeros(1,bianN); rrb(1) = 2/sqrt(3)*ra; % 边长
lb = zeros(1,bianN); lb(1) =6*rrb(1); % 周长
ppib=zeros(1,bianN); ppib(1) =lb(1)/(2*ra);% 周长/直径
jingdu = zeros(1,bianN);jingdu(1) = ppib(1)-ppi(1);% 计算精度
fprintf('当前边数为:%d 圆周率为:%.17f与%.17f之间,精度为:%.17f\n',rn(1),ppi(1),ppib(1),jingdu(1))
for i = 2:bianN
% 内接正rn边形
rn(i) = 3*2^i;
rr(i) = sqrt((rr(i-1)/2)^2+(ra-sqrt(ra^2-(rr(i-1)/2)^2))^2);
l(i) = rn(i)*rr(i);
ppi(i) = l(i)/(2*ra);
% 外接rn边形
rrb(i) = (rrb(i-1)^2 - 4*(sqrt(ra^2+(rrb(i-1)/2)^2)-ra)^2)/(2*rrb(i-1));
lb(i) = rn(i)*rrb(i);
ppib(i) = lb(i)/(2*ra);
jingdu(i) = ppib(i) - ppi(i);
if jingdu(i)<0
break;
end
% % 经查,1.公式是错的,2.此方法中因为存在rb = rbi(1)的固定关系,所以无法计算,
% % 【但是,对于这种逆向思想还是灵光一闪,值得收藏的】
% % % 外接正rn边形--因为r(i)表达式不好得到,而r(i-1)的表示式(用r(i)表示)则好找,所以此处为反向求解的思维
% % rb = 1; %
% % rbi(i) = rb; % 边长随便设置
% % for j = i-1:-1:1
% % rbi(j) = 2*sqrt((rbi(j+1)/2)^2+(sqrt(rbi(j+1)/2+rb)-rb)^2)+rbi(j+1);
% % end
% % rB2A = ra/rbi(1); % 外接的放大倍率,需要用来转换数据会正常比例
% % lb0 = rn(i)*rbi(i); % 外接数据下的周长,非内接比例
% % lb(i) = lb0/rB2A; % 恢复正常比例--变为与内接等比例
% % ppib(i) = lb(i)/(2*ra);
% % jingdu(i) = abs(ppib(i)-ppi(i));
fprintf('当前边数为:%d 圆周率为:%.17f与%.17f之间 精度为:%.17f\n',rn(i),ppi(i),ppib(i),jingdu(i))
end