# Problem converting double integral to PDE and solving with NDSolve

I'm trying to use the well-known "trick" of converting an integral to a differential equation, which can be much faster if the integral has to be evaluated repeatedly (in my case millions of times), but it seems there is a problem.

The integral is the following: $$\xi(r)=\frac{1}{2\pi^2}\int_{k_{min}}^{k_{max}} P(k) \frac{sin(kr)}{kr}k^2 dk,$$ where $$k_{min}=10^{-4}$$, $$k_{max}=100$$, the parameter $$r$$ runs in the range $$r\in[10,150]$$ and $$P(k)$$ is a complicated function that can only be determined numerically, but for the sake of this example we can assume it is reasonably close to a log-normal of the form: $$P(k)=\frac{108400 \sqrt{\frac{2}{\pi }} e^{-\frac{5000 \left(\log (k)+\frac{283}{100}\right)^2}{14641}}}{121 k}$$

Calculating the integral and converting it to an ODE is simple:

(* Some parameters *)
kmin = 10^-4;
kmax = 100;
rmin = 10;
rmax = 150;

(* The complicated function *)
P[k_] := A E^(-((-m + Log[k])^2/(2 s^2)))/(
k Sqrt[2 \[Pi]] s) /. {m -> -283/100, s -> 121/100, A -> 2168};

(* Implementation of the integral with NIntegrate *)
xi0[r_] :=
NIntegrate[1/(2 Pi^2) P[k] Sin[k r]/(k r) k^2, {k, kmin, kmax},
Method -> Automatic, WorkingPrecision -> 20, PrecisionGoal -> 15,
AccuracyGoal -> 15, MaxRecursion -> 15]

(* Converting the integral to an ODE and solving with NDSolve *)
xi1[r_] :=
xi[kmax] /.
NDSolve[{xi'[k] == 1/(2 Pi^2) P[k] Sin[k r]/(k r) k^2,
xi[kmin] == 0}, xi, {k, kmin, kmax}, WorkingPrecision -> 20,
PrecisionGoal -> 15, AccuracyGoal -> 15, MaxSteps -> 10^5][[1]]


Both NIntegrate and NDSolve give results in good agreement for $$r=50$$:

(* Result from NIntegrate *)
In[39]:= xi0[50]
Out[39]= 0.0079866930495540006015


and

(* Result from NDSolve *)
In[40]:= xi1[50]
Out[40]= 0.0079866930482492951647


However, I want to evaluate them for million values of the parameter $$r$$, so why not hit the ODE with a derivative with respect to $$r$$, convert it to a PDE and try NDSolve again. Then I have:

(* Convert to PDE by hitting the previous equation with derivative wrt r *)
In[48]:= sol=NDSolve[{D[xi[k, r], k, r] ==
1/(2 Pi^2) P[k] (Cos[k r]/r - Sin[k r]/(k r^2)) k^2,xi[kmin, r]
== 0, xi[k, rmax] == 0},xi,{k,kmin,kmax},{r,rmin,rmax}];


Which however for $$r=50$$ gives the wrong result:

(* Result from PDE *)
In[47]:= xi[kmax, 50] /. sol[[1]]
Out[47]= -0.0362161


Any ideas what went wrong in this case? Any help is much appreciated!

The boundary conditions are specified incorrectly and the size of the grid does not correspond to a rapidly oscillating function. However, even with correctly defined boundary conditions and a fine grid, the accuracy leaves much to be desired. Example

(*Some parameters*)kmin = 10^-4;
kmax = 100;
rmin = 10;
rmax = 150;

(*The complicated function*)
P[k_] := A E^(-((-m + Log[k])^2/(2 s^2)))/(k Sqrt[
2 \[Pi]] s) /. {m -> -283/100, s -> 121/100, A -> 2168};

(*Implementation of the integral with NIntegrate*)
xi0[r_] :=
NIntegrate[1/(2 Pi^2) P[k] Sin[k r]/(k r) k^2, {k, kmin, kmax},
Method -> Automatic, WorkingPrecision -> 20, PrecisionGoal -> 15,
AccuracyGoal -> 15, MaxRecursion -> 15]

(*Converting the integral to an ODE and solving with NDSolve*)
XI = ParametricNDSolveValue[{xi'[k] ==
1/(2 Pi^2) P[k] Sin[k r]/(k r) k^2, xi[kmin] == 0},
xi, {k, kmin, kmax}, {r}, WorkingPrecision -> 20,
PrecisionGoal -> 15, AccuracyGoal -> 15, MaxSteps -> 10^6];
(*Convert to PDE by hitting the previous equation with derivative wrt \
r*)sol = NDSolveValue[{D[xi2[k, r], k, r] ==
1/(2 Pi^2) P[k] (Cos[k r]/r - Sin[k r]/(k r^2)) k^2,
xi2[kmin, r] == 0, xi2[k, rmax] == XI[rmax][k]},
xi2, {k, kmin, kmax}, {r, rmin, rmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 10000, "MaxPoints" -> 100000,
"DifferenceOrder" -> "Pseudospectral"}}, MaxSteps -> 10^6];
In[9]:= xi0[50]
XI[50][kmax]

Out[9]= 0.0079866930495540006015

Out[10]= 0.0079866930493438685592

In[11]:= sol[kmax, 50]

Out[11]= 0.00781333

• Many thanks for the response! I'm worried a bit as the solution of the PDE is quite slower than NIntegrate in this case (probably due to using the correct boundary conditions which require knowing xi(r)) and it seems to take up lots of RAM. I'll see if it can be improved somehow. In any case, thanks again! – Rebel-Scum Nov 10 '18 at 22:13
• I recommend using a parametric function XI[r][k] instead of solving PDE. – Alex Trounev Nov 10 '18 at 22:30