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!

up vote 2 down vote accepted

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! – Loki Nov 10 at 22:13
  • 2
    I recommend using a parametric function XI[r][k] instead of solving PDE. – Alex Trounev Nov 10 at 22:30

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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