最优化方法课程实验报告

下载可编辑

for(k=0;k

for(j=0;j

c[k]+=a[k][j]*b[j]; } }

cout<<\矩阵的逆矩阵和梯度向量相乘，两矩阵相乘结果为:\for(i=0;i

cout<

cout<<\两个矩阵阶数不符，不能相乘\}

//开始计算minf(Xk+tPk)时的步长t的值，由于这是n元二次函数所以minf(t)是关于t>0的二次函数，

先求二次方程a,b,c系数,用一阶导为零。

void abc(double x[n],double p[n],double f_xs[n+n+1+(n-1)*n/2],double t[3])//t[3]中返回的是

a,b,c的系数值 {

int i,j;

t[0]=t[1]=t[2]=0;

t[0]=fun(p,f_xs)-f_xs[2*n];

for(i=n;i<2*n;i++) {

t[0]-=f_xs[i]*p[i%n]; }

for(i=0;i

t[1]+=2*f_xs[i]*x[i]*p[i]; }

for(;i<2*n;i++) {

t[1]+=f_xs[i]*p[i%n]; }

for(i=0;i

for(j=i+1;j

t[1]+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*(x[i]*p[j]+x[j]*p[i]); } }

t[2]=fun(x,f_xs);

cout<

\}

//二次函数一阶导为零计算t的值，t>0 double t_c(double t[3]) {

cout<<\用一阶导求得步长为：t=\return -t[1]/t[0]/2; }

void main() {

double f_xs[n+n+1+(n-1)*n/2]={4,2,9,-3,16,0};//正定二次函数的各项系数，第一个为：X1^2系数，

第二个为：X2^2系数，第三个为：X1系数，第四个为:X2系数，第五个为：常数项，第六个为：X1X2系数；

.专业.整理.

下载可编辑

//n元的系数存放类推

double x[n]={0,0};//函数自变量初值 double f0;//函数值

double g0[n];//梯度的值

double Q[n][n+1];//正定二次函数对应的实对称矩阵 double c=0.01;//H准则限值

double Hesse[n][n];//Hesse矩阵

double inv_Hesse[n][n];//Hesse矩阵的逆

double t[3];//返回求minf()时t的二次函数的a,b,c的系数值 double t_bc;//步长

double p[n];//保存矩阵相乘结果_加负号后变成下降方向 int i,tap=0;//迭代次数

Q_fun(f_xs,Q);//计算正定二次函数对应的实对称矩阵ok! f0=fun(x,f_xs);//求函数初值ok! D_fun(x,Q,g0);//返回梯度的初值 while(H(g0,c)==0) {

Hesse_fun(Q,Hesse);//返回Hesse矩阵的值？？一个二次函数的Hesse矩阵是变的还是不变 inv(Hesse,inv_Hesse);//返回Hesse矩阵的逆矩阵 xiang_cheng(inv_Hesse,g0,p); for(i=0;i

abc(x,p,f_xs,t);//开始计算minf(Xk+tPk)时的步长t的值， t_bc=t_c(t);//求一阶导来计算t for(i=0;i

x[i]=x[i]+t_bc*p[i]; }

tap++;

f0=fun(x,f_xs); D_fun(x,Q,g0); }

cout<<\修正newton法\

cout<<\函数f(x1,x2)=4(X1+1)^2+2(X2-1)^2+X1+X2+10.的极小点为：\cout<<\自变量取值为：\for(i=0;i

cout<<\}

cout<<\迭代次数为：\}

附录9

#include #include #define N 10

#define eps pow(10,-6)

double f(double x[],double p[],double t) { double s;

s=pow(x[0]+t*p[0],2)+4*pow(x[1]+t*p[1],2); return s; }

void sb(double *a,double *b,double x[],double p[]) {

double t0,t1,t,h,alpha,f0,f1; int k=0;

t0=2.5; /*初始值*/ h=1; /*初始步长*/ alpha=2; /*加步系数*/ .专业.整理.

下载可编辑

f0=f(x,p,t0);

t1=t0+h; f1=f(x,p,t1); while(1) { if(f1

else { if(k==0) {h=-h;t=t1;} else { *a=tt1?t:t1; break; } }

t1=t0+h; f1=f(x,p,t1); } }

double hjfg(double x[],double p[]) {

double beta,t1,t2,t; double f1,f2; double a=0,b=0;

double *c,*d; c=&a,d=&b;

sb(c,d,x,p);/*调用进退法搜索区间*/

printf(\ printf(\ beta=(sqrt(5)-1.0)/2; t2=a+beta*(b-a); f2=f(x,p,t2); t1=a+b-t2; f1=f(x,p,t1); while(1) {

if(fabs(t1-t2)

{ if(f1

{ t=(t1+t2)/2; b=t2;

t2=t1; f2=f1;

t1=a+b-t2; f1=f(x,p,t1); }

else { a=t1; t1=t2;

f1=f2;

t2=a+beta*(b-a); f2=f(x,p,t2); } } }

t=(t1+t2)/2; return t; }

void gtd() {

double x[N],g[N],p[N],t=0,f0,mod1=0,mod2=0,nanda=0; int i,k,n;

printf(\请输入函数的元数值n=\

.专业.整理.

下载可编辑

scanf(\

printf(\请输入初始值:\\n\

for(i=0;i

mod1=sqrt(pow(g[0],2)+pow(g[1],2));/*求梯度的长度*/ if(mod1>eps) { p[0]=-g[0]; p[1]=-g[1]; k=0;

while(1) { t=hjfg(x,p);/*调用黄金分割法求t的值*/ printf(\ x[0]=x[0]+t*p[0]; x[1]=x[1]+t*p[1]; g[0]=2*x[0]; g[1]=50*x[1]; /*printf(\ mod2=sqrt(pow(g[0],2)+pow(g[1],2)); /*求梯度的长度*/ if(mod2<=eps) break; else { if(k+1==n) { g[0]=2*x[0]; g[1]=50*x[1]; p[0]=-g[0]; p[1]=-g[1]; k=0; } else { nanda=pow(mod2,2)/pow(mod1,2); printf(\ p[0]=-g[0]+nanda*p[0]; p[1]=-g[1]+nanda*p[1]; mod1=mod2; k++; } } printf(\ } }

printf(\最优解为x1=%lf,x2=%lf\ printf(\最终的函数值为%lf\}

main() {

gtd(); }

附录10

#include #include #include

int const n=2;//正定二次函数的自变量个数

double fun(double x[n],double f_xs[n+n+1+(n-1)*n/2]);//输入变量为函数自变量初值 void Q_fun(double f_xs[n+n+1+(n-1)*n/2],double Q[n][n+1]); void D_fun(double x[n],double Q[n][n+1],double g0[n]); .专业.整理.