home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMERICA.PAK / LISTINTE.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  4.5 KB  |  150 lines
  1.  
  2. (* :Name: NumericalMath`ListIntegrate` *)
  3.  
  4. (* :Title:  Approximate Integration of Functions That Are Known Only at a Few
  5.         Distinct Points *)
  6.  
  7. (* :Author: Jerry B. Keiper *)
  8.  
  9. (* :Summary:
  10.     Given a function in the form of a sorted list of ordered pairs
  11.     (or simply as a list of ordinates, assuming equally spaced
  12.     abscissae), ListIntegrate gives an approximation to the definite
  13.     integral of that function from the first to the last abscissa.
  14. *)
  15.  
  16. (*:Context: NumericalMath`ListIntegrate` *)
  17.  
  18. (*:Package Version: Mathematica 2.0 *)
  19.  
  20. (* :Copyright: Copyright 1990  Wolfram Research, Inc.
  21.     Permission is hereby granted to modify and/or make copies of
  22.     this file for any purpose other than direct profit, or as part
  23.     of a commercial product, provided this copyright notice is left
  24.     intact.  Sale, other than for the cost of media, is prohibited.
  25.  
  26.     Permission is hereby granted to reproduce part or all of
  27.     this file, provided that the source is acknowledged.
  28. *)
  29.  
  30. (* :History:
  31.     Revised by Jerry B. Keiper, November, 1990.
  32.     Original version by Jerry B. Keiper, March, 1989.
  33. *)
  34.  
  35. (* :Keywords: integration, numerical integration, interpolation *)
  36.  
  37. (* :Source: any elementary numerical analysis textbook. *)
  38.  
  39. (* :Warnings:
  40.     ListIntegrate[ ] should not be used if the data contain
  41.     significant error.  In such cases the error should be smoothed
  42.     using Fit[ ], and the fitted function should be integrated
  43.     using Integrate[ ] or NIntegrate[ ].
  44. *)
  45.  
  46. (* :Limitations:
  47.     The list {{x0,y0},{x1,y1},...,{xn,yn}} is assumed to be in sorted
  48.     order (i.e., the xi's are assumed to be monotonically increasing);
  49.     no checking is done.  It does not allow integration over an
  50.     arbitrary interval; you must integrate from the initial point to
  51.     the final point.
  52. *)
  53.  
  54. BeginPackage["NumericalMath`ListIntegrate`"]
  55.  
  56. ListIntegrate::usage =
  57. "ListIntegrate[{y0,y1,...,yn}, h, k] gives an approximation to the integral of
  58. a function from x = x0 to x = xn with constant step size h between the
  59. consecutive x values.  The function integrated is piecewise of degree k-1 and
  60. is the interpolating polynomial through the k/2 points to the left and k/2
  61. points to the right of x.  If k is odd it is increased by 1.  k must be less
  62. than or equal to the total number of data points.  The default value of k is 4.
  63. ListIntegrate[{{x0,y0},{x1,y1},...,{xn,yn}},k] can be used for variable
  64. stepsize data.  If the data are known to contain errors, you may be better off
  65. performing Integrate[] on the result of Fit[] applied to the data."
  66.  
  67.  
  68. Begin["NumericalMath`ListIntegrate`Private`"]
  69.  
  70. ListIntegrate[cl_?VectorQ, h_, k_Integer:4] :=
  71.     Module[{ans},
  72.         ans /; (ans = ListIntegrate0[h, cl, k]) =!= $Failed
  73.     ] /; k > 0;
  74.  
  75. ListIntegrate[cl_?MatrixQ, k_Integer:4] :=
  76.     Module[{ans, dim = Length[Dimensions[cl]]},
  77.         ans /; ((dim == 2) &&
  78.             (Length[cl[[1]]] == 2) &&
  79.             ((ans = ListIntegrate1[cl, k]) =!= $Failed))
  80.     ] /; k > 0;
  81.  
  82.  
  83. ListIntegrate::degerr = "Number of data points (`1`) is insufficient for the
  84. degree of interpolant requested (`2`)."
  85.  
  86. ListIntegrate0[h_, cl_, k_] := 
  87.     Module[{ke = k + If[IntegerQ[k/2], 0, 1], data},
  88.         If[Length[cl] < ke,
  89.             Message[ListIntegrate::degerr,Length[cl],ke-1];
  90.             Return[$Failed]
  91.             ];
  92.         If[ke == 2,
  93.             h (Apply[Plus,cl] - (cl[[1]] + cl[[-1]])/2),
  94.             (* else *)
  95.             data = Array[h #&,{Length[cl]}];
  96.             data = Transpose[{data,cl}];
  97.             ListInt[data,ke]
  98.             ]
  99.         ];
  100.  
  101. ListIntegrate1[cl_, k_] := 
  102.     Module[{ke = k + If[IntegerQ[k/2], 0, 1], x, y},
  103.         If[Length[cl] < ke,
  104.             Message[ListIntegrate::degerr,Length[cl],ke-1];
  105.             Return[$Failed]
  106.             ];
  107.         If[ke == 2,
  108.             y = Transpose[cl];
  109.             x = Drop[y[[1]],1]-Drop[y[[1]],-1];
  110.             y = Drop[y[[2]],1]+Drop[y[[2]],-1];
  111.             (x . y)/2,
  112.             (* else *)
  113.             ListInt[cl,ke]
  114.             ]
  115.         ];
  116.  
  117. ListInt[cl_,k_] :=
  118.     Module[{i,tmp,sum,x,m=k/2,m1=k/2+1},
  119.         If[Length[cl] == k,
  120.             IntInt[cl,1,-1],
  121.             (* else *)
  122.             sum = IntInt[Take[cl,k],1,m1];
  123.             Do[sum += IntInt[Take[cl,{i,i+k-1}],m,m1],
  124.                 {i,2,Length[cl]-k}];
  125.             tmp = Take[cl,-k];
  126.             sum + IntInt[Take[cl,-k],m,k]
  127.             ]
  128.         ];
  129.  
  130. IntInt[data_,i_,j_] :=
  131.     Module[{cl, x, k, sum, tmp},
  132.         cl = CoefficientList[InterpolatingPolynomial[data,x],x];
  133.         x = data[[j,1]];
  134.         tmp = 0;
  135.         Do[tmp = (tmp + cl[[k]]/k) x,{k,Length[cl],1,-1}];
  136.         sum = tmp;
  137.         x = data[[i,1]];
  138.         tmp = 0;
  139.         Do[tmp = (tmp + cl[[k]]/k) x,{k,Length[cl],1,-1}];
  140.         sum - tmp
  141.         ];
  142.  
  143. End[] (* NumericalMath`ListIntegrate`Private` *)
  144.  
  145. Protect[ListIntegrate];
  146.  
  147. EndPackage[] (* NumericalMath`ListIntegrate` *)
  148.  
  149.  
  150.