10
$\begingroup$

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.

$\endgroup$
9
  • $\begingroup$ Try to solve the corresponding differential equation with NDSolve. $\endgroup$ Aug 11, 2016 at 20:36
  • $\begingroup$ E[-2 t^2] is not defined. Did you mean Exp[-2 t^2]? $\endgroup$
    – QuantumDot
    Aug 11, 2016 at 20:38
  • 2
    $\begingroup$ I was thinking about the fact that if $f(x) = \int_k^x g(t) dt$ then $f(x) = G(x) - G(k)$ where $G'(x) = g(x)$. Mathematica can't do e.g. Integrate[Sin[t]^Exp[t], {t, 0, x}] but it can do NDSolve[f'[x] == Sin[x]^Exp[x] && f[0] == 0, f, {x, 0 10}] just fine. But maybe this is too simple since you have x in the integrand as well... $\endgroup$ Aug 11, 2016 at 22:12
  • 4
    $\begingroup$ I think @MariusLadegårdMeyer is talking about converting your integration into the solution of a differential equation, which is a nice trick in Mathematica. I.e. NIntegrate[f[x], {x, 0, t}] is equivalent to NDSolveValue[{g'[x] == f[x], f[0] == 0}, g, {x, 0, 50}][t] (as long as you choose 50 to be a large enough number, beyond the maximum t you use). The NDSolve version is usually much faster than the NIntegrate version. $\endgroup$
    – march
    Aug 11, 2016 at 22:12
  • 1
    $\begingroup$ @march, can you post it as an answer, like to upvote. $\endgroup$
    – bobbym
    Aug 12, 2016 at 4:12

3 Answers 3

16
$\begingroup$

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

enter image description here

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.

$\endgroup$
2
  • $\begingroup$ Beautiful job, thank you! Just out of curiosity, what does the second pair of [x] on s do in the Plot command? $\endgroup$
    – Buddhapus
    Aug 12, 2016 at 17:23
  • $\begingroup$ @Buddhapus The first [x] gives the value of the parameter {x} used to solve the ODE. The second [x] gives the value of t for which the solution is desired, which happens to be x in this case.. $\endgroup$
    – bbgodfrey
    Aug 12, 2016 at 17:28
8
$\begingroup$

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

enter image description here

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

enter image description here

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[];
$\endgroup$
6
$\begingroup$

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"}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.