powell法c/c++程序

以下代码的vs2019工程文件打包点此下载去掉运行部分的注释即可

/*powell.cpp*/
#include "stdafx.h"
#include <stdio.h>
#include <math.h>
#include "Powell.h"
#include <iostream>
#define m 10              /*数组长度m  >=  维数n    */


float powell::powell_f(float x[])
{
	float result;
	//result = (x[0] - 2)*(x[0] - 2) + (x[1] - 3)*(x[1] - 3) + (x[2] - 4)*(x[2] - 4);
	result = x[0] * x[0] + 2 * x[1] * x[1] - 2 * x[0] * x[1] - 4 * x[0];
	return result;
}

/*多维进退法子程序*/
void powell::mjtf(int n, float x0[], float h, float s[], float a[], float b[])
{
	int i;
	float x1[m], x2[m], x3[m], powell_f1, powell_f2, powell_f3;
	for (i = 0; i<n; i++)         /*计算初始两试点*/
	{
		x1[i] = x0[i];
		x2[i] = x0[i] + h*s[i];
	}
	powell_f1 = powell_f(x1);
	powell_f2 = powell_f(x2);
	if (powell_f2 >= powell_f1)               /*判断搜索方向*/
	{                       /*搜索方向为反向，转身*/
		h = (-1)*h;
		for (i = 0; i<n; i++)
			x3[i] = x1[i];
		powell_f3 = powell_f1;
		for (i = 0; i<n; i++)
			x1[i] = x2[i];
		powell_f1 = powell_f2;
		for (i = 0; i<n; i++)
			x2[i] = x3[i];
		powell_f2 = powell_f3;
	}                       /*搜索方向为正向*/

	for (i = 0; i<n; i++)         /*计算第三试点*/
		x3[i] = x2[i] + h*s[i];
	powell_f3 = powell_f(x3);
	while (powell_f3<powell_f2)             /*判断是否未完成搜索*/
	{                       /*未完成，继续搜索*/
		h = 2 * h;
		for (i = 0; i<n; i++)
			x1[i] = x2[i];
		powell_f1 = powell_f2;
		for (i = 0; i<n; i++)
			x2[i] = x3[i];
		powell_f2 = powell_f3;
		for (i = 0; i<n; i++)
			x3[i] = x2[i] + h*s[i];
		powell_f3 = powell_f(x3);
	}                       /*已完成*/
	for (i = 0; i<n; i++)         /*输出初始搜索区间*/
	{
		if (x1[i]<x3[i])
		{
			a[i] = x1[i];
			b[i] = x3[i];
		}
		else
		{
			a[i] = x3[i];
			b[i] = x1[i];
		}
	}
}

/*多维黄金分割法子程序*/
void powell::mhjfgf(int n, float a[], float b[], float flag, float x[])
{
	int i;
	float x1[m], x2[m], powell_f1, powell_f2, sum;
	for (i = 0; i<n; i++)              /*计算初始两试点*/
		x1[i] = b[i] - (float)0.618*(b[i] - a[i]);
	powell_f1 = powell_f(x1);
	for (i = 0; i<n; i++)
		x2[i] = a[i] + (float)0.618*(b[i] - a[i]);
	powell_f2 = powell_f(x2);
	do
	{
		if (powell_f1 <= powell_f2)                  /*判断消去区间*/
		{                          /*消去右*/
			for (i = 0; i<n; i++)
				b[i] = x2[i];
			for (i = 0; i<n; i++)
				x2[i] = x1[i];
			powell_f2 = powell_f1;
			for (i = 0; i<n; i++)
				x1[i] = b[i] - (float)0.618*(b[i] - a[i]);
			powell_f1 = powell_f(x1);
		}
		else
		{                          /*消去左*/
			for (i = 0; i<n; i++)
				a[i] = x1[i];
			for (i = 0; i<n; i++)
				x1[i] = x2[i];
			powell_f1 = powell_f2;
			for (i = 0; i<n; i++)
				x2[i] = a[i] + (float)0.618*(b[i] - a[i]);
			powell_f2 = powell_f(x2);
		}
		sum = 0;
		for (i = 0; i<n; i++)
			sum += (b[i] - a[i])*(b[i] - a[i]);
	} while (sqrt(sum)>flag*0.1);
	for (i = 0; i<n; i++)
		x[i] = (float)0.5*(b[i] + a[i]);
}

