7.2.3 数值随机化算法解非线性方程组

关于问题描述、问题分析和问题求解思路等可以参考《计算机算法设计与分析》7.2.3节或者参考博客：http://m.blog.csdn.net/liufeng_king/article/details/9029091和http://www.cnblogs.com/tanky_woo/archive/2010/12/10/1901780.html

我的代码也是根据来自与这两篇博客，只是加上了一些显示过程的代码和一些流程的注释，还修改了一处小bug。
//随机化算法 解线性方程组

#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 t = 1;            //用来表示函数内程序的循环执行次数
int main()
{
double  *x0,                //根初值
    *x,                 //根
    *dx0,               //增量初值
    a0 = 0.001,         //步长
    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^2-x2^2=1" << endl;
cout << "x1^2+x2=3" << endl;

cout << "*****************************************************" << endl;
cout << "第一次搜索开始" << endl;
cout << "变量的初始值:" << x0[1] << "和" << x0[2] << endl;

cout << "步进因子的初始值：" << a0 << endl;
cout << "*****************************************************" << 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);
}

cout << "此方程租的根为：" << endl;
for (int i = 1; i <= n; i++)
{
    cout << "x" << i << "=" << x[i] << " ";
}
cout << endl;
cout <<  t << "次" << "搜索" << endl;
system("pause");
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;
        if (t < 15)
        {
            cout <<"第"<<t<<"次"<< "搜索成功"  << endl;
            cout << "目标函数的值：" << fx << endl;
            cout << "*****************************************************" << endl;
            cout << "第" << t+1 << "次" << "搜索开始" << endl;

            cout << "改变随机搜索的步长：" << a << endl;

        }
    }
    else//搜索失败
    {
        mm++;
        if (mm%M==0)        //搜索失败次数大于M次后，减小步长因子，就是用更大的精度搜索，每次搜索的跨度变小，这里原本为mm>M,容易陷入搜索步长越来越小的麻烦（多次测试都是陷入死循环），所以改成了这样子。而且你还可以设置下限，防止无限减小步长。
        {
            a /= k;
            //if (a <= 0.005)
            //  a = 0.005;
            cout << "当搜索失败次数大于1000次后，改变随机搜索的步长：" << a << endl;
            cout << "*****************************************************" << endl;

        }
        success = false;
        if (t < 15)
        {

            cout << "第" << t << "次" << "搜索失败"  << endl; 
            cout << "目标函数的值：" << fx<<"  上一次目标函数值："<<min << endl;
            cout << "*****************************************************" << endl;
            cout << "第" << t + 1 << "次" << "搜索开始" << endl;


        }
    }

    if (min<epsilon)         //这一次的值小于精度，则搜索完成，当前的变量值就是方程组的解
    {
        break;
    }

    //(2)计算随机搜索方向和增量
    for (int i = 1; i <= n; i++) 
    {
        r[i] = 2.0 * rnd.fRandom() - 1;    // 对每一个变量产生一个-1到1的随机数，作为搜索方向，这个在fRandom函数参考博客里
    }


    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];
    }

    if (t < 15)
    {
        if (success)
        {
            cout << "搜索成功情况下，继续在这一层的变量基础随机搜索：" << endl;

            cout << "dx1=" << dx[1] << endl;
            if (dx[1]>0) cout << "方向为正" << endl;
            else  cout << "方向为负" << endl;
            cout << "x1=" << x[1] << endl;

            cout << "dx2=" << dx[2] << endl;
            if (dx[2]>0) cout << "方向为正" << endl;
            else  cout << "方向为负" << endl;
            cout << "x2=" << x[2] << endl;
            cout << "*****************************************************" << endl;
        }
        else
        {
            cout << "搜索失败情况下，退回到上一层变量基础随机搜索：" << endl;
            //cout << "a*r1=" << a*r[1] << endl;
            cout << "dx1=" << dx[1] << endl;
            if (dx[1]>0) cout << "方向为正" << endl;
            else  cout << "方向为负" << endl;
            cout << "x1=" << x[1] << endl;
            cout << endl;
            //cout << "a*r2=" << a*r[2] << endl;
            cout << "dx2=" << dx[2] << endl;
            if (dx[2]>0) cout << "方向为正" << endl;
            else  cout << "方向为负" << endl;
            cout << "x2=" << x[2] << endl;
                cout <<"*****************************************************"<< endl;
        }

    }

    //(4)计算目标函数值
    fx = f(x, n);
    t++;
    j++;
    }

if (fx <= epsilon)
{
    return true;
}
else
{
    return false;
}
}

double f(double *x, int n)//计算目标函数值，这里被写死了（是两个方程的平方和）。
{
return (x[1] * x[1] - x[2] * x[2] - 1)*(x[1] * x[1] - x[2] * x[2] - 1)
    + (x[1] * x[1] + x[2] - 3)*(x[1] * x[1] + x[2] - 3);
}

这段代码是参考博客里的代码，需要的函数都能在里面找见。