numerical integration of a complex function?

[thejasvi]

Hello FriCAS community,

I'm  postdoc studying bioacoustics trying to implement models to predict beam shapes. I'm also a newbie at FriCAS. The acoustic model relies on numerically integrating certain functions and then doing some linear algebra with high precision - and FriCAS seems like a nice platform to run it all fast.

Having searched the docs I found the NumericalQuadrature page. All the routines seem to accept functions which produce a Float. One term in my model produces a complex output, and the routines (romberg, simpson etc.) throw errors.

Is numerical integration with complex functions implemented? Could someone please point me to them if it exists. For reference, the jFricas notebook  with equations is here, and the .input file is attached.

thanks for any help and feedback!

regards,

Thejasvi

https://thejasvibr.github.io/

[Ralf Hemmecke]

Welcome Thejasvi,

I cannot probably give you full help, but I modified your .input file a bit.

Note that ++ is for docstrings in .spad files. The usual comment starts
with --.

The first thing that I have noticed is that legendreP requires its first
argument to be a NonNegativeInteger. Unfortunately, FriCAS does not know
that n-1 for n: PositiveInteger would be of type NonNegativeInteger. You
have to help FriCAS here, i.e.,

(n-1)::NonNegativeInteger

(See start of the notebook. I have defined a few macros that I usually use.)

Still there is a problem with the aromberg function, but that comes from
the fact that its first argument must be of type F->F.

https://fricas.github.io/api/NumericalQuadrature.html#l4e756d65726963616c51756164726174757265-61726f6d62657267

You, however, give a type DoubleFloat -> Complex(DoubleFloat).
That cannot work.

Sorry, I have no time now to analyxe further, but according to

F: FloatingPointSystem

for NumericalQuadrature(F), it will probably not work to use
Complex(DoubleFloat) since

(9) -> Complex(DoubleFloat) has FloatingPointSystem

(9) false

In other words, aromberg is not available for Complex(DF).

Sorry that I cannot be more clear.

Ralf

PS: Please subscribe to fricas-devel for further discussion.

[Ralf Hemmecke]

Oh... attachment...

[Waldek Hebisch]

On Tue, Apr 13, 2021 at 06:28:48PM +0200, thejasvi wrote:
> Having searched the docs I found the NumericalQuadrature

> <https://fricas.github.io/api/NumericalQuadrature.html> page. All the

> routines seem to accept functions which produce a Float. One term in my
> model produces a complex output, and the routines (romberg, simpson etc.)
> throw errors.

Yes, they are only for Float or DoubleFloat

> Is numerical integration with complex functions implemented? Could someone
> please point me to them if it exists.

I am not sure what accuracy you need. For DoubleFloat accuracy
and reasonably regular function routine below should be enough
(and should work faster than romberg or other routines from
NumericalQuadrature).

Note: For irregular functions this routine is somewhat adaptive,
but usualy is doing to much work and gets better accuracy then
required (but my hit limit on recusion level without need)

Example use:

f1(x : DoubleFloat) : Complex(DoubleFloat) == exp(complex(1, x))
ad_gauss(f1, 0, 1, 1.0e-12, 25)

First argument is the function, second is lower limit (lower end
of integration interval), third is upper limit, fourth is bound
on number of recursive calls (bounds above 50 rarely make sense,
25 is probably reasonable default).

--------------<cut here>-----------------------------------------

pl := [0.9681602395_0762608983, 0.8360311073_266357943,
0.6133714327_0059039731, 0.3242534234_0380892904]::List(DoubleFloat)
wl := [0.3302393550_0125976316, 0.0812743883_6157441196_6,
0.1806481606_9485740398, 0.2606106964_029354623,
0.3123470770_4000284007]::List(DoubleFloat)

gauss9(f : DoubleFloat -> Complex(DoubleFloat), x0 : DoubleFloat,
h : DoubleFloat) : Complex(DoubleFloat) ==
h2 := h/(2::DoubleFloat)
xm := x0 + h2
s : Complex(DoubleFloat) := f(xm)*first(wl)
for a in pl for w in rest(wl) repeat
hh := h2*a
s := s + w*(f(xm + hh) + f(xm - hh))
h2*s