/*鲍威尔法子程序*/
void powell::mbwef(int n, float x0[], float h, float flag, float a[], float b[], float x[])
{
	int i, j, k, r;
	float x1[m], x11[m], x2[m], powell_f0, powell_f1, powell_f2, fn[m], s[m][m], sum;
	for (i = 0; i<n; i++)//构造一个单位对角矩阵s
		for (k = 0; k<n; k++)
			if (i == k)
				s[i][k] = 1;
			else
				s[i][k] = 0;
	k = 0;
	while (1)
	{
		for (i = 0; i<n; i++)
		{
			x1[i] = x0[i];
			x11[i] = x1[i];
		}
		for (i = 0; i<n; i++)
		{
			mjtf(n, x1, h, s[i], a, b);
			mhjfgf(n, a, b, flag, x1);
			fn[i] = powell_f(x11) - powell_f(x1);
		}
		for (i = 0; i<n; i++)
			x2[i] = 2 * x1[i] - x0[i];
		for (i = 1; i<n; i++)
			if (fn[0]<fn[i])
			{
				fn[0] = fn[i];
				r = i;
			}
			else
				r = 0;
		powell_f0 = powell_f(x0);
		powell_f1 = powell_f(x1);
		powell_f2 = powell_f(x2);
		if (powell_f2 >= powell_f0 || (powell_f0 - 2 * powell_f1 + powell_f2)*(powell_f0 - powell_f1 - fn[0])*(powell_f0 - powell_f1 - fn[0]) >= 0.5*fn[0] * (powell_f0 - powell_f2)*(powell_f0 - powell_f2))
		{
			sum = 0;
			for (i = 0; i<n; i++)
				sum += (x1[i] - x0[i])*(x1[i] - x0[i]);
			if (powell_f1 <= powell_f2)
				for (i = 0; i<n; i++)
					x0[i] = x1[i];
			else
				for (i = 0; i<n; i++)
					x0[i] = x2[i];
		}
		else
		{
			for (i = r; i<n - 1; i++)
				for (j = 0; j<n; j++)
					s[i][j] = s[i + 1][j];
			for (i = 0; i<n; i++)
				s[n - 1][i] = x1[i] - x0[i];
			mjtf(n, x1, h, s[n - 1], a, b);
			mhjfgf(n, a, b, flag, x1);
			sum = 0;
			for (i = 0; i<n; i++)
				sum += (x1[i] - x0[i])*(x1[i] - x0[i]);
			for (i = 0; i<n; i++)
				x0[i] = x1[i];
		}
		if (sqrt(sum) <= flag)
			break;
		else
			k += 1;
	}
	for (i = 0; i<n; i++)
		x[i] = x1[i];
}

/*鲍威尔法主程序*/

int main(void)
{
	int i, n = 2;
	float h, flag, x0[m], a[m], b[m], x[m];
	printf("\n<鲍威尔法>\n");
	/*printf("请输入维数:\n");
	scanf_s("%d", &n);*/
	printf("请输入初始点:");
	for (i = 0; i<n; i++)
	{
		printf("\nx0[%d]=", i);
		scanf_s("%f", &x0[i]);
	}
	printf("\n请输入初始步长:\n");
	scanf_s("%f", &h);
	printf("\n请输入精度:\n");
	scanf_s("%f", &flag);
	powell one;
	one.mbwef(n, x0, h, flag, a, b, x);
	printf("\n极小点坐标为:\n");
	for (i = 0; i<n; i++)
		printf("x[%d]=%f\n", i, x[i]);
	printf("\n极小值为:\n%f\n", one.powell_f(x));
	while (1)
	{
		;
	}

}

#ifndef _powell_h_
#define _powell_h_

class powell
{
public:
	float powell_f(float x[]);
	void mjtf(int n, float x0[], float h, float s[], float a[], float b[]);
	void mhjfgf(int n, float a[], float b[], float flag, float x[]);
	void mbwef(int n, float x0[], float h, float flag, float a[], float b[], float x[]);
private:

};

#endif // !
#pragma once