2
$\begingroup$

I have an integral in which the integrand is a combination of ratio of modified bessel function of second kind. I am able to find the solution of this integral with mathematica Nintegrate module with global adaptive method but i want to use this global adaptive method code to write a separate program in fortran. Can i do that?

enter image description here

$\endgroup$
  • 1
    $\begingroup$ You say you have been able to use NIntegrate to compute these integrals. Then it would be a good idea to share your code with us, no? $\endgroup$ – Henrik Schumacher May 31 '19 at 10:48
  • 3
    $\begingroup$ So this is not essentially about Mathematica, but just how to do an integration in Fortran. $\endgroup$ – corey979 May 31 '19 at 12:15
  • 2
    $\begingroup$ An example implementation is given here. In this answer, I adapted it for vector arguents. You can also find implementations in various languages/systems/pseudocode in some introductory numerical analysis books. Then you just have to translate the algorithm into Fortran. $\endgroup$ – Michael E2 May 31 '19 at 12:32
  • 1
    $\begingroup$ Depending on your code, Mathematica would probably have chosen the Levin rule, which is considerably different than the Gauss-Kronrod rule used in the example in the tutorial I linked my previous comment. $\endgroup$ – Michael E2 May 31 '19 at 12:37
  • 1
    $\begingroup$ I see essentially two strategies. 1) The integral is evaluated in MA for a range of parameters, Pade-approximated and given in this form to Fortran. 2) MA is only used for benchmarking. All evaluations are done in Fortran. I recommend second option, using GSL library gnu.org/software/gsl/doc/html/index.html . It contains subroutines for quadratures and for modified Bessel functions. They are callable from C, but mixing C and Fortran is easy. $\endgroup$ – yarchik Jun 1 '19 at 9:02
3
$\begingroup$

That's not possible. That's not so easy. NIntegrate is probably too complex to be automagically ported to FORTRAN or any other programming language. Depending on the effort that you are willing to invest, there are several options:

  1. You may use Mathematica to find out which fixed quadrature rule works best for this family of integrands. Afterwards, you will have to implement this quadrature rule yourself in FORTRAN. Or you better use a library that is specialized on computing numerical integrals (as you will also have to use a library for evaluating the modified Bessel functions). From the fact that you are working in FORTRAN, I deduce that you are interested mostly in runtime performance. If you have to compute really many of such integrals and if these computations are time-critical, this is probably the most efficient way; in contrast adaptive methods tend to have quite an overhead.

  2. You may use WSTP to set up communication between FORTRAN and a Mathematica kernel. I do not have any experience with this. IIRC, this will require you (i) to write a wrapper library in C that sets up the WSTP connection and (ii) to call the generated library from FORTRAN. This will require a running Mathematica kernel and thus an appropriate license. (The latter might be mended by using the free Wolfram Engine for developers; I am not sure about this.) In any case, this won't be very efficient.

  3. As JimB pointed out, the commercial addon MathCode F90 might do the job. I cannot tell, because I have not used it yet. But I somewhat doubt that it will work as expected. As far as I know, only a subset of Mathematica code is actually compilable by MathCode F90 (Mathematica's Compile and FunctionCompile have similar limitations). For example, it is stated on this page that FindMinimum is not compilable. NIntegrate is similarly complex, so I do believe that it is not compilable either. You might find more details in Appendix A of the manual.

  4. If you require the integrals only for an a a priorily know, low-dimensional range of parameters, then you may try to compute the integrals once on a fixed grid of points with Mathematica, export the data to file, load it into FORTRAN (at runtime), and use interpolation to approximate the integrals for arbitrary parameters. This might be very fast, even faster than method 1 because Bessel functions are notoriously expensive to evaluate. However, this requires that the integrals depend smoothly on the parameters.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Are you saying that MathCode F90 can’t do what the OP is needing? $\endgroup$ – JimB May 31 '19 at 8:34
  • $\begingroup$ @JimB Good point. It is always a bad idea to state that something would be impossible. I haven't thought about MathCode F90 (because I do not have any experience with it). $\endgroup$ – Henrik Schumacher May 31 '19 at 10:44
  • $\begingroup$ But i think you’re still correct now that I’ve also looked again at the literature for MathCode F90. (I have to wonder if MATLAB has a similar issue.) $\endgroup$ – JimB May 31 '19 at 15:17
  • $\begingroup$ How can we link Fortran with mathematica by using WSTP? Give an example. $\endgroup$ – Prabhakar Namdev Jun 7 '19 at 10:07
  • 1
    $\begingroup$ @PrabhakarNamdev Nope, I won't do that. Because (i) I am not an expert (and I haven't tried, yet) and because (ii) you certainly won't be happy with the performance. $\endgroup$ – Henrik Schumacher Jun 7 '19 at 10:41
2
$\begingroup$

I recommend that you look at numerical libraries for FORTRAN. I haven't used FORTRAN in a quarter of a century, so I'm not familiar with what's available today, but this is the kind of problem FORTRAN was designed for. I expect you'll find something that fits your needs.

It wouldn't surprise me if NIntegrate works in this case by first symbolically analyzing the problem (very difficult to automate in FORTRAN) and then handing it to a chosen FORTRAN library function. You'll have to do that choosing yourself.

| improve this answer | |
$\endgroup$

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