Integral Too Hard For Mathematica

Question

I have a monstrous integral that I desperately want to solve with Mathematica. It takes the form of:

Integrate[Exp[-2 t^2] t^(9/5) (x-t)^(4/5) HypergeometricPFQ[{1/2, 1}, {7/5, 19/10}, -t^2], {t, 0, x}]

where the numerical parameters can vary a bit. This is just one example. I'm keeping parameters rational in case Mathematica sees a trick it can use. What can I do? I would love to use NIntegrate but one of the bounds has to be variable. I also substituted it a Series approximation to HypergeometricPFQ, and that made it integrate, except the series introduced large error.

Thanks for any and all help.

Answer 1

This problem can be solved numerically by

s = ParametricNDSolveValue[{z'[t] == Exp[-2 t^2] t^(9/5) (x - t)^(4/5) 
    HypergeometricPFQ[{1/2, 1}, {7/5, 19/10}, -t^2], z[0] == 0}, z, {t, 0, x}, {x}];
Plot[s[x][x], {x, 0, 5}]

As noted by Marius Ladegård Meyer in a comment above, NDSolve itself is not sufficient, because x appears both in the integrand and as the upper limit of integration.

Answer 2

Here's an old-school (?) way that takes focused advantage of NDSolve's capabilities. See Components and Data Structures for detailed documentation.

ParametricNDSolve[] is very convenient, relatively recent, syntactic sugar and fairly well-suited for this problem, as @bbgodfrey's answer shows. In particular, it's easy to code (except for the [x][x] thing that surprised the OP). However, it integrates the solution for each distinct x. (It caches each solution in case it is reused, up to the cache size limit.) However the OP's problem can be solved with a single integration over a domain. Further that domain can be extended (not recomputed) by NDSolve using the documented internal components and data structures. I always thought this was a neat hook provided by WRI, but I hardly ever see it used.

Let f[x] represent the OP's integral. We can differentiate it to get another integral fp[x] that is equal to f'[x]. We can integrate the differential equation f'[x] == fp[x] with NDSolve up to whatever value of x and then extend it when needed. This can be done automatically with some simple code for the numerical version of f[x], which is called fN[x] below. (See Code Dump at end.)

There is one hitch to this method. It's a good idea to first call the function on the anticipated domain, so that NDSolve can integrate the function once for the domain. While extending the domain does not require the whole domain to be reintegrated, too many small extensions can be expensive. For instance,

fN[{0, 5}] // AbsoluteTiming
DownValues@fN`df // Length
(*
  {2.11563, {0., 0.47917}}
  77  <-- more on this later
*)

Then subsequent calls are very fast:

Plot[fN[x], {x, 0, 5}] // AbsoluteTiming

If you don't do the initialization of fN[] and try to Plot fN[x], many of the Plot calls to fN[x] end up extending the integration interval: that is, NDSolve is called many times to take unnecessarily short steps up to whatever value of x that Plot needs. This ends up calling the derivative function df, which uses NIntegrate[] to compute the derivative, a large number of times. (Re-execute the code defining fN[] if you wish to test the following; otherwise, it will be fast like above.)

Plot[fN[x], {x, 0, 5}] // AbsoluteTiming
DownValues@fN`df // Length

We can see that df[] was called 642 times instead of just 76 (one of the down values is the general definition). That accounts for the difference in speed.

Code Dump

I've gotten more used to writing packages, but my early experience has permanently scarred by brain with the impression that I don't know what I'm doing. Consequently, I conservatively remove from the "Global`" context all symbols I wish localized in my local context. (I'm assuming this code will be executed in a notebook in a running session as needed. If you properly pacletize it into a .wl file, prepping the "Global`" context should be unnecessary.)

