Re: hadamard finite part
- To: mathgroup at smc.vnet.net
- Subject: [mg68023] Re: hadamard finite part
- From: "antononcube" <antononcube at gmail.com>
- Date: Thu, 20 Jul 2006 06:05:03 -0400 (EDT)
- References: <e9ibvc$52e$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Below is an implementation with examples based on the algorithm in
"Computational Integration" by Krommer & Ueberhuber, 1998, pp 18-19.
I have made the definitions to work with symbolic singular roots, and
to use one of the integrators Integrate or NIntegrate.
Anton Antonov,
Wolfram Research, Inc.
Clear[HFP];
HFP[f_, {x_, a_, b_}] :=
Module[{y},
y = f /. \!\(\((x + y1_. )\)\^\(-1\)\) -> -y1;
Which[
TrueQ[y == a], Log[b - y],
TrueQ[y == b], -Log[y - a],
True(*a<y<b*), Log[(b - y)/(y - a)]
]
] /; MatchQ[f, 1/(x + y_.)];
HFP[f_, {x_, a_, b_}] :=
Module[{y, \[Alpha]},
{y, \[Alpha]} = f /. ((x + y1_.))^\[Alpha]1_. -> {-y1, -\[Alpha]1};
\[Alpha] = \[Alpha] - 1;
Which[
TrueQ[y == a], -1/(\[Alpha] ((b - y))^\[Alpha]),
TrueQ[y == b], 1/(\[Alpha] ((y - a))^\[Alpha]),
True(*a<y<
b*), -1/\[Alpha] (1/(((b - y))^\[Alpha]) - 1/(((y - a))^\[Alpha]))
]] /; (MatchQ[f, ((x + y_.))^\[Alpha]_.]);
HFP[F_, {x_, a_, b_}, int_: Integrate] :=
Module[{f, y, \[Alpha], s, k, t},
{f, y, \[Alpha]} =
F /. f1_. ((x + y1_.))^\[Alpha]1_. -> {f1, -y1, -\[Alpha]1};
HFP[f, \!\(\((x - y)\)\^\(-\[Alpha]\)\), {x, a, b}, int]
];
HFP[f_, d_, {x_, a_, b_}, int_: Integrate] :=
Module[{y, \[Alpha], s, k, t},
{y, \[Alpha]} = d /. ((x + y1_.))^\[Alpha]1_. -> {-y1, -\[Alpha]1};
\[Alpha] = \[Alpha] - 1;
(*Print["{y,\[Alpha]}=",{y,\[Alpha]}];*)
If[\[Alpha] <= 0, Return[$Failed]];
k = Ceiling[\[Alpha]];
s = Total@
Table[((\!\(\[PartialD]\_{x, i}f\))/(i!) /.\[InvisibleSpace]x ->
y)
HFP[1/(((x - y))^(\[Alpha] + 1 - i)), {x, a, b}], {i, 0, k}];
t = Total@
Table[((\!\(\[PartialD]\_{x, i}f\))/(
i!) /.\[InvisibleSpace]x -> y) ((x - y))^i, {i, 0, k}];
int[(f - t)/(((x - y))^(\[Alpha] + 1)), {x, a, b}] + s
] /; MatchQ[
d, ((x + y_.))^\[Alpha]_.] && (int === NIntegrate || int ==
Integrate);
Clear[x, \[Lambda]]
In[124]:= HFP[1/(x - \[Lambda]), {x, \[Lambda], 1}]
(Out[124]= Log[1 - \[Lambda]]
In[125]:= HFP[1/(x - \[Pi]/4), {x, \[Pi]/4, 1}]
Out[125]= Log[1 - \[Pi]/4]
In[126]:= HFP[1/(x - \[Lambda]), {x, 0, \[Lambda]}]
Out[126]= -Log[\[Lambda]]
In[127]:= HFP[1/(((x - \[Lambda]))^(\[Alpha] + 1)), {x, 0, 1}]
Out[127]= -(\!\(\((1 - \[Lambda])\)\^\(-\[Alpha]\)\) -
\!\(\[Lambda]\^\(-\[Alpha]\)\))/\
\[Alpha]
In[128]:= HFP[1/(((x - \[Pi]/4))^2), {x, 0, 1}]
Out[128]= -1/(1 - \[Pi]/4) + 4/\[Pi]
In[129]:= HFP[1/(((x))^2), {x, 0, 1}]
(Out[129]= -1
In[130]:= HFP[(1 + x + x^2)/(x^2), {x, 0, 2}]
Out[130]= 3/2 + Log[2]
Davis & Rabinowitz example, "Methods of Numerical Integration", 2nd
ed., 1984, pp 190
In[132]:= res = HFP[(\!\(y\^\(-2\)\))/(((((y - 2))^2 + 1))^(1/2)), {y,
0, 1},
Integrate]
% // N
Out[132]= -1/(Sqrt[5]) + (3 - Sqrt[10] - 2 ArcSinh[3] + Log[100])/(5
Sqrt[5])
Out[133]= -0.375123
In[134]:= res = HFP[(\!\(y\^\(-2\)\))/(((((y - 2))^2 + 1))^(1/2)), {y,
0, 1},
NIntegrate]
Out[134]= -0.375123
In[135]:= exact = (2 Sqrt[5])/25 (Log[10 Sqrt[10] - 30] - 1) -
(Sqrt[2])/5;
In[136]:= res - exact // N
Out[136]= -2.53964Ã?\!\(10\^\(-15\)\)
dimmechan at yahoo.com wrote:
> I work in the field of applied mathematics and I am interested in the
> symbolical/numerical integration of integrals in the Hadamard sense
> (that is, the finite part of divergent integrals).
> My integrals are much more complicated but here I use some trivial
> examples to show the point.
> [For clarity, I have converted the output expressions to
> InputForm]
>
> 5.2 for Microsoft Windows (June 20, 2005)
>
> Consider first the integral,
>
> Integrate[1/x,{x,0,1}]
> 1
> Integrate::idiv: Integral of - does not converge on {-1, 2}.
> More�
> x
>
> Integrate[x^(-1), {x, -1, 2}]
>
>
> which does not exist in the Reimann sense, since it has a
> non-integrable singularity at x=0.
>
> However in the Cauchy sense the integral exists,
>
> Integrate[1/x,{x,-1,2},PrincipalValue->True]
> Log[2]
>
> N[%]
>
> 0.6931471805599453
>
> which is equivalent to
>
> Integrate[1/x,{x,-1,-e1},Assumptions0<e1<1]+Integrate[1/x,{x,e2,2},Assumptions->0<e2<2]
>
> Log[e1] + Log[2/e2]
>
> PowerExpand[%]
>
> Log[2] + Log[e1] - Log[e2]
>
> and then taking e1=e2
>
> %/.e1->e2
>
> Log[2]
>
> Limit[%,e2->0,Direction->-1]
>
> Log[2]
>
> Of course the Mathematica is also able to evaluate numerically the
> principal value of
> this integral,
>
> Needs["NumericalMath`"]
>
> CauchyPrincipalValue[1/x,{x,-1,{0},2}]
>
> NIntegrate::ploss: Numerical integration stopping due to loss of
> precision. Achieved neither the requested PrecisionGoal nor
> AccuracyGoal; suspect one of the following: highly oscillatory
> integrand or the true value of the integral is 0. If your
> integrand is oscillatory on a (semi-)infinite interval try using
> the option Method->Oscillatory in NIntegrate. More�
>
> 0.6931471805602387
>
> So everything is good up to now.
>
> Next, suppose the integral
>
> Integrate[1/x^2,{x,-1,2}]
>
> -2
> Integrate::idiv: Integral of x does not converge on {-1, 2}.
> More�
>
> Integrate[x^(-2), {x, -1, 2}]
>
> There is a again non-integrable singularity at x=0. Due to the kind of
> the singularity (double pole) the integral does not exist even in the
> Cauchy sense.
>
> Integrate[1/x2,{x,-1,-e1},Assumptions->0<e1<1]+Integrate[1/x2,{x,e2,2},Assumptions->0<e2<=2]
>
> -3/2 + e1^(-1) + e2^(-1)
>
> %/.e1->e2
>
> -3/2 + 2/e2
>
> Limit[%,e2->0,Direction->-1]
>
> Infinity
>
>
> Nevertheless, the integral exists in the Hadamard sense and is equal to
> finite part of the previous result
>
> %%%-( e1^(-1) + e2^(-1))
>
> -3/2
>
> Mathematica is able to provide previous result:
>
> Integrate[1/x^2,{x,-1,2},GenerateConditions->False]
>
> -3/2
>
> My first question now:
> Is it a way to get the finite part of a divergent integral through
> performing
> numerical integration (e.g. using NIntegrate) in Mathematica?
> I have seen some papers presenting some propper algorithms dealing with
> numerical integration of Hadamard finite part integrals but I cannot
> find any related work in connection with Mathematica.
>
> There are also some other questions.
>
> Consider the integral,
>
>
> Integrate[1/x,{x,0,2}]
>
> 1
> Integrate::idiv: Integral of - does not converge on {0, 2}. More�
> x
>
> Integrate[x^(-1), {x, 0, 2}]
>
>
> The integral does not exist neither in the Reimann sense nor in the
> Cauchy
> sense. However it does exist in the Hadamard sense and is equal to the
> finite part of
>
> Integrate[1/x,{x,e,2},Assumptions->0<x<2]//PowerExpand
>
> Log[2] - Log[e]
>
> %-Log[e]
>
> Log[2]
>
> Mathematica 3.0 and 4.0 suceeds in providing this result:
>
> Integrate[1/x,{x,0,2},GenerateConditions->False] (*version 3.0 and
> 4.0*)
>
> Log[2]
>
> However Mathematica 5.1 and 5.2 gives the result
>
> Integrate[1/x,{x,0,2},GenerateConditions->False] (*version 5.1 and
> 5.2*)
>
> Infinity
>
> Why exists this difference?
> I can trust that for divergent integrals
> Integrate[integrand,{x,a,b},GenerateConditions->False]
> provides the desirable result in the Hadamard sense?
> Is it a way to get Integrate to give always the finite part of a
> divergent integral?
> Are there any other alternative methods (such as the implementation of
> aymptotic techniques in mathematica) to get the finite part of a
> divergent integral?
> Is there any book that connect Mathematica with distribution theory?
>
> I really appreciate any assistance.
>
> P.S. The finite part of a divergent integral is of great importance in
> the area of applied mathematics. For example, the vertical displacement
> at the surface of an elastic linear isotropic half space acted upon by
> an applied line load (the so called Flamant problem) is given by the
> Fourier cosine transform of the function 1/ξ, which does not exist in
> the Reimann sense (the integral has a non integrable singularity at
> ξ=0). Indeed,
>
> Series[Cos[j]/j,{j,0,3}]//Normal
>
> \!\(TraditionalForm\`x\^3\/24 - x\/2 + 1\/x\)
>
> However in the Hadamard sense it does exist and is equal to
>
> Cancel[PowerExpand[Integrate[(1/Î?Î?)Cos[
> Î?Î? x],{Î?Î?,0,βÂ?Â?},GenerateConditions\[Rule]False]]]
>
> -EulerGamma - Log[x]
>
> or
>
> FullSimplify[Sqrt[Pi/2]FourierCosTransform[1/j,j,x],x>0]
>
> -EulerGamma - Log[x]
>
> This is a common practice to elasticity problems and the "loss" of
> infinity correcponds to rigid-body displacement (cf. e.g. Horacio A.
> Sosa and Leon Y. Bahar: "Transition from airy stress function to state
> space formulation of elasticity" Journal of the Franklin Institute,
> Volume 329, Issue 5, September 1992, Pages 817-828)