ch10.h

数值处理算法程序

          if ((i!=r)&&(f[i]>fg))
            { g=i; fg=f[i];}
        for (j=0; j<=n-1; j++)
          { xf[j]=0.0;
            for (i=0; i<=n; i++)
              if (i!=r)
                xf[j]=xf[j]+xx[j*(n+1)+i]/nn;
            xt[j]=2.0*xf[j]-xx[j*(n+1)+r];
          }
        ft=jjsimf(xt,n);
        if (ft<f[l])
          { for (j=0; j<=n-1; j++)
              xf[j]=(1.0+u)*xt[j]-u*xf[j];
            ff=jjsimf(xf,n);
            if (ff<f[l])
              { for (j=0; j<=n-1; j++)
                  xx[j*(n+1)+r]=xf[j];
                f[r]=ff;
              }
            else
              { for (j=0; j<=n-1; j++)
                  xx[j*(n+1)+r]=xt[j];
                f[r]=ft;
              }
          }
        else
          { if (ft<=f[g])
              { for (j=0; j<=n-1; j++)
                  xx[j*(n+1)+r]=xt[j];
                f[r]=ft;
              }
            else 
              { if (ft<=f[r])
                  { for (j=0; j<=n-1; j++)
                      xx[j*(n+1)+r]=xt[j];
                    f[r]=ft;
                  }
                for (j=0; j<=n-1; j++)
                  xf[j]=v*xx[j*(n+1)+r]+(1.0-v)*xf[j];
                ff=jjsimf(xf,n);
                if (ff>f[r])
                  for (i=0; i<=n; i++)
                    { for (j=0; j<=n-1; j++)
                        { xx[j*(n+1)+i]=(xx[j*(n+1)+i]+
					xx[j*(n+1)+l])/2.0;
                          x[j]=xx[j*(n+1)+i]; xe[j]=x[j];
                        }
                      fe=jjsimf(xe,n); f[i]=fe;
                    }
                else
                  { for (j=0; j<=n-1; j++)
                      xx[j*(n+1)+r]=xf[j];
                    f[r]=ff;
                  }
              }
          }
        ff=0.0; ft=0.0;
        for (i=0; i<=n; i++)
          { ff=ff+f[i]/(1.0+nn);
            ft=ft+f[i]*f[i];
          }
        ft=(ft-(1.0+n)*ff*ff)/nn;
      }
    for (j=0; j<=n-1; j++)
      { x[j]=0.0;
        for (i=0; i<=n; i++)
          x[j]=x[j]+xx[j*(n+1)+i]/(1.0+nn);
        xe[j]=x[j];
      }
    fe=jjsimf(xe,n); x[n]=fe;
    free(xt); free(xf); free(xe);
    return(kk);
}
/////////////////////////////////////////////////////////////
static double rn(double* rr)
{ 
	int m;
    double y,s,u,v;
    s=65536.0; u=2053.0; v=13849.0;
    *rr=u*(*rr)+v; m=*rr/s; *rr=*rr-m*s;
    y=*rr/s;
    return(y);
}
 
