单纯形法完全c语言程序

更新时间：2017/2/8 13:37:00  浏览量：566  手机版

单纯形法完全c语言程序,能运行

#include "math.h"

#include "stdio.h"

#define N 2

void paixu(p,n)

int n;

double p[];

{ int m,k,j,i;

double d;

k=0; m=n-1;

while (k

{ j=m-1; m=0;

for (i=k; i<=j; i++)

if (p>p[i+1])

{ d=p; p=p[i+1]; p[i+1]=d; m=i;}

j=k+1; k=0;

for (i=m; i>=j; i--)

if (p[i-1]>p)

{ d=p; p=p[i-1]; p[i-1]=d; k=i;}

}

return;

}

double mubiao(double *x)

{ double y;

y=x[1]-x[0]*x[0]; y=100.0*y*y;

y=y+(1.0-x[0])*(1.0-x[0]);

return(y);

}

main()

{ int i,j,k,l,m=0;

double c,xx[N+1][N],f0[N+1],f[N+1],x0[N]={1.2,1},x1[N],s=0.0; double a,b;

double xa[N],xb[N],xc[N],xe[N],xw[N],xr[N],xo[N];

double fr,fe,fw,fc,fo;

double aef=1.0,r=1.0,eps1=1.0e-30,eps2=1.0e-30,bt=0.5,rou=0.5; c=1.0;

b=(c/(N*sqrt(2)))*(sqrt(N+1)-1);

a=b+c/sqrt(2);

//printf("a=%13.7e b=%13.7e ",a,b);

//printf("\n");

//给xx[N][N+1]赋值,每一行构成单纯形的一个定点

//***********************

for(i=0;i

xx[0]=x0;

for(i=1;i

for(j=0;j

{if(j==i-1)

xx[j]=x0[j]+a;

else

xx[j]=x0[j]+b;

}

for (i=0;i

{for (j=0;j

printf("xx[%d][%d]%13.7e ",i,j,xx[j]);

printf("\n");

}

loop1:

//求单纯形的每个定点的函数值f0,f和x1是过渡数组

printf("\n");

printf("\n");

for(i=0;i

{for(j=0;j

x1[j]=xx[j];

f0=mubiao(x1);

f=mubiao(x1);

printf("f0[%d]=%13.7e f[%d]=%13.7e\n",i,f0,i,f);

}

printf("\n");

//比较f的大小，f[0]是最小值，f[N]是最大值

paixu(f,N+1);

for(i=N;i>=0;i--)

printf("f[%d]=%13.7e \n",i,f);

//找最好点和最坏点分别是哪一个点，即xx[][]的行数

for(i=0;i

{if(f0==f[0])

k=i;

if(f0==f[N])

l=i;

}

printf("最好点k=%d\n",k);

printf("最坏点l=%d\n",l);

//终止判断条件

printf("f[N]-f[0]=%13.7e \n",f[N]-f[0]);

if((f[N]-f[0])

{printf("迭代次数m=%d\n",m);

for(j=0;j

printf("optx[%d]=%13.7e\n",j,xx[k][j]);

printf("fmin=%13.7e\n",f[0]);

}

else

{ m=m+1;

//把xx[][]中最好点移到第一行，最坏点移到最后一行

for(j=0;j

{xb[j]=xx[k][j];

xx[k][j]=xx[0][j];

xx[0][j]=xb[j];

//

xw[j]=xx[l][j];

xx[l][j]=xx[N][j];

xx[N][j]=xw[j];

}

for (i=0;i

{for (j=0;j

printf("xx[%d][%d]=%13.7e ",i,j,xx[j]);

printf("\n");

}

//求除最坏点f[N]外其余点的中点xc[]

for(i=0;i

xa=0;

for(j=0;j

for(i=0;i

xa[j]=xa[j]+xx[j];}

xa[j]=xa[j]/N;

}

for(i=0;i

printf("xa[%d]=%13.7e xb[%d]=%13.7e xw[%d]=%13.7e \n",i,xa,i,xb,i,xw); //求xw[N]的反射点xr[N];

for(i=0;i

{xr=xa+aef*(xa-xw);

printf("xr[%d]=%13.7e ",i,xr);

}

printf("\n");

//求xr[N]的函数值fr

fr=mubiao(xr);

printf("fr=%13.7e \n",fr);

//判断xr与xb的好坏

if(fr<=f[0])

{for(i=0;i

{xe=xr+r*(xr-xa);

//printf("xe[%d]=%13.7e ",i,xe); }

printf("\n");

fe=mubiao(xe);

if(fe<=f[0])

for(j=0;j

xx[N][j]=xe[j];

else

for(j=0;j

xx[N][j]=xr[j];

goto loop1;

}

else

{ fw=f[N];

if(fr>=fw)

{for(i=0;i

xc=xa-bt*(xa-xw);

fc=mubiao(xc);

if(fc>=fw)

{for(i=1;i

for(j=0;j

xx[j]=xx[0][j]-rou*(xx[j]-xx[0][j]); goto loop1;

}

else

{for(j=0;j

xx[N][j]=xc[j];

goto loop1;

}

}

else

{if(fr>=fe)

{ for(i=0;i

xo=xa+bt*(xa-xw);

fo=mubiao(xo);

if(fo>=fr)

{for(i=1;i

for(j=0;j

xx[j]=xx[0][j]-rou*(xx[j]-xx[0][j]); goto loop1;

}