[GiNaC-devel] symbolic integration.
Chris Dams
C.Dams at science.ru.nl
Sun Sep 26 14:50:04 CEST 2004
Dear all,
Here's a patch plus two new files for GiNaC that add a class to
symbolically represent an integral. Details can be found in the patch,
since it also patches the documentation.
Best,
Chris
-------------- next part --------------
Index: doc/tutorial/ginac.texi
===================================================================
RCS file: /home/cvs/GiNaC/doc/tutorial/ginac.texi,v
retrieving revision 1.155
diff -c -r1.155 ginac.texi
*** doc/tutorial/ginac.texi 20 Aug 2004 16:14:10 -0000 1.155
--- doc/tutorial/ginac.texi 24 Sep 2004 11:36:26 -0000
***************
*** 693,698 ****
--- 693,699 ----
* Lists:: Lists of expressions.
* Mathematical functions:: Mathematical functions.
* Relations:: Equality, Inequality and all that.
+ * Integrals:: Symbolic integrals.
* Matrices:: Matrices.
* Indexed objects:: Handling indexed quantities.
* Non-commutative objects:: Algebras with non-commutative products.
***************
*** 1807,1813 ****
wrapped inside an @code{ex}.
! @node Relations, Matrices, Mathematical functions, Basic Concepts
@c node-name, next, previous, up
@section Relations
@cindex @code{relational} (class)
--- 1808,1814 ----
wrapped inside an @code{ex}.
! @node Relations, Integrals, Mathematical functions, Basic Concepts
@c node-name, next, previous, up
@section Relations
@cindex @code{relational} (class)
***************
*** 1833,1840 ****
however, that @code{==} here does not perform any simplifications, hence
@code{expand()} must be called explicitly.
! @node Matrices, Indexed objects, Relations, Basic Concepts
@c node-name, next, previous, up
@section Matrices
@cindex @code{matrix} (class)
--- 1834,1901 ----
however, that @code{==} here does not perform any simplifications, hence
@code{expand()} must be called explicitly.
+ @node Integrals, Matrices, Relations, Basic Concepts
+ @c node-name, next, previous, up
+ @section Integrals
+ @cindex @code{integral} (class)
+
+ An object of class @dfn{integral} can be used to hold a symbolic integral.
+ If you want to symbolically represent the integral of @code{x*x} from 0 to
+ 1, you would write this as
+ @example
+ integral(x, 0, 1, x*x)
+ @end example
+ The first argument is the integration variable. It should be noted that
+ GiNaC is not very good (yet?) at symbolically evaluating integrals. In
+ fact, it can only integrate polynomials. An expression containing integrals
+ can be evaluated symbolically by calling the
+ @example
+ .eval_integ()
+ @end example
+ method on it. Numerical evaluation is available by calling the
+ @example
+ .evalf()
+ @end example
+ method on an expression containing the integral. This will only evaluate
+ integrals into a number if @code{subs}ing the integration variable by a
+ number in the fourth argument of an integral and then @code{evalf}ing the
+ result always results in a number. Of course, also the boundaries of the
+ integration domain must @code{evalf} into numbers. It should be noted that
+ trying to @code{evalf} a function with discontinuities in the integration
+ domain is not recommended. The accuracy of the numeric evaluation of
+ integrals is determined by the static member variable
+ @example
+ ex integral::relative_integration_error
+ @end example
+ of the class @code{integral}. The default value of this is 10^-8.
+ The integration works by halving the interval of integration, until numeric
+ stability of the answer indicates that the requested accuracy has been
+ reached. The maximum depth of the halving can be set via the static member
+ variable
+ @example
+ int integral::max_integration_level
+ @end example
+ The default value is 15. If this depth is exceeded, @code{evalf} will simply
+ return the integral unevaluated. The function that performs the numerical
+ evaluation, is also available as
+ @example
+ ex adaptivesimpson(const ex&x, const ex &a, const ex &b, const ex &f,
+ const ex &error)
+ @end example
+ This function will throw an exception if the maximum depth is exceeded. The
+ last parameter of the function is optional and defaults to the
+ @code{relative_integration_error}. To make sure that we do not do too
+ much work if an expression contains the same integral multiple times,
+ a lookup table is used.
+
+ If you know that an expression holds an integral, you can get the
+ integration variable, the left boundary, right boundary and integrant by
+ respectively calling @code{.op(0)}, @code{.op(1)}, @code{.op(2)}, and
+ @code{.op(3)}. Differentiating integrals with respect to variables works
+ as expected. Note that it makes no sense to differentiate an integral
+ with respect to the integration variable.
! @node Matrices, Indexed objects, Integrals, Basic Concepts
@c node-name, next, previous, up
@section Matrices
@cindex @code{matrix} (class)
Index: ginac/Makefile.am
===================================================================
RCS file: /home/cvs/GiNaC/ginac/Makefile.am,v
retrieving revision 1.32
diff -c -r1.32 Makefile.am
*** ginac/Makefile.am 18 Dec 2003 22:28:01 -0000 1.32
--- ginac/Makefile.am 24 Sep 2004 11:36:26 -0000
***************
*** 5,14 ****
constant.cpp ex.cpp expair.cpp expairseq.cpp exprseq.cpp fail.cpp \
fderivative.cpp function.cpp idx.cpp indexed.cpp inifcns.cpp \
inifcns_trans.cpp inifcns_gamma.cpp inifcns_nstdsums.cpp \
! lst.cpp matrix.cpp mul.cpp ncmul.cpp normal.cpp numeric.cpp operators.cpp \
! power.cpp registrar.cpp relational.cpp remember.cpp pseries.cpp print.cpp \
! structure.cpp symbol.cpp symmetry.cpp tensor.cpp utils.cpp wildcard.cpp \
! input_parser.yy input_lexer.ll \
input_lexer.h remember.h tostring.h utils.h
libginac_la_LDFLAGS = -version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE) \
-release $(LT_RELEASE)
--- 5,14 ----
constant.cpp ex.cpp expair.cpp expairseq.cpp exprseq.cpp fail.cpp \
fderivative.cpp function.cpp idx.cpp indexed.cpp inifcns.cpp \
inifcns_trans.cpp inifcns_gamma.cpp inifcns_nstdsums.cpp \
! integral.cpp lst.cpp matrix.cpp mul.cpp ncmul.cpp normal.cpp numeric.cpp \
! operators.cpp power.cpp registrar.cpp relational.cpp remember.cpp \
! pseries.cpp print.cpp structure.cpp symbol.cpp symmetry.cpp tensor.cpp \
! utils.cpp wildcard.cpp input_parser.yy input_lexer.ll \
input_lexer.h remember.h tostring.h utils.h
libginac_la_LDFLAGS = -version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE) \
-release $(LT_RELEASE)
***************
*** 16,24 ****
ginacinclude_HEADERS = ginac.h add.h archive.h assertion.h basic.h class_info.h \
clifford.h color.h constant.h container.h ex.h expair.h expairseq.h exprseq.h \
fail.h fderivative.h flags.h function.h hash_map.h idx.h indexed.h inifcns.h \
! lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h power.h print.h \
! pseries.h ptr.h registrar.h relational.h structure.h symbol.h symmetry.h \
! tensor.h tinfos.h version.h wildcard.h
AM_LFLAGS = -Pginac_yy -olex.yy.c
AM_YFLAGS = -p ginac_yy -d
EXTRA_DIST = function.pl input_parser.h version.h.in
--- 16,24 ----
ginacinclude_HEADERS = ginac.h add.h archive.h assertion.h basic.h class_info.h \
clifford.h color.h constant.h container.h ex.h expair.h expairseq.h exprseq.h \
fail.h fderivative.h flags.h function.h hash_map.h idx.h indexed.h inifcns.h \
! integral.h lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h \
! power.h print.h pseries.h ptr.h registrar.h relational.h structure.h \
! symbol.h symmetry.h tensor.h tinfos.h version.h wildcard.h
AM_LFLAGS = -Pginac_yy -olex.yy.c
AM_YFLAGS = -p ginac_yy -d
EXTRA_DIST = function.pl input_parser.h version.h.in
Index: ginac/basic.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/basic.cpp,v
retrieving revision 1.87
diff -c -r1.87 basic.cpp
*** ginac/basic.cpp 3 Aug 2004 15:20:37 -0000 1.87
--- ginac/basic.cpp 24 Sep 2004 11:36:26 -0000
***************
*** 207,213 ****
* @see basic::dbgprinttree */
void basic::dbgprint() const
{
! this->print(std::cerr);
std::cerr << std::endl;
}
--- 207,213 ----
* @see basic::dbgprinttree */
void basic::dbgprint() const
{
! this->print(print_dflt(std::cerr));
std::cerr << std::endl;
}
***************
*** 480,485 ****
--- 480,499 ----
return *this;
else
return map(map_evalm);
+ }
+
+ /** Function object to be applied by basic::eval_integ(). */
+ struct eval_integ_map_function : public map_function {
+ ex operator()(const ex & e) { return eval_integ(e); }
+ } map_eval_integ;
+
+ /** Evaluate integrals, if result is known. */
+ ex basic::eval_integ() const
+ {
+ if(nops() == 0)
+ return *this;
+ else
+ return map(map_eval_integ);
}
/** Perform automatic symbolic evaluations on indexed expression that
Index: ginac/basic.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/basic.h,v
retrieving revision 1.69
diff -c -r1.69 basic.h
*** ginac/basic.h 25 May 2004 13:59:32 -0000 1.69
--- ginac/basic.h 24 Sep 2004 11:36:26 -0000
***************
*** 136,141 ****
--- 136,142 ----
virtual ex eval(int level = 0) const;
virtual ex evalf(int level = 0) const;
virtual ex evalm() const;
+ virtual ex eval_integ() const;
protected:
virtual ex eval_ncmul(const exvector & v) const;
public:
Index: ginac/ex.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/ex.h,v
retrieving revision 1.78
diff -c -r1.78 ex.h
*** ginac/ex.h 20 Aug 2004 16:14:12 -0000 1.78
--- ginac/ex.h 24 Sep 2004 11:36:26 -0000
***************
*** 117,122 ****
--- 117,123 ----
ex evalf(int level = 0) const { return bp->evalf(level); }
ex evalm() const { return bp->evalm(); }
ex eval_ncmul(const exvector & v) const { return bp->eval_ncmul(v); }
+ ex eval_integ() const { return bp->eval_integ(); }
// printing
void print(const print_context & c, unsigned level = 0) const;
***************
*** 806,811 ****
--- 807,815 ----
inline ex evalm(const ex & thisex)
{ return thisex.evalm(); }
+
+ inline ex eval_integ(const ex & thisex)
+ { return thisex.eval_integ(); }
inline ex diff(const ex & thisex, const symbol & s, unsigned nth = 1)
{ return thisex.diff(s, nth); }
Index: ginac/ginac.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/ginac.h,v
retrieving revision 1.22
diff -c -r1.22 ginac.h
*** ginac/ginac.h 8 Jan 2004 15:06:48 -0000 1.22
--- ginac/ginac.h 24 Sep 2004 11:36:26 -0000
***************
*** 34,39 ****
--- 34,40 ----
#include "constant.h"
#include "fail.h"
+ #include "integral.h"
#include "lst.h"
#include "matrix.h"
#include "numeric.h"
Index: ginac/indexed.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/indexed.cpp,v
retrieving revision 1.92
diff -c -r1.92 indexed.cpp
*** ginac/indexed.cpp 5 Aug 2004 17:06:16 -0000 1.92
--- ginac/indexed.cpp 24 Sep 2004 11:36:26 -0000
***************
*** 36,41 ****
--- 36,42 ----
#include "lst.h"
#include "archive.h"
#include "utils.h"
+ #include "integral.h"
namespace GiNaC {
***************
*** 520,525 ****
--- 521,533 ----
return really_free;
} else
return basis_indices;
+ }
+
+ exvector integral::get_free_indices() const
+ {
+ if(a.get_free_indices().size() || b.get_free_indices().size())
+ throw (std::runtime_error("integral::get_free_indices: boundary values should not have free indices"));
+ return f.get_free_indices();
}
/** Rename dummy indices in an expression.
Index: ginac/tinfos.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/tinfos.h,v
retrieving revision 1.26
diff -c -r1.26 tinfos.h
*** ginac/tinfos.h 29 Apr 2004 17:25:29 -0000 1.26
--- ginac/tinfos.h 24 Sep 2004 11:36:26 -0000
***************
*** 38,43 ****
--- 38,44 ----
const unsigned TINFO_function = 0x00031001U;
const unsigned TINFO_fderivative = 0x00032001U;
const unsigned TINFO_ncmul = 0x00031002U;
+ const unsigned TINFO_integral = 0x00033001U;
const unsigned TINFO_lst = 0x00040001U;
Index: ginac/wildcard.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/wildcard.cpp,v
retrieving revision 1.17
diff -c -r1.17 wildcard.cpp
*** ginac/wildcard.cpp 8 Jan 2004 15:06:53 -0000 1.17
--- ginac/wildcard.cpp 24 Sep 2004 11:36:26 -0000
***************
*** 119,122 ****
--- 119,132 ----
return is_equal(ex_to<basic>(pattern));
}
+ bool haswild(const ex&x)
+ {
+ if(is_a<wildcard>(x))
+ return true;
+ for(int i=0; i<x.nops(); ++i)
+ if(haswild(x.op(i)))
+ return true;
+ return false;
+ }
+
} // namespace GiNaC
Index: ginac/wildcard.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/wildcard.h,v
retrieving revision 1.12
diff -c -r1.12 wildcard.h
*** ginac/wildcard.h 8 Jan 2004 15:06:54 -0000 1.12
--- ginac/wildcard.h 24 Sep 2004 11:36:26 -0000
***************
*** 75,80 ****
--- 75,83 ----
return wildcard(label);
}
+ /** Check whether x has a wildcard anywhere as a subexpression. */
+ bool haswild(const ex&x);
+
} // namespace GiNaC
#endif // ndef __GINAC_WILDCARD_H__
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: integral.h
Url: http://www.cebix.net/pipermail/ginac-devel/attachments/20040924/4f229ddd/integral.h
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: integral.cpp
Url: http://www.cebix.net/pipermail/ginac-devel/attachments/20040924/4f229ddd/integral.cpp
More information about the GiNaC-devel
mailing list