单精度simpson.rute法求二维有限区间的积分.cpp

4个关于二维积分函数的C程序

#include <stdio.h>
#include <math.h>

float func(float x1,float x2);

       ///////////////////////////////////////////////////////////////
       // 单精度Simpson.Rute求二维有限区域的积分                    //
       // a,b第一变量的积分上下限;c,d第二变量的积分上下限           //
       // m1(1<=m1<=500)区间(a,b)的等分数;n1>=1区间(c,d)的等分数;   //
       // func(x1,x2)被积函数                                       //
       // 调用方式: rtn=sekbnd(func,a,b,c,d,m1,n1,&result,flag)     //
       //    rtn=0 :正常返回;                                       //
       //       =-1:程序异常:m1不满足1<=m1<=500;flag=1打印信息      //
       //       =-2:程序异常:n1<0;flag=1打印信息                    //
       ///////////////////////////////////////////////////////////////
int sekbnd(float (*func)(float,float ),float a,float b,float c,float d,int m1,int n1,float *result,int flag)
  { double swk,dr;
    float s[1001],dx,dy,x,y;
    int m,n,mm,i,j,error=0;
    
    if (m1<1||m1>500)
       { error=-1;
         if (flag)
           printf("\nM1=%d,M1 must be smaller than 500 and larger than 1!\n",m1);
         return(error);
         }
    if (n1<1)
       { error=-2;
         if (flag)
           printf("\nN1=%d,N1 must be  larger than 0!\n",n1);
         return(error);
         }
    m=m1<<1;   n=n1<<1;
    mm=m+1;
    dx=(b-a)/m;
    dy=(d-c)/n;
    for (i=1;i<=mm;i++)
       { x=a+(i-1)*dx;
         swk=(*func)(x,c)-(*func)(x,d);
         s[i-1]=swk;
         for (j=1;j<=n;j+=2)
            {  y=c+j*dy;
               swk=swk+4*(*func)(x,y)+2*(*func)(x,y+dy);
               s[i-1]=swk;
               }
         }
    dr=s[0]-s[mm-1];
    for (i=2;i<=m;i+=2)
      dr=dr+4*s[i-1]+2*s[i];
    *result=dr*dx*dy/9;
    return(error);
    }

   void main(  )
   {  float a,b,c,d;
      int m1,n1,flag,rtn;
      float result;

      a=c=0;
      b=d=1;
      flag=0;
      m1=n1=5;
      rtn=sekbnd(func,a,b,c,d,m1,n1,&result,flag);
      if (!rtn)
        printf("\nRESULT=%f\n",result);
      else
        printf("\nIt's wrong!Return code is %d!\n",rtn);

      scanf("%d",&m1);
      }

   float func(float x1,float x2)
   { float y;
     y=exp(-x1-x2);
     return(y);
     }