home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMERICA.PAK / GAUSSIAN.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  2.9 KB  |  93 lines
  1.  
  2. (* :Name: NumericalMath`GaussianQuadrature` *)
  3.  
  4. (* :Title: Coefficients for Gaussian Quadrature *)
  5.  
  6. (* :Author: Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     Calculates the weights and abscissae for the elementary Gaussian
  10.     quadrature rule with n points on the interval a to b.  Calculates
  11.     the error in the elementary Gaussian quadrature formula with n
  12.     points on the interval a to b.
  13. *)
  14.  
  15. (* :Context: NumericalMath`GaussianQuadrature` *)
  16.  
  17. (* :Package Version: Mathematica 2.0 *)
  18.  
  19. (* :Copyright: Copyright 1990,  Wolfram Research, Inc.
  20.     Permission is hereby granted to modify and/or make copies of
  21.     this file for any purpose other than direct profit, or as part
  22.     of a commercial product, provided this copyright notice is left
  23.     intact.  Sale, other than for the cost of media, is prohibited.
  24.  
  25.     Permission is hereby granted to reproduce part or all of
  26.     this file, provided that the source is acknowledged.
  27. *)
  28.  
  29. (* :History:
  30.     Revised by Jerry B. Keiper, December, 1990.
  31.     Original version by Jerry B. Keiper, May, 1990.
  32. *)
  33.  
  34. (* :Keywords: Gaussian quadrature, abscissae, weights *)
  35.  
  36. (* :Source: Any elementary numerical analysis textbook. *)
  37.  
  38. (* :Mathematica Version: 2.0 *)
  39.  
  40. (* :Limitation: *)
  41.  
  42. BeginPackage["NumericalMath`GaussianQuadrature`"]
  43.  
  44. GaussianQuadratureWeights::usage =
  45. "GaussianQuadratureWeights[n, a, b, prec] gives a list of the pairs
  46. {abscissa, weight} to prec digits precision for the elementary n-point
  47. Gaussian quadrature formula for quadrature on the interval a to b.  The
  48. argument prec is optional."
  49.  
  50. GaussianQuadratureError::usage =
  51. "GaussianQuadratureError[n, f, a, b, prec] gives the leading term in the
  52. error in the elementary n-point Gaussian quadrature formula to prec digits
  53. precision for the function f on an interval from a to b.  The argument prec
  54. is optional."
  55.  
  56. Begin["NumericalMath`GaussianQuadrature`Private`"]
  57.  
  58. GaussianQuadratureWeights[n_Integer, a_, b_, prec_Integer:$MachinePrecision] :=
  59.     GQW[n, a, b, prec] /; (n > 0 && prec > 0)
  60.  
  61. GaussianQuadratureError[n_Integer, f_, a_, b_, prec_Integer:$MachinePrecision] :=
  62.     GQE[n, f, a, b, prec] /; (n > 0 && prec > 0)
  63.  
  64.  
  65. GQW[n_, a_, b_, prec_] :=
  66.     Module[{i, w, x, m, rhs, legn, t},
  67.         legn = Numerator[LegendreP[n,t]];
  68.         If[ n == 1, Return[N[{{0, 2}}, prec]]];
  69.         x = Apply[List, Map[#[[2]]&, NRoots[legn == 0, t, prec]]];
  70.         m = Table[x^i, {i,n-1}];
  71.         PrependTo[m, Table[1, {n}]];
  72.         rhs = Table[(1+Sign[(-1)^i])/(i+1), {i,0,n-1}];
  73.         w = Expand[LinearSolve[m, rhs](b-a)/2];
  74.         x = Expand[(x (b-a) + (b+a))/2];
  75.         Transpose[{x, w}]
  76.         ]
  77.  
  78. GQE[n_, f_, a_, b_, prec_] :=
  79.     Module[{t, fint, fs, x, w, len},
  80.     fs = t^(2n)/(2n)!;
  81.     {x, w} = Transpose[GQW[n, -len, len, prec]];
  82.     fs = Expand[fs /. t->x] . w;
  83.     fint = 2 len^(2n+1)/(2n+1)!;
  84.     (Expand[fs - fint] /. len -> (b-a)/2) Derivative[2n][f]
  85.     ]
  86.  
  87. End[] (* NumericalMath`GaussianQuadrature`Private` *)
  88.  
  89. Protect[GaussianQuadratureWeights, GaussianQuadratureError];
  90.  
  91. EndPackage[] (* NumericalMath`GaussianQuadrature` *)
  92.  
  93.