Problem converting double integral to PDE and solving with NDSolve

Question

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!

Answer 1

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