| 网站首页 | 业界新闻 | 小组 | 威客 | 人才 | 下载频道 | 博客 | 代码贴 | 在线编程 | 编程论坛
欢迎加入我们,一同切磋技术
用户名:   
 
密 码:  
共有 1340 人关注过本帖
标题:够难!帮忙处理一下floating point error :domain(急用)
取消只看楼主 加入收藏
木美子
Rank: 1
等 级:新手上路
帖 子:1
专家分:0
注 册:2008-5-25
收藏
 问题点数:0 回复次数:0 
够难!帮忙处理一下floating point error :domain(急用)
出错处已经标出
#include "Stdio.h"
#include "Conio.h"
#include "stdlib.h"
#define pi 3.1415
#define kkg 9
double r0=1;
double Dealpos1(int i)
{double a[5]={ 14.5,15.8,16,15,12.5};
 return a[i];
}
double Dealpos2(int i)
{double b[5]={1.8,6,9.5,14.5,19};
  return b[i];
}
double objf(double x[])
{float f=0;
 float Lmin=1.772,dlt_L=0.25;
 float c1,c2,c3,L,d1,d2,aa,bb;
 int n;
 for(n=0;n<5;n++)
    { L=Lmin+dlt_L*n;
      c1=acos((x[0]*x[0]+2.5*2.5-L*L)/(2*2.5*x[0]));
      c2=acos(x[5]-x[0]*cos(x[1]))/sqrt(x[0]*x[0]+x[5]*x[5]-2*x[0]*x[5]*cos(x[1]));
      c3=acos((x[2]*x[2]+x[3]*x[3]-2.385*2.385)/(2*x[2]*x[3]));
      d1=c1-x[1]-pi/6;
      d2=d1-(c2+c3);
      aa=x[5]*cos(d1)+x[4]*cos(d2);
      bb=x[5]*sin(d1)+x[4]*sin(d2)+9;
      f=f+(aa-Dealpos1(n))*(bb-Dealpos2(n));
    }
    return(f);
 }
 void strain(double x[],double c[])
 {float Lmin=1.772,dlt_L=0.25;
  float L;
   int n;
   for (n=0;n<5;n++)
    {L=Lmin+dlt_L*n;
      c[1]=pi-acos((x[0]*x[0]+2.5*2.5-L*L)/(2*2.5*x[0]));
      c[2]=acos((x[0]*x[0]+2.5*2.5-L*L)/(2*2.5*x[0]));
      c[3]=-acos((x[0]*x[0]+2.5*2.5-L*L)/(2*2.5*x[0]))+acos((x[0]*x[0]+x[5]*x[5]-x[1]*x[1])/(2*x[0]*x[5]))+2*pi/3;
      c[4]=acos((x[0]*x[0]+2.5*2.5-L*L)/(2*2.5*x[0]))+acos((x[0]*x[0]+x[5]*x[5]-x[1]*x[1])/(2*x[0]*x[5]))+pi/3;
      c[5]=pi-acos((x[2]*x[2]+x[3]*x[3]-2.385*2.385)/(2*x[2]*x[3]));
      c[6]=acos((x[2]*x[2]+x[3]*x[3]-2.385*2.385)/(2*x[2]*x[3]));
      c[7]=pi/2-acos((x[1]*x[1]+x[5]*x[5]-x[0]*x[0])/(2*x[1]*x[5]));
      c[8]=acos((x[1]*x[1]+x[5]*x[5]-x[0]*x[0])/(2*x[2]*x[5]));
      c[9]=x[5]-1-x[4];
    }
 }
 double ldf(double *x)
 {int i;
  double f,sg;
  double c[kkg];
  sg=0;
  strain(x,c);
  for(i=0;i<kkg;i++)
   {if(c[i]>0)
    sg=sg+r0/c[i];
    else
    sg=sg-c[i]*1e5;
   }
   f=objf(x)+sg;
   return(f);
 }
 void ii(double *p,double a[],double b[],double s[],double h0,int n)
 {int i;
  double *x[3],h,f1,f2,f3;
  for(i=0;i<3;i++)
   x[i]=(double *)malloc(n*sizeof(double));
  h=h0;
  for(i=0;i<n;i++)
   *(x[0]+i)=*(p+i);
  f1=ldf(x[0]);
  for(i=0;i<n;i++)
   *(x[1]+i)=*(x[0]+i)+h*s[i];
  f2=ldf(x[1]);
  if(f2>=f1)
   {h=-h0;
    for(i=0;i<n;i++)
    *(x[2]+i)=*(x[0]+i);
    f3=f1;
    for(i=0;i<n;i++)
    {*(x[0]+i)=*(x[1]+i);
     *(x[1]+i)=*(x[2]+i);
    }
    f1=f2;
    f2=f3;
   }
  for(;;)
   {h=2.*h;
    for(i=0;i<n;i++)
    *(x[2]+i)=*(x[1]+i)+h*s[i];
    f3=ldf(x[2]);
    if(f2<f3)
     break;
    else
     {for(i=0;i<n;i++)
       {*(x[0]+i)=*(x[1]+i);
        *(x[1]+i)=*(x[2]+i);
       }
       f1=f2;
       f2=f3;
     }
 }
  if(h<0.)
   for(i=0;i<n;i++)
    {a[i]=*(x[2]+i);
     b[i]=*(x[0]+i);
    }
  else
   for(i=0;i<n;i++)
    {a[i]=*(x[0]+i);
     b[i]=*(x[2]+i);
    }
  for(i=0;i<3;i++)
   free(x[i]);
 }
 double gold(double a[],double b[],double xx[],int n,double eps1)
 {int i;
  double f1,f2,*x[2],ff,q,w;
  for(i=0;i<2;i++)
  x[i]=(double *)malloc(n*sizeof(double));
  for(i=0;i<n;i++)
   {*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
    *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
   }
  f1=ldf(x[0]);
  f2=ldf(x[1]);
 do
  {if(f1>f2)
   {for(i=0;i<n;i++)
     {b[i]=*(x[0]+i);
      *(x[0]+i)=*(x[1]+i);
      }
     f1=f2;
     for(i=0;i<n;i++)
      *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
      f2=ldf(x[1]);
    }
   {for(i=0;i<n;i++)
    {a[i]=*(x[1]+i);
     *(x[1]+i)=*(x[0]+i);
     }
    f2=f1;
    for(i=0;i<n;i++)
     *(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
    f1=ldf(x[0]);
   }
   q=0;
   for(i=0;i<n;i++)
   q=q+(b[i]-a[i])*(b[i]-a[i]);
   w=sqrt(q);
  }while(w>eps1);
 for(i=0;i<n;i++)
  xx[i]=0.5*(a[i]+b[i]);
 ff=ldf(xx);
 for(i=0;i<2;i++)
  free(x[i]);
 return(ff);
}
double powell(double *p,double x[],double epss,int n)
 {int i,j,m;
  double *xx[6],*ss,*s,f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
  s=(double*)malloc(n*sizeof(double));
  ss=(double*)malloc(n*(n+1)*sizeof(double));
  for(i=0;i<n;i++)
   {for(j=0;j<(n+1);j++)
   *(ss+i*(n+1)+j)=0;
   *(ss+i*(n+1)+i)=1;
   }
  for(i=0;i<6;i++)
   xx[i]=(double*)malloc(n*sizeof(double));
  for(i=0;i<n;i++)
   *(xx[0]+i)=*(p+i);
  for(;;)
   {for(i=0;i<n;i++)
    {*(xx[1]+i)=*(xx[0]+i);
     x[i]=*(xx[1]+i);
    }
    f0=f1=ldf(x);
    dlt=-1;
    for(j=0;j<n;j++)
    {for(i=0;i<n;i++)
     {*(xx[1]+i)=x[i];
      *(s+i)=*(ss+i*(n+1)+j);
      }
     ii(xx[0],xx[4],xx[5],s,0.1,n);
     f=gold(xx[4],xx[5],x,n,1e-5);
     df=f0-f;
     if(df>dlt)
     {dlt=df;
      m=j;
      }
     }
    }
  sdx=0.;
  for(i=0;i<n;i++)
   sdx=sdx+fabs(x[i]-*(xx[1]+i));
  if(sdx<epss)
   {for(i=0;i<6;i++)
    free(xx[i]);
    free(s);
    free(ss);
    return(f);
   }
  for(i=0;i<n;i++)
   *(xx[2]+i)=x[i];
  f2=f;
  for(i=0;i<n;i++)
   {*(xx[3]+i)=2.*(*xx[2]+i)-*(xx[1]+i);
    x[i]=*(xx[3]+i);
   }
  fx=ldf(x);
  f3=fx;
  q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
  d=0.5*dlt*(f1-f3)*(f1-f3);
  if((f3<f1)||(q<d))
    {if(f2<=f3)
      for(i=0;i<n;i++)
      *(xx[0]+i)=*(xx[2]+i);
     else
      for(i=0;i<n;i++)
       *(xx[0]+i)=*(xx[3]+i);
    }
  else
    {for(i=0;i<n;i++)
      {*(ss+i*(n+1)+n+1)=x[i]-*(xx[1]+i);
      *(s+i)=*(ss+i*(n+1)+n+1);
      }
     ii(xx[0],xx[4],xx[5],s,0.1,n);
     f=gold(xx[4],xx[5],x,n,1e-5);
     for(i=0;i<n;i++)
      *(xx[0]+i)=x[i];
     for(j=m+1;j<(n+1);j++)
      for(i=0;i<n;i++)
       *(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
    }
   }
 double empf(double *p,double x[],double c,double eps,int n)
{int i,Tm;
 double fom,fxo;
 fom=1e9;
 Tm=0;
 do{fxo=powell(p,x,0.01,n);
    if(fabs(fom-fxo)>eps)
     {fom=fxo;
      r0=c*r0;
      Tm=Tm+1;
      printf("Now it is the %d time of iteration,current values of variables are:\n",Tm);
      for(i=0;i<n;i++)
       {*(p+i)=x[i];
        printf("x[%d}=%f\n",i,x[i]);
       }
      printf("The value of objective function is:%f\n",fxo);
      printf("****************\n");
     }
    else
     return(fxo);
   }while(1);
}
void main()
 {
 int i;
 double z[]={2.4,7.5,3.3,1.0,6.5,9.5};
 double c,x[6];
 c=empf(z,x,0.2,1e-5,6);
 printf("\n Output the optimum results:\n");
 for(i=0;i<6;i++)
  printf("%f\n",x[i]);
  printf("The minimum objiective value is %f\n",c);
 }

[[it] 本帖最后由 木美子 于 2008-5-26 08:17 编辑 [/it]]
搜索更多相关主题的帖子: double domain point floating 
2008-05-25 23:17
快速回复:够难!帮忙处理一下floating point error :domain(急用)
数据加载中...
 
   



关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.013056 second(s), 8 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved