zoukankan      html  css  js  c++  java
  • 数值概率算法

    1、用随机投点法计算pi

      设有一半径为r的圆及其外切四边形。向该正方形随机地投掷n个点。设落入圆内的点数为k。由于所投入的点在正方形上均匀分布,因而所投入的点落入圆内的概率为(PI * pow(r,2)) / (4 * pow(r,2)) = PI / 4 。所以当n足够大时,k与n之比就逼近这一概率。从而,PI 约等于 (4*k)/n.如下图:

    实现:

    #include <iostream>
    #include <cstdlib>
    #include <limits>
    using namespace std;
    
    // 获得0-1之间的随机数
    double get_random_num ()
    {
        return (double)rand () / RAND_MAX ;
    }
    // 用随机投点法计算 PI
    double darts (int n)
    {
        int k = 0 ;
        for (int i = 0; i < n; ++ i) {
            double x = get_random_num() ;
            double y = get_random_num() ;
            if ((x * x + y * y) <= 1.0) {
                ++ k ;
            }
        }
        return (4 * k) / (double)n ;
    }
    int main()
    {
        cout << darts (200000000) << endl ;
    }
    View Code

    实现结果:

    2.计算定积分

    设f(x)是[0,1]上的连续函数,且0 <= f(x) <= 1。

    需要计算的积分为,积分I等于图中的面积G。

    在图所示单位正方形内均匀地作投点试验,则随机点落在曲线下面的概率为

    假设向单位正方形内随机地投入n个点(xi,yi)。如果有m个点落入

    G内,则随机点落入G内的概率

    假定 f(x) = x * x (0 <= x <= 1)计算定积分

    实现:

    #include <iostream>
    #include <cstdlib>
    using namespace std;
    
    // 获得0-1之间的随机数
    double get_random_num ()
    {
        return (double)rand () / RAND_MAX ;
    }
    // 用随机投点法计算 y = pow(x,2)定积分
    double darts (int n)
    {
        int k = 0 ;
        for (int i = 0; i < n; ++ i) {
            double x = get_random_num() ;
            double y = get_random_num() ;
            if (y <= x * x) {
                ++ k ;
            }
        }
        return k / (double)n ;
    }
    int main()
    {
        cout << darts (10000000) << endl ;
        return 0;
    }
    View Code

    运行结果:

    3.解非线性方程组

    求解下面的非线性方程组

     

    其中,x1,x2,…,xn是实变量,fi是未知量x1,x2,…,xn的非线性实函数。要求确定上述方程组在指定求根范围内的一组解

     

    问题分析

         解决这类问题有多种数值方法,如:牛顿法、拟牛顿法、粒子群算法等。最常用的有线性化方法和求函数极小值方法。为了求解所给的非线性方程组,构造一目标函数

         式中,x=(x1,x2,……xn)。易知,函数取得零点即是所求非线性方程组的一组解。

     求解思路

        在指定求根区域D内,选定一个随机点x0作为随机搜索的出发点。在算法的搜索过程中,假设第j步随机搜索得到的随机搜索点为xj。在第j+1步,计算出下一步的随机搜索增量dxj。从当前点xj依dxj得到第j+1步的随机搜索点。当x<时,取为所求非线性方程组的近似解。否则进行下一步新的随机搜索过程。

    实现:

    //随机化算法 解线性方程组
    #include "stdafx.h"
    #include "RandomNumber.h"
    #include <iostream>
    using namespace std;
     
    bool NonLinear(double *x0,double *dx0,double *x,double a0,
                        double epsilon,double k,int n,int Steps,int M);
    double f(double *x,int n);
     
    int main()
    {
        double  *x0,                //根初值
                *x,                    //
                *dx0,                //增量初值
                a0 = 0.0001,            //步长
                epsilon = 0.01,        //精度
                k = 1.1;            //步长变参
        int n = 2,                    //方程个数
            Steps = 10000,            //执行次数
            M = 1000;                //失败次数
     
        x0 = new double[n+1];
        dx0 = new double[n+1];
        x = new double[n+1];
     
        //根初值
        x0[1] = 0.0;
        x0[2] = 0.0;
     
        //增量初值
        dx0[1] = 0.01;
        dx0[2] = 0.01;
     
        cout<<"原方程组为:"<<endl;
        cout<<"x1-x2=1"<<endl;
        cout<<"x1+x2=3"<<endl;
     
        cout<<"此方程租的根为:"<<endl;
     
        bool flag = NonLinear(x0,dx0,x,a0,epsilon,k,n,Steps,M);
        while(!flag)
        {        
            flag = NonLinear(x0,dx0,x,a0,epsilon,k,n,Steps,M);
        }    
        for(int i=1; i<=n; i++)
        {
            cout<<"x"<<i<<"="<<x[i]<<" ";
        }
        cout<<endl;
     
        return 0;
    }
     
    //解非线性方程组的随机化算法
    bool NonLinear(double *x0,double *dx0,double *x,double a0,
                        double epsilon,double k,int n,int Steps,int M)
    {
        static RandomNumber rnd;
        bool success;            //搜索成功标志
        double *dx,*r;
     
        dx = new double[n+1];    //步进增量向量
        r = new double[n+1];    //搜索方向向量
        int mm = 0;                //当前搜索失败次数
        int j = 0;                //迭代次数
        double a = a0;            //步长因子
     
        for(int i=1; i<=n; i++)
        {
            x[i] = x0[i];
            dx[i] = dx0[i];
        }
     
        double fx = f(x,n);        //计算目标函数值
        double min = fx;        //当前最优值
     
        while(j<Steps)
        {
            //(1)计算随机搜索步长
            if(fx<min)//搜索成功
            {
                min = fx;
                a *= k;
                success = true;
            }
            else//搜索失败
            {
                mm++;
                if(mm>M)
                {
                    a /= k;
                }
                success = false;
            }
     
            if(min<epsilon)
            {
                break;
            }
     
            //(2)计算随机搜索方向和增量
            for(int i=1; i<=n; i++)
            {
                r[i] = 2.0 * rnd.fRandom()-1;
            }
     
            if(success)
            {
                for(int i=1; i<=n; i++)
                {
                    dx[i] = a * r[i];
                }
            }
            else
            {
                for(int i=1; i<=n; i++)
                {
                    dx[i] = a * r[i] - dx[i];
                }
            }
     
            //(3)计算随机搜索点
            for(int i=1; i<=n; i++)
            {
                x[i] += dx[i];
            }
     
            //(4)计算目标函数值
            fx = f(x,n);
            j++;
        }    
     
        if(fx<=epsilon)
        {
            return true;
        }
        else
        {
            return false;
        }
    }
     
    double f(double *x,int n)
    {
        return (x[1]-x[2]-1)*(x[1]-x[2]-1)
                +(x[1]+x[2]-3)*(x[1]+x[2]-3);
    }
    View Code

    运行结果:

    参考:王晓东《算法设计与分析》第二版

                https://www.cnblogs.com/chinazhangjie/archive/2010/11/11/1874924.html

                https://blog.csdn.net/liufeng_king/article/details/9029091

  • 相关阅读:
    MC9S12 硬件设计
    ESD
    选用与使用稳压二极管的介绍
    MOSFET 栅极电阻作用及其选型
    orcad常用库文件介绍
    开关电源和LDO的区别
    续流二极管的作用及选型
    为什么大电容滤低频小电容滤高频的问题
    Java常用API——时间类
    Idea问题:“marketplace plugins are not loaded”解决方案
  • 原文地址:https://www.cnblogs.com/cy0628/p/14007582.html
Copyright © 2011-2022 走看看