[fricas-devel] Simple Gaussian integral fails

[Marduk]

Dear all,

I just tried calculating this integral (I defined inf ==> %plusInfinity)

integrate(exp(-a^2*x^2),x=-inf..inf)

in FriCAS 1.3.4 and it failed. Why?

Best regards,
Marduk

[Ralf Hemmecke]

I suspect, the problem lies in the limit computation.

Ralf

(1) -> i := integrate(exp(-a^2*x^2),x)::Expression(Integer)

+--+
+---+ | 2
\|%pi erf(x\|a )
(1) -----------------
+--+
| 2
2 \|a
Type:
Expression(Integer)
(2) -> limit(erf(x),x=%plusInfinity)

(2) 1
Type:
Union(OrderedCompletion(Expression(Integer)),...)
(3) -> limit(erf(sqrt(a^2)*x),x=%plusInfinity)

(3) "failed"
Type:
Union("failed",...)

Post by Marduk
Dear all,
I just tried calculating this integral (I defined inf ==> %plusInfinity)
integrate(exp(-a^2*x^2),x=-inf..inf)
in FriCAS 1.3.4 and it failed. Why?
Best regards,
Marduk

[Marduk]

Interesting. I did some experiments and it seems that limit(erf(M*x),
x=inf) returns 1
whenever FriCAS can determine that M > 0. In particular M=a^2 and M=exp(a)
work,
whereas M=a does not. However, M = abs(a) does not work either.

Now I have some questions:

1) Why does not sqrt(a^2) return abs(a)?

2) Why does not FriCAS know that abs(a) > 0?

3) How can one tell FriCAS that a>0 in the Gaussian integral?
In Maxima this is done with assume. Maybe in FriCAS one could use
predicates.

[Ralf Hemmecke]

Post by Marduk
1) Why does not sqrt(a^2) return abs(a)?

OK, let's try with a being the imaginary unity.

(4) -> sqrt(%i^2)

(4) %i
Type:
Expression(Complex(Integer))
(5) -> abs(%i)

(5) abs(%i)

FriCAS does not simplify abs(%i). So it's left unspecified.
You would probably not agree that abs(%i)=%i

Post by Marduk
2) Why does not FriCAS know that abs(a) > 0?

Because it is just an operator and a can be anything. What would you
want abs("Hello World!") to return?

Post by Marduk
3) How can one tell FriCAS that a>0 in the Gaussian integral?
In Maxima this is done with assume. Maybe in FriCAS one could use
predicates.

I've no idea.

Ralf

[Marduk]

Post by Ralf Hemmecke
OK, let's try with a being the imaginary unity.

I thought that by default FriCAS assumes that variables are real.

Post by Ralf Hemmecke
Because it is just an operator and a can be anything.

Fair enough. Then one should be able to specify on which domain an operator
is acting
and abs(a) | a is real should know that a > 0.

I just found out that the exact same question was asked 7 years ago:
https://groups.google.com/d/topic/fricas-devel/jTeaKzcEJgs/discussion

It is surprising that nothing has been done to correct this problem. But
reading that thread
explains why. It is not a matter of mathematics or computer science, but a
matter of sociology.

To us physicists, Juanjo in that thread and myself, a CAS should be a tool
that helps us save
time by performing the calculations that come up in our work. But to the
current maintainers of
the big three: Maxima, FriCAS and Reduce, it does not matter that their
programs cannot calculate
Gaussian integrals (all three have failed my tests).

In the Axiom/FriCAS wiki it is claimed that this CAS boasts the most
complete implementation of the
Risch algorithm. And I am sure that it can calculate thousands of
integrals. But what is the purpose of
having a sophisticated integrator if it cannot calculate the simplest
integrals? In the thread above
Waldek replies: do not use it for calculating integrals you already know.
Seriously?

I must say that it is highly disappointing that to this date there is no
serious free-software alternative to
Mathematica. The three programs mentioned above have been developed for
about 40-50 years by very
capable computer scientists and mathematicians. Axiom/FriCAS is IMHO the
best designed of the three,
it is thoroughly documented and all that effort is wasted because the
developers do not care for the needs
of the users. Pretty sad.

Marduk

[Ralf Hemmecke]

Post by Marduk
I must say that it is highly disappointing that to this date there
is no serious free-software alternative to Mathematica.

It all depends on what you want to do.