int jcplx(int n,int m,double a[],double b[],double alpha,double eps,
		  double x[],double xx[],int k,double (*jcplxf)(double x[],int n),
		  void (*jcplxs)(int n,int m,double x[],double c[],double d[],double w[]))//求约束条件下n维极值的复形调优法
{ 
	//extern double jcplxf();
    //extern void jcplxs();
    //double rn();
    int r,g,i,j,it,kt,jt,kk;
    double fj,fr,fg,z,rr,*c,*d,*w,*xt,*xf;
    c=(double*)malloc(m*sizeof(double));
    d=(double*)malloc(m*sizeof(double));
    w=(double*)malloc(m*sizeof(double));
    xt=(double*)malloc(n*sizeof(double));
    xf=(double*)malloc(n*sizeof(double));
    rr=0.0;
    for (i=0; i<=n-1; i++)
      xx[i*n*2]=x[i];
    xx[n*n*2]=jcplxf(x,n);
    for (j=1; j<=2*n-1; j++)
      { for (i=0; i<=n-1; i++)
	  { xx[i*n*2+j]=a[i]+(b[i]-a[i])*rn(&rr);
            x[i]=xx[i*n*2+j];
          }
        it=1;
        while (it==1)
          { it=0; r=0; g=0;
            while ((r<n)&&(g==0))
              { if ((a[r]<=x[r])&&(b[r]>=x[r])) r=r+1;
                else g=1;
              }
            if (g==0)
              { jcplxs(n,m,x,c,d,w);
                r=0;
                while ((r<m)&&(g==0))
                  { if ((c[r]<=w[r])&&(d[r]>=w[r])) r=r+1;
                    else g=1;
                  }
              }
            if (g!=0)
              { for (r=0; r<=n-1; r++)
                  { z=0.0;
                    for (g=0; g<=j-1; g++)
                      z=z+xx[r*n*2+g]/(1.0*j);
                    xx[r*n*2+j]=(xx[r*n*2+j]+z)/2.0;
                    x[r]=xx[r*n*2+j];
                  }
                it=1;
              }
            else xx[n*n*2+j]=jcplxf(x,n);
          }
      }
    kk=1; it=1;
    while (it==1)
      { it=0;
        fr=xx[n*n*2]; r=0;
        for (i=1; i<=2*n-1; i++)
          if (xx[n*n*2+i]>fr)
            { r=i; fr=xx[n*n*2+i];}
        g=0; j=0; fg=xx[n*n*2];
        if (r==0)
          { g=1; j=1; fg=xx[n*n*2+1];}
        for (i=j+1; i<=2*n-1; i++)
          if (i!=r)
            if (xx[n*n*2+i]>fg)
              { g=i; fg=xx[n*n*2+i];}
        for (i=0; i<=n-1; i++)
          { xf[i]=0.0;
            for (j=0; j<=2*n-1; j++)
              if (j!=r)
                xf[i]=xf[i]+xx[i*n*2+j]/(2.0*n-1.0);
            xt[i]=(1.0+alpha)*xf[i]-alpha*xx[i*n*2+r];
          }
        jt=1;
        while (jt==1)
          { jt=0;
            z=jcplxf(xt,n);
            while (z>fg)
              { for (i=0; i<=n-1; i++)
                  xt[i]=(xt[i]+xf[i])/2.0;
                z=jcplxf(xt,n);
              }
            j=0;
            for (i=0; i<=n-1; i++)
              { if (a[i]>xt[i])
                  { xt[i]=xt[i]+0.000001; j=1;}
                if (b[i]<xt[i])
                  { xt[i]=xt[i]-0.000001; j=1;}
              }
            if (j!=0) jt=1;
            else
              { jcplxs(n,m,xt,c,d,w);
                j=0; kt=1;
                while ((kt==1)&&(j<m))
                  { if ((c[j]<=w[j])&&(d[j]>=w[j])) j=j+1;
                    else kt=0;
                  }
                if (j<m)
                  { for (i=0; i<=n-1; i++)
                      xt[i]=(xt[i]+xf[i])/2.0;
                    jt=1;
                  }
              }
          }
        for (i=0; i<=n-1; i++)
          xx[i*n*2+r]=xt[i];
        xx[n*n*2+r]=z;
        fr=0.0; fg=0.0;
        for (j=0; j<=2*n-1; j++)
          { fj=xx[n*n*2+j];
            fr=fr+fj/(2.0*n);
            fg=fg+fj*fj;
          }
        fr=(fg-2.0*n*fr*fr)/(2.0*n-1.0);
        if (fr>=eps)
          { kk=kk+1;
            if (kk<k) it=1;
          }
      }
    for (i=0; i<=n-1; i++)
      { x[i]=0.0;
        for (j=0; j<=2*n-1; j++)
          x[i]=x[i]+xx[i*n*2+j]/(2.0*n);
      }
    z=jcplxf(x,n); x[n]=z;
    free(c); free(d); free(w);
    free(xt); free(xf);
    return(kk);
}
/////////////////////////////////////////////////////////////
#endif
/////////////////////////////////////////////////////////////
/*用户自编函数参考
void jmax1f(double x,double y[2])
{
	y[0]=(x-1.0)*(10.0-x);
    y[1]=-2.0*x+11.0;
    return;
}

double jmaxnf(double x[],int n,int j)
{ 
	double y;
    n=n;
    switch(j)
      { case 0: y=(x[0]-1.0)*(x[0]-10.0)
                  +(x[1]+2.0)*(x[1]+2.0)+2.0;
                break;
        case 1: y=2.0*(x[0]-1.0); break;
        case 2: y=2.0*(x[1]+2.0); break;
        defaut: {}
      }
    return(y);
}

double jjsimf(double x[],int n)
{ 
	double y;
    n=n;
    y=x[1]-x[0]*x[0]; y=100.0*y*y;
    y=y+(1.0-x[0])*(1.0-x[0]);
    return(y);
}

double jcplxf(double x[],int n)
{ 
	double y;
    n=n;
    y=-(9.0-(x[0]-3.0)*(x[0]-3.0));
    y=y*x[1]*x[1]*x[1]/(27.0*sqrt(3.0));
    return(y);
}
void jcplxs(int n,int m,double x[],double c[],double d[],double w[])
{
	n=n; m=m;
    c[0]=0.0; c[1]=0.0;
    d[0]=x[0]/sqrt(3.0); d[1]=6.0;
    w[0]=x[1]; w[1]=x[0]+x[1]*sqrt(3.0);
    return;
}
/**/