share/calculus/fourie.mac

   1 define_variable(sinnpiflag,true,boolean)$
   2 define_variable(cosnpiflag,true,boolean)$
   3
   4 remfun1(fun,exp):=scanmap(lambda([q],delfun1(fun,q)),exp)$
   5
   6 delfun1(fun,exp):=if not atom(exp) and inpart(exp,0) = fun
   7          then first(args(exp)) else exp$
   8
   9 remfunn1(fun,exp):=scanmap(lambda([q],delfunn1(fun,q)),exp)$
  10
  11 delfunn1(fun,exp):=if not atom(exp) and inpart(exp,0) = fun
  12           then -first(args(exp)) else exp$
  13
  14 remfun2(fun,exp,var):=scanmap(lambda([q],delfun2(fun,q,var)),exp)$
  15
  16 delfun2(fun,exp,var):=
  17      if not atom(exp) and inpart(exp,0) = fun and member(var,listofvars(exp))
  18          then first(args(exp)) else exp$
  19
  20 remfunn2(fun,exp,var):=scanmap(lambda([q],delfunn2(fun,q,var)),exp)$
  21
  22 delfunn2(fun,exp,var):=
  23       if not atom(exp) and inpart(exp,0) = fun
  24                        and member(var,listofvars(exp))
  25           then -first(args(exp)) else exp$
  26
  27 remfun(fun,exp,[var]):=if var = [] then remfun1(fun,exp)
  28         else (if length(var) = 1 then remfun2(fun,exp,first(var))
  29                   else error("too many arguments to remfun"))$
  30
  31 remfunn(fun,exp,[var]):=if var = [] then remfunn1(fun,exp)
  32          else (if length(var) = 1 then remfunn2(fun,exp,first(var))
  33                    else error("too many arguments to remfunn"))$
  34
  35 funp1(fun,exp):=block([inflag],inflag:true,
  36       if mapatom(exp) then false
  37           else (if inpart(exp,0) = fun then true
  38                     else member(true,maplist(lambda([q],funp1(fun,q)),exp))))$
  39  /* commented out of DOE MACSYMA
  40 ; APPLY("OR",MAPLIST(LAMBDA([Q],FUNP1(FUN,Q)),EXP))))$ */
  41
  42
  43 funp2(fun,exp,var):=block([inflag],inflag:true,
  44       if mapatom(exp) then false
  45           else (if inpart(exp,0) = fun and member(var,listofvars(exp))
  46                     then true
  47                     else member(true,
  48                                 maplist(lambda([q],funp2(fun,q,var)),exp))))$
  49 /* commented out of DOE MACSYMA
  50 APPLY("OR",MAPLIST(LAMBDA([Q],FUNP2(FUN,Q,VAR)),EXP))))$ */
  51
  52 funp(fun,exp,[var]):=if var = [] then funp1(fun,exp)
  53       else (if length(var) = 1 then funp2(fun,exp,first(var))
  54                 else error("too many arguments to funp"))$
  55
  56 equalp(x,y):=block([prederror],prederror:false,
  57        if is(equal(x,y)) = true then true)$
  58
  59 parity(f%,x):=if evenfunp(f%,x) then 'even
  60         else (if oddfunp(f%,x) then 'odd else 'neither)$
  61
  62 evenfunp(f%,x):=if equalp(f%,subst(-x,x,f%)) then true$
  63
  64 oddfunp(f%,x):=if equalp(f%,-subst(-x,x,f%)) then true$
  65
  66 parint(f%,x,a,b):=
  67     if not (equalp(-a,b) or a = 'minf and b = 'inf or a = 'inf and b = 'minf)
  68         then f%
  69         else (if atom(f%) or inpart(f%,0) # "+" then parint1(f%,x)
  70                   else map(lambda([q],parint1(q,x)),f%))$
  71
  72 parint1(f%,x):=if oddfunp(f%,x) then 0 else f%$
  73
  74 adefint(f%,x,a,b):=block([asign,bsign],
  75   if equalp(a,b) then 0
  76     else (if not (freeof(%i,f%) and freeof(%i,a) and freeof(%i,b))
  77             then ldefint(f%,x,a,b)
  78             else (f%:parint(f%,x,a,b),
  79                   if equalp(f%,0) then 0
  80                     else (if not funp2('abs,f%,x)
  81                             then ldefint(f%,x,a,b)
  82                             else (asign:asksign(a), bsign:asksign(b),
  83                                   if (asign = 'neg or asign = 'zero) and (bsign = 'neg or bsign = 'zero)
  84                                     then ldefint( remfunn2('abs,f%,x),x,a,b)
  85                                     else (if (asign = 'zero or asign = 'pos) and (bsign = 'zero or bsign = 'pos)
  86                                             then ldefint( remfun2('abs,f%,x),x,a,b)
  87                                             else (if asign = 'neg
  88                                                     then ratsimp( ldefint( remfunn2( 'abs,f%,x),x,a, 0)
  89                                                                  +ldefint( remfun2( 'abs,f%,x),x, 0,b))
  90                                                     else ratsimp( ldefint( remfun2( 'abs,f%,x),x,a, 0)
  91                                                                  +ldefint( remfunn2( 'abs,f%,x),x, 0,b)))))))))$
  92
  93 indefint(f%,x,halfplane):=if halfplane = 'pos
  94           then integrate(remfun2('abs,f%,x),x)
  95           else (if halfplane = 'neg then integrate(remfunn2('abs,f%,x),x)
  96                     else (if halfplane = 'both
  97                               then append(
  98                               ldisp(integrate(remfunn2('abs,f%,x),x)),
  99                               ldisp(integrate(remfun2('abs,f%,x),x)))
 100                               else error("invalid halfplane:",halfplane)))$
 101
 102 absint(f%,x,[range]):=if range = [] then indefint(f%,x,'pos)
 103         else (if length(range) = 1 then indefint(f%,x,first(range))
 104                   else (if length(range) = 2
 105                             then adefint(f%,x,range[1],range[2])
 106                             else error("too many arguments to absint")))$
 107
 108 fourint(f%,x):=block([z],
 109 if evenfunp(f%,x) then append(fourintcos(f%,x),ldisp(b[z] = 0))
 110          else (if oddfunp(f%,x) then append(ldisp(a[z] = 0),fourintsin(f%,x))
 111                    else fourintcoeff(f%,x)))$
 112
 113 fourintcoeff(f%,x):=block([az,bz,z],assume(z > 0),
 114              az:adefint(f%*cos(z*x),x,'minf,'inf)/%pi,
 115              bz:adefint(f%*sin(z*x),x,'minf,'inf)/%pi,
 116              append(ldisp(a[z] = az),ldisp(b[z] = bz)))$
 117
 118 fourintcos(f%,x):=block([az,z],assume(z > 0),
 119            az:2*adefint(f%*cos(z*x),x,0,'inf)/%pi,ldisp(a[z] = az))$
 120
 121 fourintsin(f%,x):=block([bz,z],assume(z > 0),
 122            bz:2*adefint(f%*sin(z*x),x,0,'inf)/%pi,ldisp(b[z] = bz))$
 123
 124 fourier(f%,x,p):=block([n],
 125 if evenfunp(f%,x) then append(fourcos(f%,x,p),ldisp(b[n] = 0))
 126          else (if oddfunp(f%,x)
 127                    then append(ldisp(a[0] = 0),ldisp(a[n] = 0),
 128                                foursin(f%,x,p)) else fourcoeff(f%,x,p)))$
 129
 130 fourcoeff(f%,x,p):=block([a0,an,bn,n],assume(n > 0),
 131           a0:adefint(f%,x,-p,p)/(2*p),an:adefint(f%*cos(n*%pi*x/p),x,-p,p)/p,
 132           bn:adefint(f%*sin(n*%pi*x/p),x,-p,p)/p,
 133           append(ldisp(a[0] = a0),ldisp(a[n] = an),ldisp(b[n] = bn)))$
 134
 135 fourcos(f%,x,p):=block([a0,an,n],assume(n > 0),a0:adefint(f%,x,0,p)/p,
 136         an:2*adefint(f%*cos(n*%pi*x/p),x,0,p)/p,
 137         append(ldisp(a[0] = a0),ldisp(a[n] = an)))$
 138
 139 foursin(f%,x,p):=block([bn,n],assume(n > 0),
 140         bn:2*adefint(f%*sin(n*%pi*x/p),x,0,p)/p,ldisp(b[n] = bn))$
 141
 142 foursimp(exp):=if listp(exp)
 143           then map(lambda([q],first(ldisp(foursimp(apply('ev,[q]))))),exp)
 144           else (if not freeof("=",exp) then lhs(exp) = foursimple(rhs(exp))
 145                     else foursimple(exp))$
 146
 147 foursimple(exp):=block([n],
 148            if funp1('integrate,exp) then exp
 149                else (if sinnpiflag then exp:ratsubst(0,sin(n*%pi),exp),
 150                      if cosnpiflag then exp:ratsubst((-1)^n,cos(n*%pi),exp),
 151                      factor(exp)))$
 152
 153 fourexpand(l,x,p,nn):=block([simpsum,series,l1,lhsl1,n],
 154            if not listp(l) then error("first argument not a list")
 155                else (if l = [] then error("argument list is empty")
 156                          else (l:apply('ev,[l]),simpsum:true,series:0,
 157                                unless l = [] do
 158                                       (l1:first(l),l:rest(l),lhsl1:lhs(l1),
 159                                        if lhsl1 = a[0]
 160                                            then series:series+rhs(l1)
 161                                            else (if lhsl1 = a[n]
 162                                                      then series
 163                                                      :series
 164                                                       +apply('sum,
 165                                                              [rhs(l1)
 166                                                                *cos(
 167                                                                 n*%pi*x/p),n,
 168                                                               1,nn])
 169                                                      else (if lhsl1 = b[n]
 170                                                                then series
 171                                                                :series
 172                                                                 +apply(
 173                                                                  'sum,
 174                                                                  [
 175                                                                   rhs(l1)
 176                                                                    *sin(
 177                                                                     n*%pi*x
 178                                                                      /p),n,1,
 179                                                                   nn])
 180                                                                else error(
 181                                                                "invalid equati
 182 on in argument list:",
 183                                                                l1)))),
 184                                series)))$
 185
 186 totalfourier(f%,x,p):=fourexpand(foursimp(fourier(f%,x,p)),x,p,'inf)$