I would guess that code for x < 0 is not needed by the OP (because f[x] is real for x >= 0 and complex for x < 0). Nonetheless, how to include it is shown. You need separate NDSolve calls. NDSolve will integrate a singly infinite interval {x, 0, ±Infinity}, but not a doubly infinite interval {x, -Infinity, Infinity}. (I don't know why; it makes no sense to me, since a finite interval {x, a, b} is automatically broken up into {x, x0, a} and {x, x0, b} when the IC is at x0 with a < x0 < b.)

The derivative df[] is memoized because it makes the code run twice as fast. (NDSolve must evaluate df[x] twice for each x on average, I guess.) If you're concerned about the amount of memory used by memoization (not a problem above), then you can free the memory with fN`freeDownValuesDF[]. It won't affect the subsequent speed, since the values are not re-used after the integration by NDSolve. (One could add a line to automatically call it after each NDSolve`Iterate[..] call.)

Quiet@Remove["fN`*"];
(* Removes symbols in Global` *)
Remove[fN, state, stateNeg, df, freeDownValuesDF, (* data/funcs *)
  f, \[Xi], x, t]; (* variables *)
ClearAll[fN]; (* create Global` symbol *)
fN // Attributes = {Listable};
Begin["fN`"]; (* Localize auxiliaries *)
(* Redundant with Remove[]; but creates local symbols in fN` *)
ClearAll[state, stateNeg, (* data *)
  f, \[Xi], x, t,(* variables *)
  df, freeDownValuesDF]; (* utilities *)

(* DERIVATIVE OF DESIRED FUNCTION f *)
With[{fp = 
    D[Inactive[Integrate][
      Exp[-2 t^2] t^(9/5) (x - t)^(4/5) HypergeometricPFQ[{1/2, 
         1}, {7/5, 19/10}, -t^2], {t, 0, x}], x]},
  (* memoization doubles the speed of NDSolve in this case *)
  mem : df[x0_?NumericQ] := mem = Block[{x = x0}, NIntegrate @@ fp];
  ];
(* release memory of memoization *)
freeDownValuesDF[] := DownValues[df] = {Last@DownValues@df};
(* END OF DERIVATIVE CODE *)

(* SOLUTION FOR POSITIVE x *)
{state} = NDSolve`ProcessEquations[
   {f'[\[Xi]] == df[\[Xi]], f[0] == 0},
   f, {\[Xi], 0, Infinity}
   (* NDSolve options go here (WorkingPrecision, Method) *)];
NDSolve`Iterate[state, 0.];
fN[x_?NumericQ] /; 
   x > state@"CurrentTime"["Forward"] := (NDSolve`Iterate[state, x];
   First@state@"SolutionVector"["Forward"]);
(* Solution for already solved x *)
fN[x_?NumericQ] /; x >= 0 := f[x] /. NDSolve`ProcessSolutions[state];
fN[x_?NumericQ] /; x < 0 := Print["Negative x not implemented."];
(* END OF POSITIVE x CODE *)

(* OMIT IF NEGATIVE x NOT ALLOWED *)
{stateNeg} = NDSolve`ProcessEquations[
   {f'[\[Xi]] == df[\[Xi]], f[0] == 0},
   f, {\[Xi], 0, -Infinity}
   (* NDSolve options go here (WorkingPrecision, Method) *)];
NDSolve`Iterate[stateNeg, 0.];
fN[x_?NumericQ] /; 
   x < stateNeg@"CurrentTime"["Backward"] := (NDSolve`Iterate[
    stateNeg, x];
   First@state@"SolutionVector"["Backward"]);
(* Solution for already solved x; overwrites previous def *)
fN[x_?NumericQ] /; x < 0 := f[x] /. NDSolve`ProcessSolutions[stateNeg];
(* END OF NEGATIVE x CODE *)

End[];

Answer 3

Problem solved with approximation by sum:

$$\int_0^x \exp \left(-2 t^2\right) t^{9/5} (x-t)^{4/5} \, _2F_2\left(\frac{1}{2},1;\frac{7}{5},\frac{19}{10};-t^2\right) \, dt=\\\int_0^x \exp \left(-2 t^2\right) t^{9/5} (x-t)^{4/5} \sum _{k=0}^{\infty } \frac{\left(\left(\frac{1}{2}\right)_k (1)_k\right) \left(-t^2\right)^k}{\left(\left(\frac{7}{5}\right)_k \left(\frac{19}{10}\right)_k\right) k!} \, dt=\\\sum _{k=0}^{\infty } \int_0^x \frac{(-1)^k e^{-2 t^2} t^{\frac{9}{5}+2 k} (-t+x)^{4/5} \left(\frac{1}{2}\right)_k (1)_k}{k! \left(\frac{7}{5}\right)_k \left(\frac{19}{10}\right)_k} \, dt=\\\frac{\Gamma \left(\frac{9}{5}\right) \Gamma \left(\frac{14}{5}\right) \sum _{k=0}^{\infty } \Gamma \left(\frac{1}{2}+k\right) (-1)^k x^{\frac{18}{5}+2 k} \, _2\tilde{F}_2\left(\frac{7}{5}+k,\frac{19}{10}+k;\frac{23}{10}+k,\frac{14}{5}+k;-2 x^2\right)}{8\ 2^{3/5}}$$

f[x_, M_] := 1/(8 2^(3/5))*Gamma[9/5] Gamma[14/5] Sum[
Gamma[1/2 + k] (-1)^k x^(18/5 + 2 k)*
HypergeometricPFQRegularized[{7/5 + k, 19/10 + k}, {23/10 + k, 
  14/5 + k}, -2 x^2], {k, 0, M}]

 s = ParametricNDSolveValue[{z'[t] == 
 Exp[-2 t^2] t^(9/5) (x - t)^(4/5) HypergeometricPFQ[{1/2, 
    1}, {7/5, 19/10}, -t^2], z[0] == 0}, z, {t, 0, x}, {x}];
 Plot[{s[x][x], f[x, 5]}, {x, 0, 5}, 
 PlotStyle -> {Black, {Dashed, Red}},PlotLabels -> {"NDSolve", "Series"}]