[GiNaC-devel] patch for series expansion of integral.
Chris Dams
C.Dams at science.ru.nl
Thu Oct 28 14:10:04 CEST 2004
Dear GiNaCers,
Here's a patch. Find out below why you want this patch ;-).
The integral class has a .derivative() method, so in principle integrals
can be series expanded. However, if somewhere deep down in the integrant a
sum with a lot of terms occurs, the same kind of problems that have
troubled the series expansion of powers, happen again. Basically we get
expression explosion. To solve this I wrote an integral::series() method.
While doing this I discovered several bugs/undesirable features, namely
(1) When requesting terms of a pseries via pseries::op(), the order term
looks like O(1)^n instead of O(x^n). This looks strange and different from
the output that one gets from printing the entire series.
(2) Since pseries does not have a let_op() method for understandable
reasons, eval_integ() does not work for a power series. I have added a new
method pseries::eval_integ() to solve this.
(3) When trying to find coefficients and exponents of power series, using
pseries::op() is not convenient and not efficient. This returns a
multiplication of the coefficient with a power that the user then may
attempt to extract the coefficient and the exponent from... To make this
easier, I added pseries::coeffop() and pseries::exponop() methods.
(4) pseries::mul_series() does not know that zero times a power series is
zero.
(5) If pseries::power_const() discovers that the user is expanding to such
a low degree that only the order term needs to be returned, GiNaC
sometimes tries to save memory by allocating a negative amount of bytes in
order to store a negative number of terms in it.
Best,
Chris
-------------- next part --------------
Index: ginac/integral.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/integral.h,v
retrieving revision 1.1
diff -c -r1.1 integral.h
*** ginac/integral.h 12 Oct 2004 09:32:11 -0000 1.1
--- ginac/integral.h 27 Oct 2004 07:57:04 -0000
***************
*** 56,61 ****
--- 56,62 ----
ex eval_integ() const;
protected:
ex derivative(const symbol & s) const;
+ ex series(const relational & r, int order, unsigned options = 0) const;
// new virtual functions which can be overridden by derived classes
// none
Index: ginac/pseries.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/pseries.cpp,v
retrieving revision 1.82
diff -c -r1.82 pseries.cpp
*** ginac/pseries.cpp 6 Oct 2004 11:59:22 -0000 1.82
--- ginac/pseries.cpp 27 Oct 2004 07:57:04 -0000
***************
*** 33,38 ****
--- 33,39 ----
#include "relational.h"
#include "operators.h"
#include "symbol.h"
+ #include "integral.h"
#include "archive.h"
#include "utils.h"
***************
*** 266,271 ****
--- 267,274 ----
if (i >= seq.size())
throw (std::out_of_range("op() out of range"));
+ if(is_order_function(seq[i].rest))
+ return Order(power(var-point, seq[i].coeff));
return seq[i].rest * power(var - point, seq[i].coeff);
}
***************
*** 424,429 ****
--- 427,457 ----
return result;
}
+ ex pseries::eval_integ() const
+ {
+ epvector *newseq = 0;
+ for(epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+ if(newseq)
+ { newseq->push_back(expair(i->rest.eval_integ(), i->coeff));
+ continue;
+ }
+ ex newterm = i->rest.eval_integ();
+ if(!are_ex_trivially_equal(newterm, i->rest)) {
+ newseq = new epvector;
+ newseq->reserve(seq.size());
+ for(epvector::const_iterator j=seq.begin(); j!=i; ++j)
+ newseq->push_back(*j);
+ newseq->push_back(expair(newterm, i->coeff));
+ }
+ }
+
+ ex newpoint = point.eval_integ();
+ if(newseq || !are_ex_trivially_equal(newpoint, point))
+ return (new pseries(var==newpoint, *newseq))
+ ->setflag(status_flags::dynallocated);
+ return *this;
+ }
+
ex pseries::subs(const exmap & m, unsigned options) const
{
// If expansion variable is being substituted, convert the series to a
***************
*** 519,524 ****
--- 547,566 ----
return seq.empty() || !is_order_function((seq.end()-1)->rest);
}
+ ex pseries::coeffop(size_t i) const
+ {
+ if(i > nops())
+ throw (std::out_of_range("coeffop() out of range"));
+ return seq[i].rest;
+ }
+
+ ex pseries::exponop(size_t i) const
+ {
+ if(i > nops())
+ throw (std::out_of_range("exponop() out of range"));
+ return seq[i].coeff;
+ }
+
/*
* Implementations of series expansion
***************
*** 729,734 ****
--- 771,781 ----
nul.push_back(expair(Order(_ex1), _ex0));
return pseries(relational(var,point), nul);
}
+
+ if(seq.empty() || other.seq.empty()) {
+ return (new pseries(var==point, epvector()))
+ ->setflag(status_flags::dynallocated);
+ }
// Series multiplication
epvector new_seq;
***************
*** 880,886 ****
throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
// adjust number of coefficients
! deg = deg - (p*ldeg).to_int();
// O(x^n)^(-m) is undefined
if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative())
--- 927,940 ----
throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
// adjust number of coefficients
! int numcoeff = deg - (p*ldeg).to_int();
! if(numcoeff <= 0) {
! epvector epv;
! epv.reserve(1);
! epv.push_back(expair(Order(_ex1), deg));
! return (new pseries(relational(var,point), epv))
! ->setflag(status_flags::dynallocated);
! }
// O(x^n)^(-m) is undefined
if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative())
***************
*** 888,896 ****
// Compute coefficients of the powered series
exvector co;
! co.reserve(deg);
co.push_back(power(coeff(var, ldeg), p));
! for (int i=1; i<deg; ++i) {
ex sum = _ex0;
for (int j=1; j<=i; ++j) {
ex c = coeff(var, j + ldeg);
--- 942,950 ----
// Compute coefficients of the powered series
exvector co;
! co.reserve(numcoeff);
co.push_back(power(coeff(var, ldeg), p));
! for (int i=1; i<numcoeff; ++i) {
ex sum = _ex0;
for (int j=1; j<=i; ++j) {
ex c = coeff(var, j + ldeg);
***************
*** 906,912 ****
// Construct new series (of non-zero coefficients)
epvector new_seq;
bool higher_order = false;
! for (int i=0; i<deg; ++i) {
if (!co[i].is_zero())
new_seq.push_back(expair(co[i], p * ldeg + i));
if (is_order_function(co[i])) {
--- 960,966 ----
// Construct new series (of non-zero coefficients)
epvector new_seq;
bool higher_order = false;
! for (int i=0; i<numcoeff; ++i) {
if (!co[i].is_zero())
new_seq.push_back(expair(co[i], p * ldeg + i));
if (is_order_function(co[i])) {
***************
*** 915,921 ****
}
}
if (!higher_order)
! new_seq.push_back(expair(Order(_ex1), p * ldeg + deg));
return pseries(relational(var,point), new_seq);
}
--- 969,975 ----
}
}
if (!higher_order)
! new_seq.push_back(expair(Order(_ex1), p * ldeg + numcoeff));
return pseries(relational(var,point), new_seq);
}
***************
*** 1033,1038 ****
--- 1087,1149 ----
}
} else
return convert_to_poly().series(r, order, options);
+ }
+
+ ex integral::series(const relational & r, int order, unsigned options) const
+ {
+ if(x.subs(r)!=x)
+ throw std::logic_error("Cannot series expand wrt dummy variable");
+
+ // Expanding integrant with r substituted taken in boundaries.
+ ex fseries = f.series(r, order, options);
+ epvector fexpansion;
+ fexpansion.reserve(fseries.nops());
+ for(size_t i=0; i<fseries.nops(); ++i) {
+ ex currcoeff=ex_to<pseries>(fseries).coeffop(i);
+ fexpansion.push_back(expair(
+ currcoeff==Order(_ex1)
+ ? currcoeff
+ : integral(x, a.subs(r), b.subs(r), currcoeff),
+ ex_to<pseries>(fseries).exponop(i)
+ ));
+ }
+
+ // Expanding lower boundary
+ ex result=(new pseries(r, fexpansion))->setflag(status_flags::dynallocated);
+ ex aseries = (a-a.subs(r)).series(r, order, options);
+ fseries = f.series(x==(a.subs(r)), order, options);
+ for(size_t i=0; i<fseries.nops(); ++i) {
+ ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
+ if(is_order_function(currcoeff))
+ break;
+ ex currexpon = ex_to<pseries>(fseries).exponop(i);
+ int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
+ currcoeff=currcoeff.series(r, orderforf);
+ ex term = ex_to<pseries>(aseries)
+ .power_const(ex_to<numeric>(currexpon+1),order);
+ term=ex_to<pseries>(term).mul_const(ex_to<numeric>(-1/(currexpon+1)));
+ term=ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
+ result=ex_to<pseries>(result).add_series(ex_to<pseries>(term));
+ }
+
+ // Expanding upper boundary
+ ex bseries = (b-b.subs(r)).series(r, order, options);
+ fseries = f.series(x==(b.subs(r)), order, options);
+ for(size_t i=0; i<fseries.nops(); ++i) {
+ ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
+ if(is_order_function(currcoeff))
+ break;
+ ex currexpon = ex_to<pseries>(fseries).exponop(i);
+ int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
+ currcoeff=currcoeff.series(r, orderforf);
+ ex term = ex_to<pseries>(bseries)
+ .power_const(ex_to<numeric>(currexpon+1),order);
+ term=ex_to<pseries>(term).mul_const(ex_to<numeric>(1/(currexpon+1)));
+ term=ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
+ result=ex_to<pseries>(result).add_series(ex_to<pseries>(term));
+ }
+
+ return result;
}
Index: ginac/pseries.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/pseries.h,v
retrieving revision 1.36
diff -c -r1.36 pseries.h
*** ginac/pseries.h 4 Jan 2004 13:11:45 -0000 1.36
--- ginac/pseries.h 27 Oct 2004 07:57:04 -0000
***************
*** 56,61 ****
--- 56,62 ----
ex normal(exmap & repl, exmap & rev_lookup, int level = 0) const;
ex expand(unsigned options = 0) const;
ex conjugate() const;
+ ex eval_integ() const;
protected:
ex derivative(const symbol & s) const;
***************
*** 82,87 ****
--- 83,92 ----
/** Returns true if there is no order term, i.e. the series terminates and
* false otherwise. */
bool is_terminating() const;
+
+ /** Get coefficients and exponents. */
+ ex coeffop(size_t i) const;
+ ex exponop(size_t i) const;
ex add_series(const pseries &other) const;
ex mul_const(const numeric &other) const;
More information about the GiNaC-devel
mailing list