home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMERICA.PAK / CAUCHYPR.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  4.2 KB  |  138 lines
  1.  
  2. (*:Name: NumericalMath`CauchyPrincipalValue` *)
  3.  
  4. (*:Title: Numerical Evaluation of Cauchy Principal Values *)
  5.  
  6. (*:Author: Jerry B. Keiper *)
  7.  
  8. (*:Summary: Numerically evaluates Cauchy principal values of integrals that
  9.     may not be Riemann integrable.  Uses NIntegrate[ ] in a symmetric
  10.     way to get cancelation.
  11. *)
  12.  
  13. (*:Context: NumericalMath`CauchyPrincipalValue` *)
  14.  
  15. (*:Package Version: Mathematica 2.0 *)
  16.  
  17. (* :Copyright: Copyright 1990,  Wolfram Research, Inc.
  18.     Permission is hereby granted to modify and/or make copies of
  19.     this file for any purpose other than direct profit, or as part
  20.     of a commercial product, provided this copyright notice is left
  21.     intact.  Sale, other than for the cost of media, is prohibited.
  22.  
  23.     Permission is hereby granted to reproduce part or all of
  24.     this file, provided that the source is acknowledged.
  25. *)
  26.  
  27. (* :History:
  28.     Revised by Jerry B. Keiper, November, 1990.
  29.     Version 1.2 by Jerry B. Keiper, January, 1990.
  30. *)
  31.  
  32. (*:Keywords: Cauchy principal value, integration, divergent integral *)
  33.  
  34. (*:Source:
  35.     E. T. Whittaker & G. N. Watson, A Course of Modern Analysis,
  36.     Cambridge University Press, (4th ed.), section 4.5
  37. *)
  38.  
  39. (*:Mathematica Version: 2.0 *)
  40.  
  41. (*:Limitation:
  42.     Since CauchyPrincipalValue is a numerical function and it
  43.     relies on symmetry near the singular points, loss of digits
  44.     due to cancelation can be a problem.  The user must know
  45.     and specify EXACTLY where the singularities are, or the
  46.     result can be grossly in error.  Again this is because of
  47.     the extreme sensitivity to perturbations.  
  48. *)
  49.  
  50. BeginPackage["NumericalMath`CauchyPrincipalValue`"]
  51.  
  52. CauchyPrincipalValue::usage = 
  53. "CauchyPrincipalValue[f[x], {x, a, {b, eps}, c}, options] gives the Cauchy
  54. principal value of NIntegrate[f[x], {x, a, c}, options] with a singularity
  55. at x == b.  This is accomplished as NIntegrate[f[x], {x, a, b-eps, b, b+eps, c},
  56. options], where the portion NIntegrate[f[x], {x, b-eps, b, b+eps}, options]
  57. is evaluated as NIntegrate[f[b+t] + f[b-t], {t, 0, eps}, options].  If eps
  58. does not appear, a default value will be used.  More generally,
  59. CauchyPrincipalValue[f[x], {x, x0, ..., xi, ..., {xj}, ..., xn}] gives
  60. the integral along the path {x0, ..., xi, ..., xj, ..., xn} using the
  61. Cauchy principle value at those abscissae enclosed in {}, any of which
  62. may include a user-specified eps."
  63.  
  64. Begin["NumericalMath`CauchyPrincipalValue`Private`"]
  65.  
  66. CauchyPrincipalValue[f_, range_List, options___] :=
  67.     Module[{e}, e /; (e = CPV0[f, range, options]) =!= Fail]
  68.  
  69. CPV0[f_, range_, options___] :=
  70.     Module[{a, b, eps, epsa, epsc, i, ilim, sum=0, temp, wprec},
  71.     ilim = Length[range]-1;
  72.     wprec = WorkingPrecision /. {options} /. Options[NIntegrate];
  73.     For[ i = 2, i <= ilim, i++,
  74.         If[ i == 2, 
  75.         a = range[[2]];
  76.         If[ Head[a] === List, Return[Fail]],
  77.           (* else *)
  78.         a = c
  79.         ];
  80.         b = range[[i+1]];
  81.         If[ Head[b] =!= List, 
  82.         c = b;
  83.         If[ a =!= c,
  84.             temp = NIntegrate[f, Evaluate[{range[[1]],a,c}], options];
  85.             If[ !NumberQ[temp], Return[Fail]];
  86.             sum += temp;
  87.         ];
  88.         Continue[]  (* next i in For[ ] loop *)
  89.         ];
  90.         c = range[[i+2]];
  91.         If[ Head[c] === List,
  92.         If[ i == ilim, Return[Fail]];
  93.         c = (c[[1]] + b[[1]])/2
  94.         ];
  95.         Which[
  96.         Length[b] > 2,
  97.             Return[Fail],
  98.         Length[b] == 2,
  99.             eps = b[[2]];
  100.             b = b[[1]],
  101.         True,
  102.             b = b[[1]];
  103.             eps = 1;
  104.             epsa = N[b-a, wprec]/2;
  105.             epsc = N[c-b, wprec]/2;
  106.             If[Abs[epsa] < 1, eps = epsa];
  107.             If[Abs[epsc] < Abs[eps], eps = epsc]
  108.         ];
  109.         temp = CPV1[f,{range[[1]],a,{b,eps},c},options];
  110.         If[ !NumberQ[temp], Return[Fail]];
  111.         sum += temp
  112.     ];
  113.     sum
  114.     ]
  115.  
  116. CPV1[f_, {x_, a_, {b_, eps_}, c_}, options___] :=
  117.     Module[{sum, t, fp, fm},
  118.     sum = Abs[N[Arg[b - a]] - N[Arg[b - c]]];
  119.     If[ sum > 1, sum = 2 N[Pi] - sum];
  120.     If[ sum < .00000000001,
  121.         NIntegrate[f, {x, a, c}, options],
  122.       (* else *)
  123.         sum = NIntegrate[f, {x, a, b-eps}, options];
  124.         sum += NIntegrate[f, {x, b+eps, c}, options];
  125.         fp = f /. x -> b+t;
  126.         fm = f /. x -> b-t;
  127.         sum + NIntegrate[fp + fm, {t, 0, eps}, options]
  128.     ]
  129.     ]
  130.  
  131. End[] (* NumericalMath`CauchyPrincipalValue`Private` *)
  132.  
  133. Protect[CauchyPrincipalValue];
  134.  
  135. EndPackage[] (* NumericalMath`CauchyPrincipalValue` *)
  136.  
  137.  
  138.