Costa,
There are two difficulties in the double numerical integration:
1) The issues around correct quoting etc. of functions.
2) Numerical properties of these integrals when there is an indicator
function.
Issue (1) is easily taken care of with something like what you propose. It
is also taken care of by careful quoting.
But issue (2) is not. As Alexander says, standard numerical integration
techniques are designed for smooth functions -- the step created by an
indicator function "confuses" them. Also, doing a double integration as
the composition of two single integrations is generally going to be far
from optimal both in speed and in accuracy.
I suspect that there are numerical packages (perhaps even in Maxima)
designed specifically for cases like this, but that is beyond my competence.
-s
On Mon, Jan 16, 2012 at 03:08, Constantine Frangos <cfrangos123 at gmail.com>wrote:
> Hi Stavro,
>
> Thanks for the detailed response. I was wondering whether the
> difficutlies you mention
> can be overcome in the following manner.
>
> A new Maxima integration function, num_integration, internally calls a
> user defined Maxima function, say integrandf, and passes the values of
> the integration variables, x1, x2, together with the values of other
> parameters, p1,p2,..., needed to compute the integrand, as follows,
>
> ........
> I1 : apply(integrandf,[x1,x2,p1,p2,...]);
> .........
>
> The user would call num_integration with something like,
>
> I : num_integration(integrandf,list_of_limits,p1,p2,....);
>
> Any suggestions would be appreciated.
>
> Regards,
> Constantine.
>
> On 1/12/12, Stavros Macrakis <macrakis at alum.mit.edu> wrote:
> > As Alexander says, you don't need the temporary variable etc.
> >
> > Also, if you explicitly quote everything, you are less likely to have bad
> > surprises. Something like this:
> >
> > (%i27) (cnt:0,quad_qags('(quad_qags('(cnt:cnt+1,if x>y then 1 else
> > 0),x,0,1)[1]),y,0,1));
> > (%o27) [0.49999763949292, 3.4177098928722671E-9, 735, 0]
> > (%i28) cnt;
> > (%o28) 621243
> >
> > Note the number of function evaluations for a very simple case! Note
> also
> > that the result is not very accurate, and that the error estimate is way
> > off. I imagine there are better routines for this kind of integration.
> > Obviously in this case it is easy to transform to a single integral
> (even
> > analytical), but presumably you want to handle more complicated
> integrands.
> >
> > -s
> >
> > On Thu, Jan 12, 2012 at 05:38, Alexander Klimov <alserkli at inbox.ru>
> wrote:
> >
> >> Hi.
> >>
> >> On Thu, 12 Jan 2012, Constantine Frangos wrote:
> >> > (1) Its not clear why the function definition I am using does not
> >> > work. According to the Maxima 5.24 manual, paragraph 20.4, the
> >> > integrand can be the name of a Maxima function, presumably defined by
> >> > the user in a file, for example, integrand.mac
> >>
> >> Your function does not always return a number. Consider your original
> >> definition:
> >>
> >> (%i1) ind(x,y) := block([s],if is(x*x + y*y <= 1) then (s : 1) else (s :
> >> 0),return(s))$
> >>
> >> Let us see how it works
> >>
> >> (%i2) trace(ind);
> >> (%o2) [ind]
> >>
> >> If one gives numerical arguments, then your function returns the
> >> number which is assigned to "s"
> >>
> >> (%i3) ind(0,1);
> >> 1 Enter ind [0, 1]
> >> 1 Exit ind 1
> >> (%o3) 1
> >>
> >> but with symbols as arguments nothing is assigned to "s" and thus the
> >> function returns just a symbol "s"
> >>
> >> (%i4) ind(x1,x2);
> >> 1 Enter ind [x1, x2]
> >> 1 Exit ind s
> >> (%o4) s
> >>
> >> and this is exactly what happens when you use it in integration
> >>
> >> (%i5) I:quad_qags(quad_qags(ind(x1,x2),x1,-1,1)[1],x2,-1,1);
> >> 1 Enter ind [x1, x2]
> >> 1 Exit ind s
> >> COERCE-FLOAT-FUN: no such Lisp or Maxima function: s
> >> -- an error. To debug this try: debugmode(true);
> >>
> >> > (2) Are there perhaps other Maxima numerical integration functions
> >> > that may perform better ?
> >>
> >> As I said, it is better to
> >>
> >> > > use a simple integral:
> >> > >
> >> > > (%i4) quad_qags(2*sqrt(1-y^2),y,-1,1);
> >> > > (%o4) [3.141592653589797, 2.000470900043183e-9, 399, 0]
> >>
> >> the problem with indicators is that while they may be suitable for
> >> symbolic manipulations, they are bad functions for Gaussian quadrature
> >> rule which is only good if the function is similar to a polynomial.
> >>
> >> --
> >> Regards,
> >> ASK
> >> _______________________________________________
> >> Maxima mailing list
> >> Maxima at math.utexas.edu
> >> http://www.math.utexas.edu/mailman/listinfo/maxima
> >>
> >
>