ad_gauss(f : DoubleFloat -> Complex(DoubleFloat), x0 : DoubleFloat,
h : DoubleFloat, eps : DoubleFloat, max_level : Integer
) : Complex(DoubleFloat) ==
val0 := gauss9(f, x0, h)
h2 := h/(2::DoubleFloat)
val1 := gauss9(f, x0, h2)
val2 := gauss9(f, x0 + h2, h2)
real(abs(val0 - val1 - val2)) < eps => val0
max_level = 0 =>
print("max_level too small")
val0
eps2 := eps/(2::DoubleFloat)
ad_gauss(f, x0, h2, eps2, max_level - 1) +
ad_gauss(f, x0 + h2, h2, eps2, max_level - 1)

--------------<cut here>-----------------------------
--
Waldek Hebisch

[Thejasvi Beleyur]

Hello Ralf,

Thanks so much for your reply and edits to the code, appreciate it a lot!

Reading through the FriCAS book again about the types, I also realised using DoubleFloat instead of Float was an act of premature optimisation on my part. Either way,

-> Complex(Float) has FloatingPointSystem

also gives

-> false

, which means none of the NumericalQuadrature routines would work as you mentioned.

I'd be happy to hear about any workarounds/alternatives  to integrate a function with complex outputs in FriCAS.

In the meanwhile I'm thinking of seeing if FriCAS can actually symbolically integrate the involved term (even though the textbook says there is no analytical solution).

best,

Thejasvi

[riccard...@gmail.com]

Hi! I would say that the third argument in Waldek's code is upper-lower limit, not the upper limit. Ric

[Waldek Hebisch]

On Tue, Apr 13, 2021 at 08:43:34PM +0200, riccard...@gmail.com wrote:
> Hi! I would say that the third argument in Waldek's code is upper-lower limit, not the upper limit. Ric

Hmm, how do you call b in

\int_a^b f(x)dx

Looking at my post I see that I forgot that routine has five parameters,
fifth (and not fourth) is limit on recursion level, fourth is upper
limit on estimated error.

--
Waldek Hebisch

[riccard...@gmail.com]

@Waldek: I mean, you said that I = \int_a^b f(x)dx = ad_gauss(f, a, b, 1.0e-12, 25), but looking at the use of h (third parameter) in the code I would say that I = ad_gauss(f, a, b-a, 1.0e-12, 25). Ric

[Waldek Hebisch]

On Tue, Apr 13, 2021 at 09:11:55PM +0200, riccard...@gmail.com wrote:
> @Waldek: I mean, you said that I = \int_a^b f(x)dx = ad_gauss(f, a, b, 1.0e-12, 25), but looking at the use of h (third parameter) in the code I would say that I = ad_gauss(f, a, b-a, 1.0e-12, 25). Ric

Right, thanks for correction.

--
Waldek Hebisch

[thejasvi]

Hello Waldek and Ric

Thanks so much for the function and comments you sent.

The integration works pretty fast for the example function you showed.

