Si, Ci, Ei, also more pattern matching examples.

Am Donnerstag, den 11.12.2008, 17:10 -0800 schrieb Richard Fateman:
> Here's a set of rules for integrating some stuff with Si, Ci, and Ei.
> Is there another place where they are defined etc?
> 
> /* Integration of some forms including cos(x)/x^n and sin(x)/x^n in
> terms of Ci(x)=integrate(cos(x)/x,x) and Si(x)=integrate(sin(x)/x,x)*/
> 
> (matchdeclare(w,true,nn, nonnegintegerp),
> tellsimp('integrate(cos(w)/w^nn,w), cosint[nn](w)),
> tellsimp('integrate(sin(w)/w^nn,w), sinint[nn](w)),
> /*cosint[n](w) is integral of cos(w)/w^n */
> cosint[n](w):= block([k:1-n], (cos(w)*w^k +sinint[-k](w))/k),
> sinint[n](w):= block([k:1-n], (sin(w)*w^k -cosint[-k](w))/k),
> cosint[0](w):= sin(w),
> sinint[0](w):= -cos(w),
> cosint[1](w):= Ci(w),
> sinint[1](w):= Si(w),
> gradef(Ci(w), cos(w)/w),
> gradef(Si(w), sin(w)/w))
> 
> /* you may be surprised [I was] that this will find closed forms for
> 
> integrate (a*cos(x)*x^4+b*sin(x))/x^2,x);
> integrate(cos(sqrt(x))/x^2,x);
> integrate(cos(x+1)/(x+1)^2,x);
> integrate(cos(x^4)/x,x);
> integrate(x*cos(sqrt(x^2+2))/(x^2+2), x);
> integrate(cos(x)*cos(sin(x))/sin(x), x);
> integrate(-sin(x)*cos(cos(x))*cos(sin(cos(x)))/sin(cos(x)),x);
> 
> integrate(Ci(x)*cos(x)/x,x);
> integrate(Ci(x^2)*cos(x^2)/x,x);
> 

David Billinghurst has already mentioned that we have some new code
concerning the Exponential Integrals. Here are the integrals for the
Exponential Integral Ci shown above:

(%i42) integrate(expintegral_ci(x)*cos(x)/x,x);
(%o42) expintegral_ci(x)^2/2
(%i43) integrate(expintegral_ci(x^2)*cos(x^2)/x,x);
(%o43) expintegral_ci(x^2)^2/4

Only the elementary integral of the function is implemented. The results
with the additional cos function Maxima gets because of the implemented
derivative differentiation algorithm.

But the other integrals Maxima can not solve with this extensions up to
now. So it is very interesting for me to look at the implementation of
the other integrals too.

> 
> */
> 
> /* more rules, for x^n*Ci(x), x^n*Si(x), n is pos or neg but not -1. */ 
> (matchdeclare (notm1, lambda([r],is (r # -1))),
>  tellsimp('integrate(Ci(w)*w^notm1,w),
>           block([k:1+notm1],(Ci(w)*w^k -integrate(w^(k-1)*cos(w),w))/k)),
> (tellsimp('integrate(Si(w)*w^notm1,w),
>           block([k:1+notm1],(Si(w)*w^k -integrate(w^(k-1)*sin(w),w))/k))))
> 
> 
> 
> /*exponential integral, Ei(x)= integral(exp(x)/x,x)
>   matches don't work so easily, since the matcher has to look inside
>   both %e^x and x^n to see which is which, and back up. Here's a
>   way around.  Another hazard is that 'integrate(z,w) simplifies to z*w.
>   Is that a bug? To get around this we introduce mintegrate(). */
> 
> (matchdeclare([z,w,v],true,nnn, negintegerp),
> 
> negintegerp(n):=is (n<0) and integerp(n),
> 
> tellsimpafter(mintegrate(z,v), intwrtf(v,z)),
> 
> defmatch(exp_mat_,exp(w)*w^nnn,w),  /* if we know w, we can find exp(w), 
> w^n */
> 
> intwrtf(v,z):=block([ans:exp_mat_(z,v)],print(ans),
>      if (ans= false) then integrate(z,v) else
>                           ev(expint[-nnn](v),ans) ),
> /*expint[n](w) is integral of exp(w)*w^n */
> expint[n](w):= block([k:1-n], (exp(w)*w^k -expint[-k](w))/k),
> expint[1](w):=Ei(w),
> expint[0](w):=exp(w),
> gradef(Ei(w), exp(w)/w))
> 

Now some results for the Exponential Integral E:

(%i45) integrate(expintegral_e(n,x),x);
(%o45) -expintegral_e(n+1,x)

If we set an expansion flag the Exponential Integral E is represented
through related functions:

(%i46) expintexpand:true;
(%o46) true

(%i47) integrate(expintegral_e(n,x),x);
(%o47) -expintegral_e(n+1,x)
(%i48) integrate(expintegral_e(0,x),x);
(%o48) -gamma_incomplete(0,x)

The log functions in the following expansion take care of the most
general case including complex values and the negative real axis:

(%i49) integrate(expintegral_e(1,x),x);
(%o49) -x*log(x)-(-x*log(-x)+log(-1/x)*x+2*x)/2-x*expintegral_ei(-x)-%
e^-x+x

This is the expansion in terms of the power function for a negative
parameter:

(%i53) expintegral_e(-1,x);
(%o53) (x+1)*%e^-x/x^2

Because we have extended the integrator to handle more power functions
we get the integral in terms of the gamma incomplete function:

(%i50) integrate(expintegral_e(-1,x),x);
(%o50) -gamma_incomplete(0,x)-gamma_incomplete(-1,x)

Here are some further examples for the expansion:

(%i51) expintegral_e(0,x);
(%o51) %e^-x/x
(%i52) expintegral_e(1,x);
(%o52) -log(x)-(log(-1/x)-log(-x))/2-expintegral_ei(-x)

For half integral values we get an expansion in therms of the Error
function:

(%i54) expintegral_e(1/2,x);
(%o54) sqrt(%pi)*erfc(sqrt(x))/sqrt(x)
(%i55) expintegral_e(-1/2,x);
(%o55) %e^-x/x+sqrt(%pi)*erfc(sqrt(x))/(2*x^(3/2))

The biggest problem with all this code is that I have not written the
documentation. Perhaps I have more time in the next weeks to do more
work on the documentation.

Dieter Kaiser