Post by Marduk
The three programs mentioned above have been developed for about
40-50 years by very capable computer scientists and mathematicians.
Axiom/FriCAS is IMHO the best designed of the three, it is thoroughly
documented and all that effort is wasted because the developers do
not care for the needs of the users. Pretty sad.

Surely this is sad. But in contrast to Mathematica developers, FriCAS
developers don't see a single cent for their efforts. People do that in
their spare time. It's just natural that solving their own problems
comes first. If that helps other people then fine.

Clearly, we could setup a channel where people would be able to devote
some money to FriCAS developers in order to get a certain problem solved
(implemented).

It's not an excuse, but also open source development costs time and
effort. The sad thing is that FriCAS does not have enough developers,
but a lot of places where it needs improvements.

So everyone is welcome to contribute, even if it is just bug reports.

Success to whatever you want to achieve.

Ralf

[Riccardo GUIDA]

Post by Marduk
3) How can one tell FriCAS that a>0 in the Gaussian integral?

Write a=A^2 :

(1) -> inf ==> %plusInfinity; integrate(exp(-A^4*x^2),x=-inf..inf)

┌───┐
\│%pi
(1) ──────
2
A
Type: Union(f1: OrderedCompletion(Expression(Integer)),...)


another sad physicist...

[Riccardo GUIDA]

Post by Marduk
1) Why does not sqrt(a^2) return abs(a)?

IIUC:

When operating on a symbolic expression x in Expression R,
the operator sqrt "represents an arbitrary (but always the same through the code) solution y of the algebraic equation y^2=x".
Even when R = Integer, there are still 2 real roots y that satisfy the equation.

Similar semantics for the n-th root nthRoot:
"an arbitrary (but always the same through the code) solution of the algebraic equation y^n=x".

Note that sqrt(x) may not be imagined to be a set with two elements
because the (pointwise) product of two sets would be a set including also
non-diagonal products of the elements.
(IE sqrt(4)={+2,-2} ==> sqrt(4)*sqrt(4) = {(+2)*(+2),...} = {+4, -4}).

So, FriCAS deliberately tries not to make a choice of roots
and consequently sqrt has a different semantics (on expressions)
as compared to all other CAS in my knowledge.

As for the (sometimes) weird behavior of abs in FriCAS
this should be related to the fact that using a smart abs one might construct zero divisors,
(nonzero x1,x2 such that x1 * x2 = 0). This would violate the fact that Expression Integer
is meant to be a differential field, which is essential for the indefinite integration algorithms.
To recover more information I advice you to track the above keywords in the past messages:
many of your (present and future) questions have been already asked many times in the last 10 years.

riccardo

[oldk1331]

I think it's doable.

First, by default, "limit" is for real expression (there are
accompanying "complexLimit"), so it knows "sqrt(a)" is positive:

(8) -> limit(sqrt(a)*x,x=%plusInfinity)

(8) + infinity

So in theory, it should also compute for "erf(sqrt(a)*x)", but:

(9) -> limit(erf(sqrt(a)*x),x=%plusInfinity)

(9) "failed"

This should be very simple to compute, but we still use Gruntz
algorithm to compute and we fail at there. This may indicate
a bug in Gruntz algorithm, and raise another question: shall
we add a short path to handle simple cases, like "specialLimit"
in limitps.spad?

Another somewhat related problem:

(1) -> limit(sqrt(a^2)*x,x=%plusInfinity)

(1) "failed"

This case fails because "mrv_normalize" will destroy the structure
of "sqrt(a^2)", so maybe a "assume" system will still be needed.
But such a system is invasive and not always work in Maxima.

Post by Marduk
Dear all,
I just tried calculating this integral (I defined inf ==> %plusInfinity)
integrate(exp(-a^2*x^2),x=-inf..inf)
in FriCAS 1.3.4 and it failed. Why?
Best regards,
Marduk
--
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
Visit this group at https://groups.google.com/group/fricas-devel.
For more options, visit https://groups.google.com/d/optout.

[Waldek Hebisch]

Post by Marduk

Post by Ralf Hemmecke
OK, let's try with a being the imaginary unity.

I thought that by default FriCAS assumes that variables are real.

Users requently try to use complex expressions and when possible
FriCAS tries to accomodate this.

Post by Marduk
https://groups.google.com/d/topic/fricas-devel/jTeaKzcEJgs/discussion
It is surprising that nothing has been done to correct this problem. But
reading that thread
explains why. It is not a matter of mathematics or computer science, but a
matter of sociology.

