Re: AiryAi
- To: mathgroup at smc.vnet.net
- Subject: [mg75817] Re: AiryAi
- From: antononcube <antononcube at gmail.com>
- Date: Fri, 11 May 2007 05:28:04 -0400 (EDT)
- References: <f1hlg0$p7q$1@smc.vnet.net><f1rvv6$8ma$1@smc.vnet.net>
In Mathematica 6.0 one can define a general integration strategy that
handles inifinite range oscillatory integrals through extrapolation of
the partial sums. Here is an example how to use it (the definition is
below):
In[109]:= NIntegrate[Sin[x]/x, {x, 1, Infinity}, Method ->
{Extrapolating, "ZeroesFunction" -> (#1*Pi & )}]
Out[109]= 0.6247132564277149
In[8]:= NIntegrate[Sin[x]/x, {x, 1, Infinity}, Method ->
{Extrapolating, "ZeroesFunction" -> (#1*Pi & )}, WorkingPrecision ->
30]
Out[8]= 0.62471325642771360552793845678`17.4518826436565
The strategy is called Extrapolating. It takes as an option a function
that gives the zeroes of the integrand. (Actually, it is not necessary
to integrate between the zeroes, we can integrate between the
consecutive points of any sequence, and then use extrapolation.)
Below are examples with the AiryAi function. (Note that Mathematica
has the function AiryAiZero function):
(* this uses Dimitris code *)
In[158]:= dat = Range[-40, -2, 0.5];
zer = Reverse[
Append[Union[(o /. FindRoot[AiryAi[o] == 0, {o, #1}] &) /@ dat,
SameTest -> (#1 - #2 < 10^(-12) &)], 0]];
In[160]:= NIntegrate[AiryAi[x], {x, 0, -Infinity},
Method -> {Extrapolating, "ZeroesFunction" -> (dat[[#]] &)}]
During evaluation of In[160]:= SequenceLimit::seqlim: The general \
form of the sequence could not be determined, and the result may be \
incorrect. >>
Out[160]= -0.665436
(* this uses Mathematica's AiryAiZero *)
In[161]:= NIntegrate[AiryAi[x], {x, 0, -Infinity}, Method ->
{Extrapolating,
"ZeroesFunction" -> (AiryAiZero[Rationalize@#] &)}]
Out[161]= -0.666667
The strategy`s code below can be improved in various ways, and I am
sure it might have some problems with inputs it is not intended for.
Anton Antonov,
Wolfram Research, Inc.
=====================
Extrapolating strategy code
=====================
Options[Extrapolating] = {"ZeroesFunction" -> None,
"SymbolicProcessing" -> 0};
ExtrapolatingProperties = Options[Extrapolating][[All, 1]];
Extrapolating /:
NIntegrate`InitializeIntegrationStrategy[Extrapolating, nfs_,
ranges_, strOpts_, allOpts_] :=
Module[{t, zerofunc, symbproc, a, b},
t = NIntegrate`GetMethodOptionValues[Extrapolating,
ExtrapolatingProperties, strOpts];
If[t === $Failed, Return[$Failed]];
{zerofunc, symbproc} = t;
If[zerofunc === None,(*message*)Return[$Failed]];
If[Length[ranges] > 1,(*message*)Return[$Failed]];
{a, b} = Rest[ranges[[1]]];
If[! (NumberQ[a] && Head[b] === DirectedInfinity),(*message*)
Return[$Failed]];
Extrapolating[{None, a, Infinity, zerofunc}]
];
Extrapolating[{ruleArg_, a_, Infinity, zerofunc_}]
["Algorithm"[regionsArg_, opts___]] :=
Module[{t, term, integral, func, var, nf, wprec =
WorkingPrecision /. opts, nopts},
nf = regionsArg[[1]]["Integrand"];
func = nf["FunctionExpression"];
var = nf["ArgumentNames"][[1]];
nopts = {Method -> {Automatic, "SymbolicProcessing" -> 0},
WorkingPrecision -> wprec};
integral = NIntegrate[func, {var, a, N[zerofunc[1], wprec]} //
Evaluate, Evaluate[Sequence @@ nopts]];
term[i_?NumberQ] :=
NIntegrate[
func, {var, N[zerofunc[i], wprec], N[zerofunc[i + 1], wprec]} //
Evaluate, Evaluate[Sequence @@ nopts]];
integral =
integral +
NSum[term[i], {i, 1, Infinity}, Method -> "WynnEpsilon",
WorkingPrecision -> wprec];
integral
];