ch6.h

数值处理算法程序

    }
    return(g);
}

double fsim2(double a,double b,double eps,double (*fsim2f)(double x,double y),
		void (*fsim2s)(double x,double y[2]))
{
	//double simp1();
    int n,j;
    double h,d,s1,s2,t1,x,t2,g,s,s0,ep;
    n=1; h=0.5*(b-a);
    d=fabs((b-a)*1.0e-06);
    s1=simp1(a,eps,fsim2f,fsim2s); s2=simp1(b,eps,fsim2f,fsim2s);
    t1=h*(s1+s2);
    s0=1.0e+35; ep=1.0+eps;
    while (((ep>=eps)&&(fabs(h)>d))||(n<16))
    {
		x=a-h; t2=0.5*t1;
		for (j=1;j<=n;j++)
        {
			x=x+2.0*h;
			g=simp1(x,eps,fsim2f,fsim2s);
			t2=t2+h*g;
		}
        s=(4.0*t2-t1)/3.0;
        ep=fabs(s-s0)/(1.0+fabs(s));
        n=n+n; s0=s; t1=t2; h=h*0.5;
    }
    return(s);
}
/////////////////////////////////////////////////////////////
double fgaus(int n,int js[],double (*fgausf)(int n,double x[]),
		void (*fgauss)(int j,int n,double x[],double y[]))
{ 
	//extern void fgauss();
    //extern double fgausf();
    int m,j,k,q,l,*is;
    double y[2],p,s,*x,*a,*b;
    static double t[5]={-0.9061798459,-0.5384693101,0.0,
                         0.5384693101,0.9061798459};
    static double c[5]={0.2369268851,0.4786286705,0.5688888889,
                        0.4786286705,0.2369268851};
    is=(int*)malloc(2*(n+1)*sizeof(int));
    x=(double*)malloc(n*sizeof(double));
    a=(double*)malloc(2*(n+1)*sizeof(double));
    b=(double*)malloc((n+1)*sizeof(double));
    m=1; l=1;
    a[n]=1.0; a[2*n+1]=1.0;
    while (l==1)//l==l 可以用一个非零常量代替，使条件永远为真
    { 
		for (j=m;j<=n;j++)
        { 
			fgauss(j-1,n,x,y);
            a[j-1]=0.5*(y[1]-y[0])/js[j-1];
            b[j-1]=a[j-1]+y[0];
            x[j-1]=a[j-1]*t[0]+b[j-1];
            a[n+j]=0.0;
            is[j-1]=1; is[n+j]=1;
		}
        j=n; q=1;
        while (q==1)
		{ 
			k=is[j-1];
            if (j==n) p=fgausf(n,x);
            else p=1.0;
            a[n+j]=a[n+j+1]*a[j]*p*c[k-1]+a[n+j];
            is[j-1]=is[j-1]+1;
            if(is[j-1]>5)
              if(is[n+j]>=js[j-1])
              { 
				  j=j-1; q=1;
                  if(j==0)
                  { 
					  s=a[n+1]*a[0];
					  free(is); free(x);free(a); free(b);
					  return(s);//这里是唯一的出口
                  }
              }
              else
              { 
				  is[n+j]=is[n+j]+1;
                  b[j-1]=b[j-1]+a[j-1]*2.0;
                  is[j-1]=1; k=is[j-1];
                  x[j-1]=a[j-1]*t[k-1]+b[j-1];
                  if(j==n)  q=1;
                  else  q=0;
              }
			else
            { 
				k=is[j-1];
                x[j-1]=a[j-1]*t[k-1]+b[j-1];
                if (j==n) q=1;
                else q=0;
			}
		}
        m=j+1;
     }
//	return 0;//此行仅仅为了消除警告: not all control paths return a value
}//VC可能给出一个警告，不要紧，出口是唯一的，也可以消除它，见上
/////////////////////////////////////////////////////////////
static double pqg1(double x,double eps,double (*fpqg2f)(double x,double y),
				   void (*fpqg2s)(double xs,double ys[]))
{
	//extern void fpqg2s();
    //extern double fpqg2f();
    int m,n,k,l,j;
    double b[10],h[10],y[2],hh,t1,t2,s0,yy,g,ep,s;
    m=1; n=1;
    fpqg2s(x,y);
    hh=y[1]-y[0]; h[0]=hh;
    t1=0.5*hh*(fpqg2f(x,y[0])+fpqg2f(x,y[1]));
    s0=t1; b[0]=t1; ep=1.0+eps;
    while ((ep>=eps)&&(m<=9))
      { t2=0.5*t1;
        for (k=0;k<=n-1;k++)
          { yy=y[0]+(k+0.5)*hh;
            t2=t2+0.5*hh*fpqg2f(x,yy);
          }
        m=m+1;
        h[m-1]=h[m-2]/2.0;
        g=t2; l=0; j=2;
        while ((l==0)&&(j<=m))
          { s=g-b[j-2];
            if (fabs(s)+1.0==1.0) l=1;
            else g=(h[m-1]-h[j-2])/s;
            j=j+1;
          }
        b[m-1]=g;
        if (l!=0) b[m-1]=1.0e+35;
        s=b[m-1];
        for (j=m;j>=2;j--) s=b[j-2]-h[j-2]/s;
        ep=fabs(s-s0)/(1.0+fabs(s));
        n=n+n; t1=t2; s0=s; hh=0.5*hh;
      }
    return(s);
}
double fpqg2(double a,double b,double eps,double (*fpqg2f)(double x,double y),
		void (*fpqg2s)(double x,double y[]))
{
	int m,n,k,l,j;
    //double pqg1();
    double bb[10],h[10],hh,s1,s2,t1,t2,x,g,s0,ep,s;
    m=1; n=1;
    hh=b-a; h[0]=hh;
    s1=pqg1(a,eps,fpqg2f,fpqg2s); s2=pqg1(b,eps,fpqg2f,fpqg2s);
    t1=hh*(s1+s2)/2.0;
    s0=t1; bb[0]=t1; ep=1.0+eps;
    while ((ep>=eps)&&(m<=9))
	{
		t2=0.5*t1;
        for (k=0;k<=n-1;k++)
		{ 
			x=a+(k+0.5)*hh;
            s1=pqg1(x,eps,fpqg2f,fpqg2s);
            t2=t2+0.5*s1*hh;
		}
        m=m+1;
        h[m-1]=h[m-2]/2.0;
        g=t2; l=0; j=2;
        while ((l==0)&&(j<=m))
		{ 
			s=g-bb[j-2];
            if (fabs(s)+1.0==1.0) l=1;
            else g=(h[m-1]-h[j-2])/s;
            j=j+1;
		}
        bb[m-1]=g;
        if (l!=0) bb[m-1]=1.0e+35;
        s=bb[m-1];
        for(j=m;j>=2;j--) 
			s=bb[j-2]-h[j-2]/s;
        ep=fabs(s-s0)/(1.0+fabs(s));
        n=n+n; t1=t2; s0=s; hh=hh/2.0;
      }
    return(s);
}
/////////////////////////////////////////////////////////////
double fmtml(int n,double a[],double b[],double (*fmtmlf)(int n,double x[]))
{ 
	int m,i;
    double r,s,d,*x;
    //extern double mrnd1();//ch13.h
    //extern double fmtmlf();
    x=(double*)malloc(n*sizeof(double));
    r=1.0; d=10000.0; s=0.0;
    for (m=0; m<=9999; m++)
    {
		for (i=0; i<=n-1; i++)
          x[i]=a[i]+(b[i]-a[i])*mrnd1(&r);
        s=s+fmtmlf(n,x)/d;
    }
    for (i=0; i<=n-1; i++)
      s=s*(b[i]-a[i]);
    free(x); return(s);
}
#endif
/////////////////////////////////////////////////////////////
/*用户自己编的函数定义参考:
double ffftsf(double x)
{
	double y;
    y=exp(-x*x);
    return(y);
}

double fsimpf(double x)
{ 
	double y;
    y=log(1.0+x)/(1.0+x*x);
    return(y);
}
  
double ffptsf(double x)
{ 
	double y;
    y=1.0/(1.0+25.0*x*x);
    return(y);
}

double frombf(double x)
{ 
	double y;
    y=x/(4.0+x*x);
    return(y);
}

double ffpqgf(double x)
{
	double y;
    y=exp(-x*x);
    return(y);
}

double flrgsf(double x)
{
	double y;
    y=x*x+sin(x);
    return(y);
}

double flagsf(double x)
{ 
	double y;
    y=x*exp(-x);
    return(y);
}

double fhmgsf(double x)
{ 
	double y;
    y=x*x*exp(-x*x);
    return(y);
}

double fcbsvf(double x)
{ 
	double y;
    y=x*x+sin(x);
    return(y);
}

double fmtclf(double x)
{ 
	double y;
    y=x*x+sin(x);
    return(y);
}

void fsim2s(double x,double y[2])//double x,y[2];
{
	y[0]=-sqrt(1.0-x*x);
    y[1]=-y[0];
    return;
}
double fsim2f(double x,double y)
{
	double z;
    z=exp(x*x+y*y);
    return(z);
}

void fgauss(int j,int n,double x[],double y[])
{ 
	double q;
    n=n;
    switch (j)
    { 
	case 0: { y[0]=0.0; y[1]=1.0; break;}
	case 1: 
		{ 
			y[0]=0.0; y[1]=sqrt(1.0-x[0]*x[0]); 
			break;
		}
	case 2: 
		{ 
			q=x[0]*x[0]+x[1]*x[1]; y[0]=sqrt(q);
			y[1]=sqrt(2.0-q); 
			break;
		}
    default: { }
	}
    return;
}

double fgausf(int n,double x[])
{ 
	double z;
    n=n;
    z=x[2]*x[2];
    return(z);
}

void fpqg2s(double x,double y[])
{ 
	y[1]=sqrt(1.0-x*x);
    y[0]=-y[1];
    return;
}

double fpqg2f(double x,double y)
{ 
	double z;
    z=exp(x*x+y*y);
    return(z);
}

double fmtmlf(int n,double x[])
{
	int i;
    double f;
    f=0.0;
    for (i=0; i<=n-1; i++)
      f=f+x[i]*x[i];
    return(f);
}
/**/