Well, concerning "nothing has been done" it depends what "this problem"
is. If you mean exact expression appearing in the report, it is
just one of several thousends of possible inputs. If you mean problem
Juanjo raised (better support of special functions), then a lot
happened. In particular FriCAS is now quite good at _indefinite_
integrals with result involving 'erf'. Concerning sociology,
any specific integral could be added via pattern matching (FriCAS
has special pattern matcher for definite integrals), that is
about 40 lines of code. That is great for benchmarketing
(you can show 100% success rate on any small testsuite).
However, if you want comprehensive support in that way you would
need to add hundreds or thousends of integrals meaning several
thoushends of lines of code. Tedious and still it would be
easy to find integrals which are not handled. I prefer
algorithms that handle _all_ functions in given relatively
large class. To put is differently, if I thought that handling
of this specific Gaussian integral is most imporant feature
for FriCAS I would implement it. But fixing via special
cases several such "most imporant features" would soak
all time that I can spent on FriCAS developement and IMO
would result in much weaker system compared to current FriCAS.

Post by Marduk
To us physicists, Juanjo in that thread and myself, a CAS should be a tool
that helps us save
time by performing the calculations that come up in our work. But to the
current maintainers of
the big three: Maxima, FriCAS and Reduce, it does not matter that their
programs cannot calculate
Gaussian integrals (all three have failed my tests).

World is not black and white. You write "Gaussian integrals" in
the plural here, so maybe you mean things like (indefinite)

integrate(x^16*exp(-x^2), x)

In the past FriCAS could not do them, now it is handled by
general routine. So it matters. But fixing _one_ definite
integral is not very significant when there are thousends
of similar cases. Now, solutions that handle large classes
take time to develop. IMO seeking general solutions pays
in longer time.

Extra remark about "does not matter": if I would serioulsy
worry about every similar problem/bug I would quickly
end up in a madhouse with serious mental ilness. And
the same applies to other sofware. Mathematica and Maple
are very secretive about their problems, but there are
enough reports to know that they have a lot of problems.

Post by Marduk
In the Axiom/FriCAS wiki it is claimed that this CAS boasts the most
complete implementation of the
Risch algorithm. And I am sure that it can calculate thousands of
integrals. But what is the purpose of
having a sophisticated integrator if it cannot calculate the simplest
integrals? In the thread above
Waldek replies: do not use it for calculating integrals you already know.
Seriously?

Read carefuly, it is not what I wrote. Hint: I was discussing possible
_pragmatic_ approches given current FriCAS limitations.

And in this disscussion you were given IMHO reasonable workaround (using
fourth power of parameter), so this problem is not a "killer" in
sense that stops you from doing serious work. Annoying, yes,
so I understand why you are annoyed...

[Waldek Hebisch]

Post by oldk1331
I think it's doable.
First, by default, "limit" is for real expression (there are
(8) -> limit(sqrt(a)*x,x=%plusInfinity)
(8) + infinity
(9) -> limit(erf(sqrt(a)*x),x=%plusInfinity)
(9) "failed"
This should be very simple to compute, but we still use Gruntz
algorithm to compute and we fail at there. This may indicate
a bug in Gruntz algorithm, and raise another question: shall
we add a short path to handle simple cases, like "specialLimit"
in limitps.spad?

Gruntz seem to be close to computing this limit, so I do not
see special reason to add alternate code path. More generaly,
various "short paths" or "special paths" may seem attractive
at first, but they are pain in longer time.

Post by oldk1331
(1) -> limit(sqrt(a^2)*x,x=%plusInfinity)
(1) "failed"
This case fails because "mrv_normalize" will destroy the structure
of "sqrt(a^2)",

Do not attach too much importance to 'sqrt' in integral that
started this discussion:

+---+
erf(a x)\|%pi
(4) --------------
2 a
Type: Expression(Integer)

is a valid indefinite integral and we should be generating it
instead of version with 'sqrt(a^2)'. Note that for this
expression limit splits into two cases: positive and negative
and _both_ give the same limit. So we would need version of
Gruntz that can split and combine cases.

And for treating 'sqrt(a^2)' as 'abs(x)' we need version of
'normalize' that splits cases.

Post by oldk1331
so maybe a "assume" system will still be needed.
But such a system is invasive and not always work in Maxima.

Yes, in longer term we will need "assume" and handling of
conditions.