ch6.h

/// 科学与工程数值算法 （Java版）

/************************************************
 Expect bugs!
 Please use and enjoy, and let me know of any bugs/mods/improvements 
 that you have found/implemented and I will fix/incorporate them into 
 this file. Thank Mr. Xushiliang once again.

					hujinshan@2002.11.3
				Airforce Engineering University
************************************************/

/*****  #include "CH6.h"  数值积分*****/
#ifndef CH6_H_
#define CH6_H_

#include "stdlib.h"
#include "math.h"
#include "stdio.h"

#include "CH13.h"
//*******************************************************************
//用户自编函数，在具体应用中必须予以确定 (给函数指针(全局变量)赋值)，函数定义参见文件末 
//*******************************************************************
double fffts(double a,double b,double eps,double (*ffftsf)(double x));//变步长梯形求积法
double fsimp(double a,double b,double eps,double (*fsimpf)(double x));//变步长辛卜生求积法
double ffpts(double a,double b,double eps,double d,double (*ffptsf)(double x));//自适应梯形求积法
double fromb(double a,double b,double eps,double (*frombf)(double x));//龙贝格求积法
double ffpqg(double a,double b,double eps,double (*ffpqgf)(double x));//计算一维积分的连分式法
void fpart(double a,double b,int m,int n,double fa[],double fb[],double s[]);//高振荡函数求积法
double flrgs(double a,double b,double eps,double (*flrgsf)(double x));//勒让德-高斯求积法
double flags(double (*flagsf)(double x));//拉盖尔-高斯求积法
double fhmgs(double (*fhmgsf)(double x));//埃尔米特-高斯求积法
double fcbsv(double a,double b,double eps,double (*fcbsvf)(double x));//切比雪夫求积法
double fmtcl(double a,double b,double (*fmtclf)(double x));//计算一维积分的蒙特卡洛法
double fsim2(double a,double b,double eps,double (*fsim2f)(double x,double y),
		void (*fsim2s)(double x,double y[2]));//变步长辛卜生二重积分法
double fgaus(int n,int js[],double (*fgausf)(int n,double x[]),
		void (*fgauss)(int j,int n,double x[],double y[]));//计算多重积分的高斯方法
double fpqg2(double a,double b,double eps,double (*fpqg2f)(double x,double y),
		void (*fpqg2s)(double x,double y[]));//计算二重积分的连分式法
double fmtml(int n,double a[],double b[],double (*fmtmlf)(int n,double x[]));//计算多重积分的蒙特卡洛法

