MATLAB中的圆周率计算新方法

直接复刻一下，matlab代码如下：

先测试一下公式：

an=@(n)2*pi*(log2(n+3)-log2(n+2));
pai=0;
for i=1:3
pai=pai+sin(an(i))
end
m=2;
 ans=0;
while ans<=3.1415926
pai=pai+sin(an(2*m))+sin(an(2*m+1))-sin(an(m-1));
ans=0.5*pai;
m=m+1;
end
m
ans

发现：当n=20840的时候才能算出pi=3.1415926，所以说这是个笨办法！

详细代码如下：

%Author@https://space.bilibili.com/51126650
%% canvas
format long
fig=figure();
screenSize=get(0,'ScreenSize');
fig.Position=[screenSize(1,[3,4])./10,screenSize(4).*0.6,screenSize(4).*0.6];
fig.Name='Pi';
fig.NumberTitle='off';
ax=axes(fig);
ax.Position=[0 0 1 1];
hold(ax,'on');
ax.XLim=[-1.5 1.5];
ax.YLim=[-1.5 1.5];
ax.XTick=[];
ax.YTick=[];
ax.XColor='none';
ax.YColor='none';
ax.Color=[249,247,246]./255;
%% Core
%画个圆
x = 0;
y = 0;
r = 1;
rectangle('Position',[x-r,y-r,2*r,2*r],'Curvature',[1,1],'EdgeColor','m',LineStyle='--')
an=@(n)2*pi*(log2(n+3)-log2(n+2));                %an公式
sum_an=0;                                                          %an的累加，初始值为0
sec_Hdl=plot([0,0],[0,1],'r--','LineWidth',1.5);     %秒针样式设计
cha_Hdl=plot([0,0],[0,1],'k--','LineWidth',1.5);    %针头痕迹设计
plot([0,0],[0,1],'k','LineWidth',1.5);                      %绘制第一个三角形的第一个边
x=0;
y=1;                                                                    %起始针头位置
pai=0;                                                                 %派的累加，初始值为0
%绘制可变化的文本框及其样式
ax1=annotation('textbox',FontSize=20,Position=[0.4 0.15 0 0],String=['\pi≈',num2str(0)]);
%前三个三角形
dtn=30;     %数值越大扫描转速越慢
for j=1:3
dt=linspace(sum_an,sum_an+an(j),dtn);
for i=1:dtn
    sec_t=pi/2-dt(i);
    sec_Hdl.XData=[0,cos(sec_t)];
    sec_Hdl.YData=[0,sin(sec_t)];            %通过改变XYdata属性实现秒针动态变化
    cha_Hdl.XData=[x,cos(sec_t)];
    cha_Hdl.YData=[y,sin(sec_t)];
    if sec_t <=(sum_an+an(j))
       drawnow
       pause(0.02)
    end
end
plot([0,cos(sec_t)],[0,sin(sec_t)],'k','LineWidth',2);%绘制第j个三角形的第某个边
maderand=[rand(),rand(),rand()];                          %生成随机数用于随机颜色填充
fill([0 x  cos(sec_t) ],[0 y sin(sec_t)],maderand)      %填充三角形以随机颜色
x=cos(sec_t);
y=sin(sec_t);                                                           %特殊针头位置记录
sum_an=sum_an+an(j);                                          %an的累加
pai=pai+sin(an(j));                                                 %派的计算
ax1.String=['\pi≈',num2str(0.5*pai,'%4.7f')];          %修改string属性以实现动态文字显示
end
%后边无穷多三角形
cha_x=0;
cha_y=1;    %针头位置初始化以纠正极度微小误差
x=0;
y=1;            %针头位置初始化以纠正极度微小误差
sum_an=0;   %an的累加
dtn=20;     %数值越大扫描转速越慢
for m=2:1000
    uistack(sec_Hdl,'top') %秒针的图层置顶显示
    dt=linspace(sum_an,sum_an+an(2*m)+an(2*m+1),dtn);
for n=1:dtn
    sec_t=pi/2-dt(n);
    sec_Hdl.XData=[0,cos(sec_t)];
    sec_Hdl.YData=[0,sin(sec_t)];
    cha_Hdl.XData=[cha_x,cos(sec_t)];
    cha_Hdl.YData=[cha_y,sin(sec_t)];
    if n==(dtn/2+1)  %修改dtn过小时注意此项
        xx=cos(pi/2-(sum_an+an(2*m)));
        yy=sin(pi/2-(sum_an+an(2*m)));
        plot([cha_x,xx],[cha_y,yy],'k','LineWidth',1.5);
        cha_x=xx;
        cha_y=yy;
        drawnow
        pause(0.01) %注释掉此项可以加快秒针转速
    elseif sec_t <=(sum_an+an(2*m)+an(2*m+1))
        drawnow
        pause(0.01)%注释掉此项可以加快秒针转速
    end
end
    plot([xx,cos(sec_t)],[yy,sin(sec_t)],'k','LineWidth',1.5);
    fill([xx x  cos(sec_t) ],[yy y sin(sec_t)],[rand(),rand(),rand()])
    x=cos(sec_t);
    y=sin(sec_t);
    cha_x=x;
    cha_y=y;
sum_an=sum_an+an(2*m)+an(2*m+1);
pai=pai+sin(an(2*m))+sin(an(2*m+1))-sin(an(m-1));
ax1.String=['\pi≈',num2str(0.5*pai,'%4.7f')];
end
delete(sec_Hdl); %不再显示秒针

代码结果如下GIF：