archive/share/trash/differ.mac

   1                         /* -*-Macsyma-*- */
   2 eval_when([batch],ttyoff:true)$
   3 define_variable(normalize,false,boolean)$
   4 define_variable(differverbose,false,boolean)$
   5 define_variable(index,'index,any)$
   6 define_variable(u,'u,any)$
   7 define_variable(array,'array,any)$
   8 /*  Linear difference equation package  */
   9
  10 magnitude(vector):=
  11 block([scalarmatrix:true],
  12       if listp(vector) or length(vector)=1
  13        then return (sqrt(vector.transpose(vector))),
  14       if length(transpose(vector))=1
  15        then return (sqrt(transpose(vector).vector)),
  16       print ("magnitude:  not a vector --",vector), false)$
  17
  18 /*  Returns a list of the eigenvalues of a matrix.
  19     Does not handle repeated eigenvalues.  */
  20
  21 eigenvalues(mx):=
  22 block([loadprint:false,programmode:true,charpoly,lambda,result],
  23       if length(mx) # length(transpose(mx))
  24       then (print ("eigenvalues:  not a square matrix --",mx), return(false)),
  25       charpoly:apply('charpoly,[mx,lambda]),
  26       result:solve(charpoly,lambda),
  27       makelist(rhs(x),x,result))$
  28
  29
  30 /*  Returns a list which is a transpose of the eigenvector  */
  31
  32 eigenvector(mx,eigenvalue):=
  33 block([loadprint:false,programmode:true,ttyoff:true,
  34        degree,xlist,xvector,eqnlist,result],
  35       degree:length(mx),
  36       if degree # length(transpose(mx))
  37       then (print("eigenvector:  not a square matrix --",mx), return(false)),
  38       xlist:makelist(concat('x,i),i,1,degree),
  39       xvector:transpose(matrix(xlist)),
  40       mx:mx.xvector-eigenvalue*xvector,
  41       eqnlist:makelist(mx[i,1]=0,i,1,degree),
  42       if member(0=0,eqnlist)
  43          then eqnlist:delete(0=0,eqnlist)
  44          else eqnlist:rest(eqnlist),
  45       result:first(solve(cons('x1=1,eqnlist),xlist)),
  46       result:makelist(part(x,2),x,result),
  47       if normalize=true then result/magnitude(result) else result)$
  48
  49 /* Puts equations in the form x[n+1] = a x[n] + b y[n] + c z[n] + ...        */
  50 /* Be careful for the case where eqn is x[n+1] = x[n+1] and var is x[n]      */
  51 /* SOLVE returns ALL for SOLVE(X=X, X), i.e. arbitrarily many solutions.     */
  52
  53 standardize(eqn,var):=
  54 block([x,y],
  55       x:part(var,0)[part(var,1)+1],
  56       y:solve(eqn,x),
  57       if y = [] then (print("difference:  no",x,"term --",eqn),
  58                       throw('missing_term)),
  59       if y = 'all then eqn else first(y)
  60       )$
  61 \f
  62 /*  Solves single first order equations of the form f[n+1] = a f[n] + b  */
  63
  64 first_order_difference(eqn,var):=
  65 block([a,b,simpsum:true],
  66       a:coeff(rhs(eqn),var),
  67       b:rhs(eqn)-a*var,
  68       var=a^index*array[0]+b*sum(a^k,k,0,index-1))$
  69
  70 /*  Solves a single higher order equation by converting to a first order system  */
  71
  72 second_order_difference(eqn,var):= block([],
  73       part(system([eqn,array[index+1]=array[index+1]],[array[index+1],var]),2,1))$
  74
  75 /*  Solves a system of first order equations  */
  76
  77 system(eqnlist,varlist):=
  78 block([listarith:true,normalize:false,a,u,lambdas,s,sinverse,d],
  79       eqnlist:map(lambda([x,y],rhs(standardize(x,y))),eqnlist,varlist),
  80       a:makelist(makelist(coeff(eqn,var),var,varlist),eqn,eqnlist),
  81       a:apply('matrix,a),
  82       u[0]:transpose(makelist(apply('ev,[var,solve(index=0)]),var,varlist)),
  83       u[index]:transpose(varlist),
  84       if differverbose=true then (apply('ldisplay,['u[index],'u[0],'a]),
  85                                   apply('ldisp,['u[index+1]='a*'u[index],
  86                                                 'u[index]=
  87                                                 '(s.lambda^n.s^^(-1).u[0])])),
  88       lambdas:eigenvalues(a),
  89       s:makelist(eigenvector(a,lambda),lambda,lambdas),
  90       s:transpose(apply('matrix,s)),
  91       sinverse:s^^-1,
  92       d:lambdas*ident(length(lambdas)),
  93       u[index]:s.d^index.sinverse.u[0],
  94       transpose(map("=",varlist,part(transpose(u[index]),1))))$
  95
  96 /*  General toplevel function  */
  97
  98 difference(eqn,var):=
  99 block([loadprint:false,programmode:true,array,index,higherorder],
 100       if listp(eqn) then (array:makelist(part(x,0),x,var),
 101                           index:part(first(var),1),
 102                           return (system(eqn,var))),
 103       array:part(var,0),
 104       index:part(var,1),
 105       higherorder:makelist(array[index+n],n,2,5),
 106       eqn:catch(standardize(eqn,var)),
 107       if eqn = 'missing_term then return(done),
 108 /*      if apply("AND",map(lambda([x],freeof(x,eqn)),higherorder))
 109 This code does not work in DOE MACSYMA due to operand error, so
 110 replaced by equivalent test. Leo Harten, Aug 22, 1984 */
 111       if not member('false,map(lambda([x],freeof(x,eqn)),higherorder))
 112          then first_order_difference(eqn,var)
 113          else second_order_difference(eqn,var))$
 114
 115 /* Bugs, comments to CWH@MC
 116    Things to do:
 117    Rewrite this in Lisp
 118    Put a catch around the call to standardize from difference.
 119    Become more intelligent about cases like f[n]=f[n-1]+f[n-2]. */
 120
 121 eval_when([batch],ttyoff:false)$