For my own use-case, I expect to need very high precision (eg. the Mathematica code I'm trying to port uses  300 digits for the same model).

Is DoubleFloat the correct Type to use?

Given that the target is 300 digits precision,  I tried to run the adaptive integration with lower and lower errors (< 1.0e-16).

It ended up taking a much longer time (didn't complete in 3mins). Do any of you have ideas on how to overcome this.

Attached are the current versions of the 'beamshapes.input' and the 'a_gaussianquad.input' which holds the integration functions.

best,

Thejasvi

--
You received this message because you are subscribed to a topic in the Google Groups "FriCAS - computer algebra system" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/fricas-devel/WLoTCKr-I3M/unsubscribe.
To unsubscribe from this group and all its topics, send an email to fricas-devel...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/20210413223424.GD16648%40math.uni.wroc.pl.

[Tobias Neumann]

> For my own use-case, I expect to need very high precision (eg. the Mathematica code I'm trying to port uses  300 digits for the same model).

> Is DoubleFloat the correct Type to use?

> Given that the target is 300 digits precision,  I tried to run the adaptive integration with lower and lower errors (< 1.0e-16).

> It ended up taking a much longer time (didn't complete in 3mins). Do any of you have ideas on how to overcome this.

Are you sure that Mathematica gives you 300 digits precision (PrecisionGoal), and doesn't just require 300 digits working precision (WorkingPrecision)? Do you need 300 digits precision? If yes, then there likely won't be an out-of-the-box solution,

especially with an oscillatory integrand.

*If* Mathematica is able to get you 300 digits (PrecisionGoal) by just using NIntegrate, you should try to figure out which

method it uses. By default it uses heuristics to choose a good method ( https://reference.wolfram.com/language/tutorial/NIntegrateIntegrationStrategies.html ). This is especially relevant

for those oscillatory integrands. I would try to replicate your success with Mathematica by explicitly specifying the integration method with NIntegrate. Once you have done so, you can either find a library that implements that method (or implement it yourself). For example

GSL has some routines for oscillatory integrands 

Best wishes,

Tobias

[Waldek Hebisch]

On Wed, Apr 14, 2021 at 02:50:49PM +0200, thejasvi wrote:
> Hello Waldek and Ric
> Thanks so much for the function and comments you sent.
> The integration works pretty fast for the example function you showed.
>
> For my own use-case, I expect to need very high precision (eg. the
> Mathematica code I'm trying to port uses 300 digits for the same model).
> Is DoubleFloat the correct Type to use?

No, DoubleFloat is limited to about 15 digits. You need to use Float
and use something like 'digits(300)' to set precision. Note: if you
really need 300 digits of final precision, you need to set bigger
precision for intermediate calculation.

>
> Given that the target is 300 digits precision, I tried to run the adaptive
> integration with lower and lower errors (< 1.0e-16).
> It ended up taking a much longer time (didn't complete in 3mins). Do any of
> you have ideas on how to overcome this.

This routine is for DoubleFloat accuracy, using it as-is you can not
get better than what DoubleFloat allows. For better accuracy you
need the following:

- replace DoubleFloat by Float.
- provide pl and wl of sufficient accuracy (that is 300 digits for
300 results)
- probably increase order: current version uses Gaussian kwadrature
on 9 point, which gives order 17. At first glance 75 points
looks like reasonable choice for 300 digits accuracy.

pl are nodes of Gaussian quadrature, they are zeros of appropriate
Legendre polynomial. wl are weights. My code assumes that number
of points is odd, then 0 is one of nodes, other nodes are symmetric.
I took adantage of symmetry to remember only half of nodes and
woights. Below is code that I used to compute nodes and weights
(this one for 45 points):

digits(70)
pn := legendreP(45, x)
sols := solve(%, 1.0e-60)
solp := map(x +-> retract(rhs(x))@Float, last(sols, 23))
dpn := D(pn, x)
wpp := [2/((1 - xi^2)*eval(dpn, x = xi)^2) for xi in solp]
wlf := map(retract, wpp)@List(Float)

plf := rest(solp)

wl := wlf::List(DoubleFloat)
pl := plf::List(DoubleFloat)

Of course, for Float the last two lines should be replaced by
assignments. Note that 'solve' needs substantial time: it
uses robust method which however is slower than Newton when
you need large accuracy. So, using 'solve' with limited
accuracy and then Newton to improve approximation may be
faster. Extra remark: you should compute nodes and weights
once and then use for all integrations with given accuracy.

Extra remarks:
- For highly oscilatory integrals you should use specialized
method. I did not look carefuly enough at your code to know
if you can use general method (Gaussian quadrature) or
need special one
- ATM in FriCAS Bessel functions have only DoubleFloat accuracy.
Bessel functions of half-integer order are elementary, so
for fixed (and hopefuly moderate) order you can use explicit
formulas to provide your own version


--
Waldek Hebisch

[Thejasvi Beleyur]

Hello Tobias,

Thanks a lot for your reply. Apologies for the delay in response.

'Are you sure that Mathematica gives you 300 digits precision (PrecisionGoal), and doesn't just require 300 digits working precision (WorkingPrecision)? Do you need 300 digits precision? ' -- the relevant part of the Mathematica notebook says {$MaxExtraPrecision, prec} = {500, 40} ,  not sure whether this is the actual precision or the 'working precision' as this is the first time I'm handling such involved numerical calculations.

Great idea about digging deeper into how Mathematica itself is dealing with the problem  (I guess you meant Gnu Scientific ilbrary by GSL, this also looks like a cool package . ) - this is a generic approach which I could then use to test how the same algorithm implemented in another package performs.

best,

Thejasvi

[Thejasvi Beleyur]

Hello Waldek,

Thanks a lot for your response. Apologies for the delay in my reply.

I will replace Float and not DoubleFloat and see if I'm still getting working solutions.

I will take a look into modifying my current code to modify it. As Tobias suggested, right now I also perhaps

need to take a deep dive into various methods of numerical integration to get a hang of which one might be

the best for my functions.

best regards,

Thejasvi