Matlab仿真两种方法求圆周率π

本文用到的matlab源码下载地址

1.1 蒲丰投针法求π

1.1.1 蒲丰投针算法介绍

1.1.2 混合同余算法介绍

1.1.3 程序及注释

混合同余法产生随机数——Mixrand()函数文件Mixrand.m

function [q] = Mixrand(n)%A是取余函数(Ax+n)的参变量a=1;A=4*a+1;%循环实验x=1;%当满足%x=rem(A*x+n,2^k)%A=4*a+1%y=x/2^k%T>2^m(优化后的混合同举法）for i = 1:nx=rem(A*x+i,2^20);d=x/(2^20);X(i)= d;A=4*a+1;%该值舍去(提高精度)%进行两次该算法的原因：因为发现y的离散性更好,故猜测多次迭代可提高精度x=rem(A*x+i,2^20);     d=x/(2^20);X(i)= d;a=a+1;x=rem(A*x+i,2^20);d=x/2^20;Y(i)= d;a=a+1;end%第二次混合同余法,取得更精准的数据求平均for i = 1:nx=rem(A*x+i,2^20);d=x/(2^20);X(i)= d;A=4*a+1;%该值舍去(提高精度)%进行两次该算法的原因：因为发现y的离散性更好,故猜测多次迭代可提高精度x=rem(A*x+i,2^20);     d=x/(2^20);X(i)= d;a=a+1;x=rem(A*x+i,2^20);d=x/2^20;Y(i)= d;a=a+1;endq=Y;
end

模拟投针过程文件.test6.m

clear all
d=1;%两条边的距离
l=0.6;%针的长度%总的实验次数
n = input('请输入n:');
N=n;
m=0;
mpi=0;
x=Mixrand(3*N);
for i=1:n
x1(i)=x(i);
y1(i)=x(i+N);
Mixpi(i)= pi*x(i+2*N);%mpi为接收到的(0-2pi/)的随机数
if x1(i)<l*sin(Mixpi(i))/2m=m+1;
end
end
rectangle('Position',[-0.6 -0.6 2.2 2.2],'EdgeColor','r','LineWidth',3);%画矩形
hold on;
rectangle('Position',[0 0 1 1],'EdgeColor','m','LineWidth',7);
axis([-0.7 1.7 -0.7 1.7]);
for i=1:nP(i)=cos(2*Mixpi(i));Q(i)=sin(2*Mixpi(i));x11(i)=x1(i)+l*P(i);y11(i)=y1(i)+l*Q(i);
hold on;
end
plot(x1,y1,'X');
title('模拟蒲丰投针');
xlabel('X轴');
ylabel('Y轴');
X=[x1;x11];
Y=[y1;y11]
line(X,Y);
q=l*n/(d*m);

主函数文件tz.m

d=1;%两条边的距离%总的实验次数
N=n;
m=0;
mpi=0;
x=Mixrand(3*n);
for i=1:n
x1(i)=x(i);
y1(i)=x(i+N);
Mixpi(i)= pi*x(i+2*N);%mpi为接收到的(0-2pi/)的随机数
if x1(i)<l*sin(Mixpi(i))/2m=m+1;
end
end
q=l*n/(d*m);

变量分析文件anlgnse.m

clear all;
t=10^4:10^4:1*10^6;
M=length(t);
for i=1:length(t)       %采集Pi的个数m(i)=tz(t(i),0.6);       %1000-10000中选择
endplot(m,'X');hold on;for i=1:Mx(i)=i;y(i)=3.1415926;endplot(x,y,'g');
title('Pi随着N变化的变化图');
xlabel('N点');
ylabel('Pi值');

分析l/a（针长/边长）对实验结果的影响文件anlgnse1.m

clear all;
%t=10^4:10^4:1*10^6;
I=0.01:0.01:1
M=length(I);
for i=1:length(I)       %采集Pi的个数m(i)=tz(100000,I(i));       %1000-10000中选择
end
for i=1:Mx(i)=i;y(i)=pi;endplot(I,m,'X');hold on;axis([0 1 3 3.3]);
title('Pi随着l/a变化的变化图');
xlabel('l/a');
ylabel('Pi值');

1.1.4 仿真结果

1.2 蒙特卡洛法求π

1.2.1 蒙特卡洛算法介绍

1.2.2 关键程序

算法2——蒙特卡洛方法：

%每次重置所有变量clear all;m = 0;n=input('请输入n:');%使用的圆的半径r = 1;%A是取余函数(Ax+n)的参变量a=1;A=4*a+1;%循环实验x=1;d1=Mixrand(2*n);
N=n;for i = 1:nX(i)= d1(i);Y(i)= d1(i+N);if (X(i)^2 + Y(i)^2 <= r^2)X1(i)=X(i);Y1(i)=Y(i);m = m + 1;elseX2(i)=X(i);Y2(i)=Y(i);endendp1=4 * m / n;%显示结果X0=0;Y0=0;t=0:pi/40:2*pi;r=1;XX=abs(X0+r*cos(t));YY=Y0+r*sin(t);plot(X1,Y1,'.',X2,Y2,'.',XX,YY,'-');
fprintf('当总实验次数n = %d时计算出来的圆周率:Pi = %d\n',n, p1);

anlgnse函数

clear all;
t=5*10^4:10^4:1*10^5;
M=length(t);
for i=1:length(t)       %采集Pi的个数m(i)=mt(t(i));       %1000-10000中选择
endplot(m,'X');hold on;for i=1:Mx(i)=i;y(i)=3.1415926;endplot(x,y,'g');
title('Pi随着N变化的变化图');
xlabel('N点');
ylabel('Pi值');

1.2.3 仿真结果

当投针N变大蒙特卡洛演示结果如下：

Matlab仿真两种方法求圆周率π相关推荐

牛客 Tree(最小深度总和)（两种方法求重心）难度⭐⭐⭐
题目链接牛妹有一张连通图,由n个点和n-1条边构成,也就是说这是一棵树,牛妹可以任意选择一个点为根,根的深度deprootdep_{root}deproot为0,对于任意一个非根的点,我们将他到 ...
centos 7安装matlab的两种方法（桌面安装和命令行安装）
matlab安装说明安装之前一直以为命令行安装(静默安装)完就是命令行界面,安装成功后才发现还是有桌面的,还跟桌面安装的一模一样.所以,个人建议对linux不太熟悉的还是用桌面版安装,虽然会有点卡顿 ...
来手把手教你通过Matlab用两种方法实现图像压缩与解压（附超详细代码），赶紧点赞收藏吧
图像压缩方法 DCT图像压缩 DCT原理介绍 DCT和它解压时的反运算的具体算法详细实现代码结果展示行程编码压缩与解压读入图像图像转为矩阵行程编码压缩行程编码解压显示图像完整代码附录 ...
求逆元的两种方法+求逆元的O(n)递推算法
到国庆假期都是复习阶段..所以把一些东西整理重温一下. gcd(a,p)=1,ax≡1(%p),则x为a的逆元.注意前提:gcd(a,p)=1; 方法一:拓展欧几里得 gcd(a,p)=1,ax≡1( ...
C语言两种方法求圆的面积与周长编程
方法一:程序如下: #define _CRT_SECURE_NO_WARNINGS #include<stdio.h> int main() {const float pi = 3.14; ...
两种方法求最大公约数和最小公倍数
#include<iostream> using namespace std;// 最小公倍数 = 两数相乘/最大公约数 /* //辗转相除法 int gcd(int a, int b){ ...
matlab中三种方法计算 Ax b,在MATLAB中，方程Ax=B的解可以用哪个命令求得？ matlab 求助解方程组...
matlab中解方程组还是很方便的,例如,对于代数方程组Ax=b(A为系数矩阵,非奇异)的求解,MATLAB中有两种方法: (1)x=inv(A)*b - 采用求逆运算解方程组: (2)x=A\B - ...
matlab 求矩阵秩,求矩阵秩的两种方法及MATLAB的应用
摘要: 高等代数是一门逻辑思维比较强和理论知识比较深的学科, 它具有丰富的数学知识, 涉及许多重要的数学思想, 其在数学领域的应用很广泛, 如行列式.矩阵的相关计算和求解线性方程组的解方面的应 ...
用matlab计算稳态误差,利用Matlab求稳态误差的两种方法.
利用Matlab求稳态误差的两种方法摘要:稳态误差是系统控制精度或抗扰动能力的一种度量,它是稳态性能的一个重要指标.本文介绍利用Matlab的控制系统工具箱和Simulink工具箱求取系统误差稳态的 ...
matlab中给信号添加高斯白噪声的两种方法，awgn计算过程,randn函数
y=awgn(x,snr,px_dBW) 给信号x添加噪声功率为某个值的高斯白噪声. snr为信噪比,单位dB. px_dBW为信号x的指定功率(注意,是指定功率,而不是x本身的功率),单位dBW. ...

Matlab仿真两种方法求圆周率π

目录:

1.1 蒲丰投针法求π

1.1.1 蒲丰投针算法介绍

1.1.2 混合同余算法介绍

1.1.3 程序及注释

1.1.4 仿真结果

1.2 蒙特卡洛法求π

1.2.1 蒙特卡洛算法介绍

1.2.2 关键程序

1.2.3 仿真结果

Matlab仿真两种方法求圆周率π相关推荐

最新文章

热门文章