share/integration/delta.mac

   1
   2 delint(expr,var,[limits]):=block([arg,arglist,isd],
   3 arglist:cons(var,limits),
   4 if length(limits)#0 and length(limits)#2 then error(
   5 "wrong number of arguments to delint"),
   6 if not atom(expr) then
   7 if part(expr,0)="+" or (part(expr,0)="-" and length(expr)>1) then
   8         return((apply('delint,cons(part(expr,1),arglist))+apply
   9         ('delint,cons(expr-part(expr,1),arglist)))),
  10 isd:isdeltap(expr),
  11 if isd=false then return(apply('integrate,cons(expr,arglist))),
  12 arg:rhs(part(isd,2)),
  13 if freeof(var,arg) then return(apply('integrate,cons(expr,arglist))),
  14 apply('delintd,cons(expr,cons(arg,arglist))))$
  15
  16 /* loadfile(utils,fasl,mps)$ */
  17
  18 matchdeclare([dum1,dum2],true)$
  19
  20
  21 defmatch(isdeltap,dum1*delta(dum2))$
  22 delintd(grand,arg,var,[l]):=
  23 block([prederror:false,programmode:true,rest,nsol,solns,jac,ans],
  24 solns:map(rhs,solve(arg,var)),
  25 if (nsol:length(solns))=0 then return(0)
  26 ,solns:removeimag(solns),
  27 if length(l)=2 then solns:rem_out_of_range(solns,l[1],l[2]),
  28 if (nsol:length(solns))=0 then return(0),
  29 rest:grand/delta(arg),
  30 jac:abs(diff(arg,var)),
  31 if length(l)=2 then sum(subst((solns[n]),var,rest/jac),n,1,nsol)
  32         else sum(subst(solns[n],var,rest/jac)*theta(var-solns[n]),n,1,nsol))$