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 ];