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!