随机模拟算法求解圆周率

圆周率(π)这个东西是从小学开始一直陪伴我们的,这里使用使用蒙特卡洛算法来产生大量的随机数求解π的近似值。

计算方式

首先我们知道 正方形的面积公式是S1 = a * a,圆形的面积S2 = π * r * r;
所以以圆的直径为正方形边长,可以得出π的表达式。

π = 4 * S2 / S1 

这样一来,重点就是求解正方形和圆形的面积,这里使用在一正方形区域内圆内产生相应的随机点,
为了便于可视化分析,在圆内和圆外的点分别用不同的颜色来表示,最后圆内产生的点数就近似记作圆的面积,区域内的点数记作正方形的面积。

可视化分析-小数据集

这里用Java swing来求解。

首先看看Circle的model类

public class Circle {

    private int x, y, r;

    public Circle(int x, int y, int r){
        this.x = x;
        this.y = y;
        this.r = r;
    }

    public int getX(){ return x; }
    public int getY(){ return y; }
    public int getR(){ return r; }
  // 判断点是否包含在圆内
    public boolean contain(Point p){
        return Math.pow(p.x - x, 2) + Math.pow(p.y - y, 2) <= r*r;
    }
}

这里定义了圆的坐标,以及圆的半径。
然后看看点的model类

public class MonteCarloPiData {

    private Circle circle;
    private LinkedList<Point> points;
    private int insideCircle = 0;

    public MonteCarloPiData(Circle circle){
        this.circle = circle;
        points = new LinkedList<Point>();
    }

    public Circle getCircle(){
        return circle;
    }
    // 得到点的数量
    public int getPointsNumber(){
        return points.size();
    }
  
    public Point getPoint(int i){
        if(i < 0 || i >= points.size()) {
            throw new IllegalArgumentException("out of bound in getPoint!");
        }

        return points.get(i);
    }
  //添加点到链表中
    public void addPoint(Point p){
        points.add(p);
        // 计算圆内点的数量
        if(circle.contain(p)) {
            insideCircle++;
        }
    }
    // π值的计算 
    public double estimatePi(){
        if(points.size() == 0) {
            return 0.0;
        }
        int circleArea = insideCircle;
        int squareArea = points.size();
        return (double)circleArea * 4 / squareArea;
    }
}

接下来就是辅助类

public class AlgoVisHelper {
     
    public static final Color Red = new Color(0xF44336);
    public static final Color Green = new Color(0x4CAF50);

    private AlgoVisHelper() {
    
    // 绘制圆
    public static void strokeCircle(Graphics2D g, int x, int y, int r){

        Ellipse2D circle = new Ellipse2D.Double(x-r, y-r, 2*r, 2*r);
        g.draw(circle);
    }
   // 填充圆
    public static void fillCircle(Graphics2D g, int x, int y, int r){

        Ellipse2D circle = new Ellipse2D.Double(x-r, y-r, 2*r, 2*r);
        g.fill(circle);
    }
   // 绘制矩形
    public static void strokeRectangle(Graphics2D g, int x, int y, int w, int h){

        Rectangle2D rectangle = new Rectangle2D.Double(x, y, w, h);
        g.draw(rectangle);
    }
  // 填充矩形
    public static void fillRectangle(Graphics2D g, int x, int y, int w, int h){

        Rectangle2D rectangle = new Rectangle2D.Double(x, y, w, h);
        g.fill(rectangle);
    }
  // 设置颜色
    public static void setColor(Graphics2D g, Color color){
        g.setColor(color);
    }
   // 设置绘制宽度
    public static void setStrokeWidth(Graphics2D g, int w){
        int strokeWidth = w;
        g.setStroke(new BasicStroke(strokeWidth, BasicStroke.CAP_ROUND, BasicStroke.JOIN_ROUND));
    }
  // 延长动画时间
    public static void pause(int t) {
        try {
            Thread.sleep(t);
        }
        catch (InterruptedException e) {
            System.out.println("Error sleeping");
        }
    }
}

然后是渲染过程类

public class AlgoFrame extends JFrame{

   private int canvasWidth;
   private int canvasHeight;

   public AlgoFrame(String title, int canvasWidth, int canvasHeight){

       super(title);

       this.canvasWidth = canvasWidth;
       this.canvasHeight = canvasHeight;
        // 实例化画笔  刷新视图
       AlgoCanvas canvas = new AlgoCanvas();
       setContentPane(canvas);
       pack();

       setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
       setResizable(false);

       setVisible(true);
   }

   public AlgoFrame(String title){

       this(title, 1024, 768);
   }

   public int getCanvasWidth(){return canvasWidth;}
   public int getCanvasHeight(){return canvasHeight;}

   private MonteCarloPiData data;
   public void render(MonteCarloPiData data){
       this.data = data;
       repaint();
   }

   private class AlgoCanvas extends JPanel{

       public AlgoCanvas(){
           // 双缓存
           super(true);
       }

       @Override
       public void paintComponent(Graphics g) {
           super.paintComponent(g);

           Graphics2D g2d = (Graphics2D)g;

           // 抗锯齿
           RenderingHints hints = new RenderingHints(
                   RenderingHints.KEY_ANTIALIASING,
                   RenderingHints.VALUE_ANTIALIAS_ON);
           hints.put(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY);
           g2d.addRenderingHints(hints);

           // 具体绘制
           Circle circle = data.getCircle();
           AlgoVisHelper.setStrokeWidth(g2d, 3);
           AlgoVisHelper.strokeCircle(g2d, circle.getX(), circle.getY(), circle.getR());
           
           /**
             * 点在圆内 就绘制成红色
             * 在圆外 绘制成绿色
             */ 
           for(int i = 0 ; i < data.getPointsNumber() ; i ++){
               Point p = data.getPoint(i);
               if(circle.contain(p)) {
                   AlgoVisHelper.setColor(g2d, AlgoVisHelper.Red);
               } else {
                   AlgoVisHelper.setColor(g2d, AlgoVisHelper.Green);
               }
               // 将点填充到圆内
               AlgoVisHelper.fillCircle(g2d, p.x, p.y, 3);
           }
       }