//*******************************************************************
double fffts(double a,double b,double eps,double (*ffftsf)(double x))
{ 
	//extern double ffftsf();
    int n,k;
    double fa,fb,h,t1,p,s,x,t;
    fa=ffftsf(a); fb=ffftsf(b);
    n=1; h=b-a;
    t1=h*(fa+fb)/2.0;
    p=eps+1.0;
    while (p>=eps)
      { s=0.0;
        for (k=0;k<=n-1;k++)
          { x=a+(k+0.5)*h;
            s=s+ffftsf(x);
          }
        t=(t1+h*s)/2.0;
        p=fabs(t1-t);
        t1=t; n=n+n; h=h/2.0;
      }
    return(t);
}
/////////////////////////////////////////////////////////////
double fsimp(double a,double b,double eps,double (*fsimpf)(double x))
{ 
	//extern double fsimpf();
    int n,k;
    double h,t1,t2,s1,s2,ep,p,x;
    n=1; h=b-a;
    t1=h*(fsimpf(a)+fsimpf(b))/2.0;
    s1=t1;
    ep=eps+1.0;
    while (ep>=eps)
      { p=0.0;
        for (k=0;k<=n-1;k++)
          { x=a+(k+0.5)*h;
            p=p+fsimpf(x);
          }
        t2=(t1+h*p)/2.0;
        s2=(4.0*t2-t1)/3.0;
        ep=fabs(s2-s1);
        t1=t2; s1=s2; n=n+n; h=h/2.0;
      }
    return(s2);
}
/////////////////////////////////////////////////////////////
static void ppp(double x0,double x1,double h,double f0,
				double f1,double t0,double eps,double d,double t[],
				double (*ffptsf)(double x))
{ 
	//extern double ffptsf();
    double x,f,t1,t2,p,g,eps1;
    x=x0+h/2.0; f=ffptsf(x);
    t1=h*(f0+f)/4.0; t2=h*(f+f1)/4.0;
    p=fabs(t0-(t1+t2));
    if ((p<eps)||(h/2.0<d))
      { t[0]=t[0]+(t1+t2); return;}
    else
      { g=h/2.0; eps1=eps/1.4;
        ppp(x0,x,g,f0,f,t1,eps1,d,t,ffptsf);
        ppp(x,x1,g,f,f1,t2,eps1,d,t,ffptsf);
        return;
      }
}
double ffpts(double a,double b,double eps,double d,double (*ffptsf)(double x))
{ 
	//extern double ffptsf();
    double h,t[2],f0,f1,t0,z;
    //void ppp();
    h=b-a; t[0]=0.0;
    f0=ffptsf(a); f1=ffptsf(b);
    t0=h*(f0+f1)/2.0;
    ppp(a,b,h,f0,f1,t0,eps,d,t,ffptsf);
    z=t[0]; return(z);
}
/////////////////////////////////////////////////////////////
double fromb(double a,double b,double eps,double (*frombf)(double x))
{ 
	//extern double frombf();
    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;
    h=b-a;
    y[0]=h*(frombf(a)+frombf(b))/2.0;
    m=1; n=1; ep=eps+1.0;
    while ((ep>=eps)&&(m<=9))
    {
		p=0.0;
        for (i=0;i<=n-1;i++)
        {
			x=a+(i+0.5)*h;
            p=p+frombf(x);
        }
        p=(y[0]+h*p)/2.0;
        s=1.0;
        for (k=1;k<=m;k++)
        {
			s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
        }
        ep=fabs(q-y[m-1]);
        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
	}
    return(q);
}
/////////////////////////////////////////////////////////////
double ffpqg(double a,double b,double eps,double (*ffpqgf)(double x))
{ 
	//extern double ffpqgf();
    int m,n,k,l,j;
    double h[10],bb[10],hh,t1,s1,ep,s,x,t2,g;
    m=1; n=1;
    hh=b-a; h[0]=hh;
    t1=hh*(ffpqgf(a)+ffpqgf(b))/2.0;
    s1=t1; bb[0]=s1; ep=1.0+eps;
    while ((ep>=eps)&&(m<=9))
      { s=0.0;
        for (k=0;k<=n-1;k++)
          { x=a+(k+0.5)*hh;
            s=s+ffpqgf(x);
          }
        t2=(t1+hh*s)/2.0;
        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;
        g=bb[m-1];
        for (j=m;j>=2;j--)
           g=bb[j-2]-h[j-2]/g;
        ep=fabs(g-s1);
        s1=g; t1=t2; hh=hh/2.0; n=n+n;
      }
    return(g);
}
/////////////////////////////////////////////////////////////
void fpart(double a,double b,int m,int n,double fa[],double fb[],double s[])
{ 
	int mm,k,j;
    double sa[4],sb[4],ca[4],cb[4],sma,smb,cma,cmb;
    sma=sin(m*a); smb=sin(m*b);
    cma=cos(m*a); cmb=cos(m*b);
    sa[0]=sma; sa[1]=cma; sa[2]=-sma; sa[3]=-cma;
    sb[0]=smb; sb[1]=cmb; sb[2]=-smb; sb[3]=-cmb;
    ca[0]=cma; ca[1]=-sma; ca[2]=-cma; ca[3]=sma;
    cb[0]=cmb; cb[1]=-smb; cb[2]=-cmb; cb[3]=smb;
    s[0]=0.0; s[1]=0.0;
    mm=1;
    for (k=0;k<=n-1;k++)
      { j=k;
        while (j>=4) j=j-4;
        mm=mm*m;
        s[0]=s[0]+(fb[k]*sb[j]-fa[k]*sa[j])/(1.0*mm);
        s[1]=s[1]+(fb[k]*cb[j]-fa[k]*ca[j])/(1.0*mm);
      }
    s[1]=-s[1];
    return;
}
/////////////////////////////////////////////////////////////
double flrgs(double a,double b,double eps,double (*flrgsf)(double x))
{
	//extern double flrgsf();
    int m,i,j;
    double s,p,ep,h,aa,bb,w,x,g;
    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};
    m=1;
    h=b-a; s=fabs(0.001*h);
    p=1.0e+35; ep=eps+1.0;
    while ((ep>=eps)&&(fabs(h)>s))
      { g=0.0;
        for (i=1;i<=m;i++)
          { aa=a+(i-1.0)*h; bb=a+i*h;
            w=0.0;
            for (j=0;j<=4;j++)
              { x=((bb-aa)*t[j]+(bb+aa))/2.0;
                w=w+flrgsf(x)*c[j];
              }
            g=g+w;
          }
        g=g*h/2.0;
        ep=fabs(g-p)/(1.0+fabs(g));
        p=g; m=m+1; h=(b-a)/m;
      }
    return(g);
}
/////////////////////////////////////////////////////////////
double flags(double (*flagsf)(double x))
{ 
	int i;
    double x,g;
    //extern double flagsf();
    static double t[5]={0.26355990,1.41340290,
            3.59642600,7.08580990,12.64080000};
    static double c[5]={0.6790941054,1.638487956,
             2.769426772,4.315944000,7.104896230};
    g=0.0;
    for (i=0; i<=4; i++)
      { x=t[i]; g=g+c[i]*flagsf(x); }
    return(g);
}
/////////////////////////////////////////////////////////////
double fhmgs(double (*fhmgsf)(double x))
{ 
	int i;
    double x,g;
    //extern double fhmgsf();
    static double t[5]={-2.02018200,-0.95857190,
                     0.0,0.95857190,2.02018200};
    static double c[5]={1.181469599,0.9865791417,
           0.9453089237,0.9865791417,1.181469599};
    g=0.0;
    for (i=0; i<=4; i++)
      { x=t[i]; g=g+c[i]*fhmgsf(x); }
    return(g);
}
/////////////////////////////////////////////////////////////
double fcbsv(double a,double b,double eps,double (*fcbsvf)(double x))
{
	//extern double fcbsvf();
    int m,i,j;
    double h,d,p,ep,g,aa,bb,s,x;
    static double t[5]={-0.8324975,-0.3745414,0.0,
                         0.3745414,0.8324975};
    m=1;
    h=b-a; d=fabs(0.001*h);
    p=1.0e+35; ep=1.0+eps;
    while ((ep>=eps)&&(fabs(h)>d))
      { g=0.0;
        for (i=1;i<=m;i++)
          { aa=a+(i-1.0)*h; bb=a+i*h;
            s=0.0;
            for (j=0;j<=4;j++)
              { x=((bb-aa)*t[j]+(bb+aa))/2.0;
                s=s+fcbsvf(x);
              }
            g=g+s;
          }
        g=g*h/5.0;
        ep=fabs(g-p)/(1.0+fabs(g));
        p=g; m=m+1; h=(b-a)/m;
      }
    return(g);
}
/////////////////////////////////////////////////////////////
double fmtcl(double a,double b,double (*fmtclf)(double x))
{ 
	int m;
    double r,d,x,s;
    //extern double fmtclf();
    //extern double mrnd1();//ch13.h
    r=1.0; s=0.0; d=10000.0;
    for (m=0; m<=9999; m++)
      { x=a+(b-a)*mrnd1(&r);
        s=s+fmtclf(x)/d;
      }
    s=s*(b-a);
    return(s);
}
/////////////////////////////////////////////////////////////
static double simp1(double x,double eps,
					double (*fsim2f)(double xf,double yf),
					void (*fsim2s)(double xs,double ys[2]))
{ 
	//extern void fsim2s();
    //extern double fsim2f();
    int n,i;
    double y[2],h,d,t1,yy,t2,g,ep,g0;
    n=1;
    fsim2s(x,y);
    h=0.5*(y[1]-y[0]);
    d=fabs(h*2.0e-06);
    t1=h*(fsim2f(x,y[0])+fsim2f(x,y[1]));
    ep=1.0+eps; g0=1.0e+35;
    while (((ep>=eps)&&(fabs(h)>d))||(n<16))
    { 
		yy=y[0]-h;
        t2=0.5*t1;
        for (i=1;i<=n;i++)
        { 
			yy=yy+2.0*h;
            t2=t2+h*fsim2f(x,yy);
        }
        g=(4.0*t2-t1)/3.0;
        ep=fabs(g-g0)/(1.0+fabs(g));
        n=n+n; g0=g; t1=t2; h=0.5*h;