(* ::Package:: *) (*:Name: NumericalMath`CauchyPrincipalValue` *) (*:Title: Numerical Evaluation of Cauchy Principal Values *) (*:Author: Jerry B. Keiper *) (*:Summary: This package introces a function for numerically evaluating Cauchy principal values of integrals that may not be Riemann integrable. The built-in function NIntegrate is used in a symmetric way to get cancelation. *) (*:Context: NumericalMath`CauchyPrincipalValue` *) (*:Package Version: Mathematica 2.0 *) (* :Copyright: Copyright 1990-2007, Wolfram Research, Inc. *) (* :History: Version 1.2 by Jerry B. Keiper, January 1990. Revised by Jerry B. Keiper, November 1990. *) (*:Keywords: Cauchy principal value, integration, divergent integral *) (*:Source: E. T. Whittaker & G. N. Watson, A Course of Modern Analysis, Cambridge University Press, (4th ed.), Section 4.5. *) (*:Mathematica Version: 2.0 *) (*:Limitation: Since CauchyPrincipalValue is a numerical function and it relies on symmetry near the singular points, loss of digits due to cancelation can be a problem. The user must know and specify EXACTLY where the singularities are, or the result can be grossly in error. Again this is because of the extreme sensitivity to perturbations. *) Message[General::obspkg, "NumericalMath`CauchyPrincipalValue`"] BeginPackage["NumericalMath`CauchyPrincipalValue`"] CauchyPrincipalValue::usage = "CauchyPrincipalValue[f[x], {x, a, {b, eps}, c}, options] gives the Cauchy \ principal value of NIntegrate[f[x], {x, a, c}, options] with a singularity \ at x == b. This is accomplished as NIntegrate[f[x], {x, a, b-eps, b, b+eps, c}, \ options], where the portion NIntegrate[f[x], {x, b-eps, b, b+eps}, options] \ is evaluated as NIntegrate[f[b+t] + f[b-t], {t, 0, eps}, options]. If eps \ does not appear, a default value will be used. More generally, \ CauchyPrincipalValue[f[x], {x, x0, ..., xi, ..., {xj}, ..., xn}] gives \ the integral along the path {x0, ..., xi, ..., xj, ..., xn} using the \ Cauchy principle value at those abscissas enclosed in {}, any of which \ may include a user-specified eps." Begin["NumericalMath`CauchyPrincipalValue`Private`"] issueObsoleteFunMessage[fun_, context_] := (Message[fun::obspkgfn, fun, context]; ) CauchyPrincipalValue[f_, range_List, options___] := (issueObsoleteFunMessage[CauchyPrincipalValue,"NumericalMath`CauchyPrincipalValue`"]; Module[{e}, e /; (e = CPV0[f, range, options]) =!= Fail]) CPV0[f_, range_, options___] := Module[{a, b, eps, epsa, epsc, i, ilim, sum=0, temp, wprec}, ilim = Length[range]-1; wprec = WorkingPrecision /. {options} /. Options[NIntegrate]; For[ i = 2, i <= ilim, i++, If[ i == 2, a = range[[2]]; If[ Head[a] === List, Return[Fail]], (* else *) a = c ]; b = range[[i+1]]; If[ Head[b] =!= List, c = b; If[ a =!= c, temp = NIntegrate[f, Evaluate[{range[[1]],a,c}], options]; If[ !NumberQ[temp], Return[Fail]]; sum += temp; ]; Continue[] (* next i in For[ ] loop *) ]; c = range[[i+2]]; If[ Head[c] === List, If[ i == ilim, Return[Fail]]; c = (c[[1]] + b[[1]])/2 ]; Which[ Length[b] > 2, Return[Fail], Length[b] == 2, eps = b[[2]]; b = b[[1]], True, b = b[[1]]; eps = 1; epsa = N[b-a, wprec]/2; epsc = N[c-b, wprec]/2; If[Abs[epsa] < 1, eps = epsa]; If[Abs[epsc] < Abs[eps], eps = epsc] ]; temp = CPV1[f,{range[[1]],a,{b,eps},c},options]; If[ !NumberQ[temp], Return[Fail]]; sum += temp ]; sum ] arg[z_?(#!=0&)]:=Arg[z] arg[z_?(#==0&)]:=0 CPV1[f_, {x_, a_, {b_, eps_}, c_}, options___] := Module[{sum, t, fp, fm}, sum = Abs[N[arg[b - a]] - N[arg[b - c]]]; If[ sum > 1, sum = 2 N[Pi] - sum]; If[ sum < .00000000001, NIntegrate[f, {x, a, c}, options], (* else *) sum = NIntegrate[f, {x, a, b-eps}, options]; sum += NIntegrate[f, {x, b+eps, c}, options]; fp = f /. x -> b+t; fm = f /. x -> b-t; sum + NIntegrate[fp + fm, {t, 0, eps}, options] ] ] End[] (* NumericalMath`CauchyPrincipalValue`Private` *) Protect[CauchyPrincipalValue]; EndPackage[] (* NumericalMath`CauchyPrincipalValue` *)