数值概率算法

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

原文地址:https://www.cnblogs.com/cy0628/p/14007582.html