[GiNaC-list] Indefinite integration support

[Jorge Cardona]

Hi,

I recently wrote a code that will let define a entity that represents
a indefinite integral, i called iintegral, basically i took the
integral class, and remove the integration range, then i add to the
function.pl extra code to define a new option on function, called
iintegral_func, with that is possible to define the indefinite
integral of a function when we define the function itself.

I want to know if you like that idea and if you have some extra
comment on this part.

Look a definition of a function on inifcns_trans :

static ex log_iinteg(const ex & x, unsigned iinteg_param)
{
    GINAC_ASSERT(iinteg_param==0);

    // int {log(x)} -> x * (log (x) - 1)
    return x * (log (x) - _ex1);
}

....

REGISTER_FUNCTION(log, eval_func(log_eval).
                       evalf_func(log_evalf).
                       derivative_func(log_deriv).
                       iintegral_func(log_iinteg).
                       series_func(log_series).
                       real_part_func(log_real_part).
                       imag_part_func(log_imag_part).
                       latex_name("\\ln"));

...

Right now i call the iintegral constructor with 2 arguments, the
integrand expression and the integration variable (the integration
variable remains after the indefinite integral, is not as dummy as in
the definite integral), and this iintegral instance remains until the
eval_integ is called, which in a future (near future i hope) will have
a better algorithm to integrate the expression, but right now just
have the code that integrate a polynomial expression and a code that
check if the expression is a functions, and if that function has an
integral defined it also check if the argument is equal to the
integration variable, and then call that function, if any of this
don't match it return the iintegral object.

I want to know what do you think about that way, or if you have some
extra ideas, basically i want to implement the integration algorithm,
but don't want to have some code that wont be added to GiNaC, so,
please let me know your comments.


Code on function::iintegral :

ex function::iintegral(const symbol & s) const
{
    GINAC_ASSERT(serial<registered_functions().size());
    const function_options &opt = registered_functions()[serial];

    // No iintegral defined? Then return abstract iintegral object
    if (opt.iintegral_f == NULL)
        return GiNaC::iintegral(s,this->eval());

    size_t iinteg_param = -1;
    size_t num = seq.size();
    for (size_t i=0; i<num; i++) {
        ex a = s;
        ex b = seq[i];
        if(a.is_equal(b))
        {
            iinteg_param = i;
            break;
        }
    }

    current_serial = serial;
    if (opt.iintegral_use_exvector_args)
        return ((iintegral_funcp_exvector)(opt.iintegral_f))(seq, iinteg_param);


    if(iinteg_param < 0)
        throw(std::logic_error("Composed functions on integral."));

    switch (opt.nparams) {
        // the following lines have been generated for max. 14 parameters
    case 1:
        return ((iintegral_funcp_1)(opt.iintegral_f))(seq[1-1],iinteg_param);
    case 2:
        return ((iintegral_funcp_2)(opt.iintegral_f))(seq[1-1],
seq[2-1],iinteg_param);
    case 3:
        return ((iintegral_funcp_3)(opt.iintegral_f))(seq[1-1],
seq[2-1], seq[3-1],iinteg_param);
    case 4:
        return ((iintegral_funcp_4)(opt.iintegral_f))(seq[1-1],
seq[2-1], seq[3-1], seq[4-1],iinteg_param);
...(Go with more arguments)


Code on iintegral::eval_integ():

ex iintegral::eval_integ() const
{
    if (!(flags & status_flags::expanded))
        return this->expand().eval_integ();

    if (f==x)
        return x*x/2;
    if (is_a<power>(f) && f.op(0)==x) {
        if (f.op(1)==-1)
            return log(x);
        if (!f.op(1).has(x)) {
            ex primit = power(x,f.op(1)+1)/(f.op(1)+1);
            return primit;
        }
    }

    if (is_a<function>(f)) {
        function fun = ex_to<function>(f);
        // TODO: Change the signature.
        return fun.iintegral(ex_to<symbol>(x));
    }

    return *this;
}


There is a patch at:
http://wiki.scilab.org/Gsoc2009Symbolic?action=AttachFile&do=get&target=ExtendingIntegrals.patch

This is a little code test:

#include "ginac/ginac.h"

using namespace std;
using namespace GiNaC;

int main()
{
    symbol x("x");
    ex f,g,h;

    f = 2*x;
    cout << f << endl;

    g = iintegral(x,f);
    cout << g << endl;

    h = g.eval_integ();
    cout << h << endl;

    f = 2*sin(x) + exp(x) + log(x);
    cout << f << endl;

    g = iintegral(x,f);
    cout << g << endl;

    h = g.eval_integ();
    cout << h << endl;
}


That returns:

jcardona at terminus:~/stuff/personal/symbolic/spikes/Integrals$ ./test
2*x
iintegral(x,2*x)
x^2
exp(x)+log(x)+2*sin(x)
iintegral(x,exp(x)+log(x)+2*sin(x))
exp(x)-2*cos(x)+(-1+log(x))*x

Bye.


-- 
Jorge Eduardo Cardona
jorgeecardona at gmail.com
jorgeecardona.blogspot.com
------------------------------------------------
Linux registered user  #391186
Registered machine    #291871
------------------------------------------------