数值概率算法

1、用随机投点法计算pi值

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

实现：

#include 
#include 
#include 
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 
#include 
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 
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<<"原方程组为："<    cout<<"x1-x2=1"<    cout<<"x1+x2=3"< 
    cout<<"此方程租的根为："< 
    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"<    }
    cout< 
    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    {
        //(1)计算随机搜索步长
        if(fx        {
            min = fx;
            a *= k;
            success = true;
        }
        else//搜索失败
        {
            mm++;
            if(mm>M)
            {
                a /= k;
            }
            success = false;
        }
 
        if(min        {
            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

运行结果：

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