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
--
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/fricas-devel/20210413182353.GA40312%40math.uni.wroc.pl.
