share/numeric/dblint.mac

   1 define_variable(dblint_y,10,fixnum,"number of y divisions");
   2 define_variable(dblint_x,10,fixnum,"number of x divisions");
   3 dblint(f,c,d,a,b):=
   4      (mode_declare([function(f,c,d),a,b],float),
   5                   block([m2,n2,h,j1,j2,j3,x,dox,cox,hx,k1,k2,k3,y,z,l],
   6                         mode_declare([m2,n2,h,j1,j2,j3],float),
   7                                 n2:.5/dblint_x,
   8         m2:.5/dblint_y,
   9                                 h:(b-a)*n2,
  10                                 j1:0.,j2:0.,j3:0.,
  11                                 for i:0 thru 2*dblint_x step 1 do
  12                                 (mode_declare(i,fixnum,[x,dox,cox,hx,k1,k2,k3],float),
  13                                              x:a+float(i)*h,
  14                                              dox:apply(d,[x]),cox:apply(c,[x]),
  15                                              hx:(dox-cox)*m2,
  16                                              k1:apply(f,[x,cox])+
  17                                                 apply(f,[x,dox]),
  18                                              k2:0.,k3:0.,
  19                                     for j:1 thru 2*dblint_y-1 step 1 do
  20                                     (mode_declare(j,fixnum,[y,z,l],float),
  21                                              y:cox+float(j)*hx,
  22                                              z:apply(f,[x,y]),
  23                                              if evenp(j) then k2:k2+z else k3:k3+z),
  24                                     l:(k1+2.*k2+4.*k3)*hx/3.,
  25                                     if (i=0 or i=2*dblint_x) then j1:j1+l else
  26                                        if evenp(i) then j2:j2+l else j3:j3+l),
  27                                 return((j1+2.*j2+4.*j3)*h/3.)));