       @Override
       public Dimension getPreferredSize(){
           return new Dimension(canvasWidth, canvasHeight);
       }
   }
}

最后就是视图类了

public class AlgoVisualizer {

    private static int DELAY = 40;
    private MonteCarloPiData data;
    private AlgoFrame frame;
    private int N;

    public AlgoVisualizer(int sceneWidth, int sceneHeight, int N){

        if(sceneWidth != sceneHeight) {
            throw new IllegalArgumentException("This demo must be run in a square window!");
        }
        this.N = N;

        Circle circle = new Circle(sceneWidth/2, sceneHeight/2, sceneWidth/2);
        data = new MonteCarloPiData(circle);

        // 初始化视图
        EventQueue.invokeLater(() -> {
            frame = new AlgoFrame("Get Pi with Monte Carlo", sceneWidth, sceneHeight);

            new Thread(() -> {
                run();
            }).start();
        });
    }
     // 动画主要逻辑
    public void run(){

        for(int i = 0 ; i < N ; i ++){
            // 每次绘制100个点
            if (i % 100 == 0) {
                frame.render(data);
                // 设置每次产生点的时间
                AlgoVisHelper.pause(DELAY);
                System.out.println(data.estimatePi());
            }
            int x = (int)(Math.random() * frame.getCanvasWidth());
            int y = (int)(Math.random() * frame.getCanvasHeight());
            data.addPoint(new Point(x,y));
        }
}
    public static void main(String[] args) {

        int sceneWidth = 800;
        int sceneHeight = 800;
        // 设置点的数量
        int N = 20000;

        AlgoVisualizer vis = new AlgoVisualizer(sceneWidth, sceneHeight, N);
    }
}

运行后,就可以看到点的绘制过程,最终结果如下


捕获.PNG

很显然 可以看到结果 精确到小数点后一位,点就无法绘制了,要想更加精确,试试更大的数据集。

无视图-大数据集

public class MonteCarloExeperiment {
    private int squareSide;
    private int N;
    private int outputInterval = 100;

    public MonteCarloExeperiment(int squareSide, int N) {
        if (squareSide <= 0 || N <= 0) {
            throw new IllegalArgumentException("squareSide and N must large than 0!");
        }
        this.squareSide = squareSide;
        this.N = N;
    }
    public void setOutputInterval(int interval) {
        if (interval <= 0) {
            throw new IllegalArgumentException("interval must be larger than 0!");
        }
        this.outputInterval = interval;
    }
    public void run() {
        Circle circle = new Circle(squareSide/2,squareSide/2,squareSide/2);
        MonteCarloPiData data = new MonteCarloPiData(circle);

        for (int i=0;i < N;i++) {
            if (i % outputInterval == 0) {
                System.out.println(data.estimatePi());
            }
            int x = (int)(Math.random() * squareSide);
            int y = (int)(Math.random() * squareSide);
            data.addPoint(new Point(x,y));
        }
    }
    public static void main(String[] args) {

        int squareSide = 800;
        int N = 10000000;

        MonteCarloExeperiment exp = new MonteCarloExeperiment(squareSide,N);
        exp.setOutputInterval(100000);
        exp.run();
    }
}

这里设置了每次产生点的数量,其它逻辑差不多,去掉了视图的绘制, 这次用产生1千万随机数的方法来测试。
run一下,可以看到


捕获.PNG

这次的结果,基本能精确到小数点后面4位、5位了,至此,计算过程就到这里了。

总结

这里使用蒙特卡洛算法,也就是产生更大的数据量是估算的值达到更精确。并且整个过程使用的是MVC的设计,再一次见识了计算机的魅力。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,591评论 6 501
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,448评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,823评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,204评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,228评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,190评论 1 299
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,078评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,923评论 0 274
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,334评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,550评论 2 333
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,727评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,428评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,022评论 3 326
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,672评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,826评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,734评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,619评论 2 354

推荐阅读更多精彩内容

  • 今日第二篇 淘金归淘金。 一代宗师里一句话说的有味道,宁在一思进,莫在一思停。 我理解的可能不知道歪在哪里去了。 ...
    aldenYang阅读 309评论 0 0
  • 十二要一份带汤的桂林米粉,而我在打着字,心里想着“桂林米粉汤,这样老板要是给我带一份汤怎么办,桂林米汤粉,什么东西...
    爱天蓝QG阅读 237评论 0 2
  • 朋友发来短信说,领导要把他调到北京,他在考虑要不要去。 我回道,去!为什么不去。 人们不是都说在一个城市就会渐...
    路中间的姑娘阅读 542评论 0 1
  • 不知道我有几个三年耗在网络的世界,但是我知道的是我再也不会开了电脑就爬上游戏没日没夜,这个曾经让我当来愉悦的乐子也...
    赵任慈i阅读 98评论 0 1
  • 嗨,你今年也32了呢!时间过的好快呀,还记得小时候,你为了让爸爸早点回家,嫌日子过的慢,将月份牌子撕到冬月,现在,...
    妮娜在灯火阑珊处阅读 926评论 0 0