共轭梯度法的C+=语言实现

共轭梯度法的C+=语言实现
---------------------基于成都理工大学最优化教材#include <iostream.h>
#include "Matrix.h"
#include<LIMITS.H>
#define MAX_M 2048
#include<math.h>
#define beta (sqrt(5)-1)/2
#define eclipse 0.01
CMatrix hesse()
{
CMatrix temp(2,2);
temp.A[0][0] = 2;
temp.A[0][1] = 0;
temp.A[1][0] = 0;
temp.A[1][1] = 8;
return temp;
}
double fun(double t,CMatrix &X1,CMatrix &X2)
{
double x1 = X1.A[0][0] + t*X2.A[0][0];
double x2 = X1.A[1][0] + t*X2.A[1][0];
return x1*x1 + 4*x2*x2;
}
double Max_Num(double a,double b)
{
if(a>b)
{
return a;
}
else
{
return b;
}
}
double Min_Num(double a,double b)
{
if(a<b)
{
return a;
}
else
{
return b;
}
}
void StepAdding_Search(double &a,double &b,CMatrix &mt1,CMatrix &mt2) {
double t[MAX_M]={0};
double h[MAX_M]={0};
double f[MAX_M]={0};
double result=0;
double p=2.1;
t[0]=0;
h[0]=1;
f[0]=fun(t[0],mt1,mt2);
for(int k=0;k<MAX_M-1;k++)
{
t[k+1]=t[k]+h[k];
f[k+1]=fun(t[k+1],mt1,mt2);
if(f[k+1]<f[k])
{
h[k+1]=p*h[k];
result=t[k];
t[k]=t[k+1];
f[k]=f[k+1];
}
else if(k == 0)
{
h[k]=-h[k];
result=t[k+1];
}
else
{
a=Min_Num(result,t[k+1]);
b=Max_Num(result,t[k+1]);
}
}
}
double Fibonacci(double &a,double &b,CMatrix mt1,CMatrix mt2)
{
double t2,t1,result_1,result_2;
t2 = a + beta*(b-a);
result_2 = fun(t2,mt1,mt1);
while(true)
{
t1 = a+b-t2;
result_1 = fun(t1,mt1,mt2);
if(fabs(t1-t2)<eclipse)
{
return (t1+t2)/2;
}
else if(result_1 < result_2)
{
b = t1;
t1 = t2;
result_1 = result_2;
t2 = a+beta*(b-a);
result_2 = fun(t2,mt1,mt2);
}
else if(result_1 == result_2)
{
a = t1;
b = t2;
result_2 = fun(a,mt1,mt2);
}
else
{
a = t2;
t2 = t1;
result_2 = result_1;
}
}
}
double distance(CMatrix &X1,CMatrix &X2) {
double x1 = X1.A[0][0] - X2.A[0][0];
double x2 = X1.A[1][0] - X2.A[1][0];
return x1*x1 + x2*x2;
}
CMatrix diff_fun(CMatrix &mt)
{
CMatrix temp(2,1);
temp.A[0][0] = 2*mt.A[0][0];
temp.A[1][0] = 8*mt.A[1][0];
return temp;
}
double fanshu(CMatrix mt)
{
return sqrt(mt.A[0][0]*mt.A[0][0] + mt.A[1][0]*mt.A[1][0]);
}
void main()
{
int i =0;
int n = 2000;
CMatrix temp_X1;
CMatrix X1(2,1);
X1.A[0][0] = 1;
X1.A[1][0] = 1;
CMatrix m_temp = hesse();
CMatrix n_temp = diff_fun(X1);
CMatrix X2 = (m_temp.Inverse()* n_temp).Matrix_Multiply(-1);
double a,b;
StepAdding_Search(a,b,X1,X2);
double m = Fibonacci(a,b,X1,X2);
temp_X1 = X1;
CMatrix X3 = X1 + X2.Matrix_Multiply(m);
while(distance(X1,X3) > eclipse)
{
X1 = X3;
if(i+1 == n)
{
i = 0;
X2 = m_temp.Inverse()* diff_fun(X1);
}
else
{
m = (fanshu(temp_X1)*fanshu(temp_X1))/(fanshu(X1)*fanshu(X1));
X2 = diff_fun(X1).Matrix_Multiply(-1) + X2.Matrix_Multiply(m);
}
X1 = X3;
X2 = m_temp.Inverse()* diff_fun(X1);
StepAdding_Search(a,b,X1,X2);
m = Fibonacci(a,b,X1,X2);
X3 = X1 + X2.Matrix_Multiply(m);
i++;
}
cout<<"求解的结果是:"<<endl;
cout<<"("<<X1.A[0][0]<<","<<X1.A[1][0]<<")"<<endl;
cout<<"最优解是:"<<endl;
cout<<fun(m,X1,X2)<<endl;
}
（注：本资料素材和资料部分来自网络，仅供参考。

请预览后才下载，期待您的好评与关注！）。