Re: non-polynomial differential operator

From: Carl Devore (devore_at_math.udel.edu)
Date: 12/18/04 

Date: Fri, 17 Dec 2004 22:12:31 -0500

On Thu, 16 Dec 2004, Alain Verghote wrote:
> On Wed, 15 Dec 2004 13:48:35 -0500, Carl Devore wrote:
> >> psi:= proc(DiffOp, f, deg)
> >> local x,p,d,_D,k;
> >> try d:= degree(f(x),x) catch: d:= FAIL end try;
> >> if not d::nonnegint then d:= `if`(nargs=3, deg, Order) fi;
> >> p:= convert(taylor(subs([D= _D, I= 1], DiffOp), _D, d+2),
> 'polynom');
> >> unapply(simplify(add(coeff(p,_D,k)*D@@k, k= 0..d+1)(f)(x)), x)
> >> end:
>
> However we must care about negative powers of D or integrals
> Example:the sum of 8th powers 1^8+ 2^8+3^8+ ..x^8 has to be computed
> p(x)= (D/(I-exp(-D)))(x^9/9) instead of (I/(I-exp(-D)))(x^8)
> to avoid an error message:..not have a taylor expansion..
> giving the right Bernoulli polynom:
> x^9/9+x^8/2+2*x^7/3-7*x^5/15+2*x^3/9-x/30 !!!

Okay, it just needs a minor modification:

> INT:= F-> unapply(int(F, 0.._X), _X):
> unprotect(invfunc): invfunc[D]:= INT: protect(invfunc):
> psi:= proc(DiffOp, f, deg)
> local x,p,d,_D,k;
> try d:= degree(f(x),x) catch: d:= FAIL end try;
> if not d::nonnegint then d:= `if`(nargs=3, deg, Order) fi;
> p:= convert(series(subs([D= _D, I= 1], DiffOp), _D, d+2), 'polynom');
> unapply(simplify(add(coeff(p,_D,k)*D@@k, k= -d-1..d+1)(f)(x)), x)
> end:

> psi(I/(I-exp(-D)), x-> x^8);

                9 8 7 5 3
      x -> 1/9 x + 1/2 x + 2/3 x - 7/15 x + 2/9 x - 1/30 x

Relevant Pages