share/calculus/pade2.mac

   1
   2 /*
   3 This subroutine computes pade approximants of multivariable functions
   4 Copyright (C) 1999  Dan Stanger
   5
   6 This library is free software; you can redistribute it and/or modify it
   7 under the terms of the GNU Library General Public License as published
   8 by the Free Software Foundation; either version 2 of the License, or (at
   9 your option) any later version.
  10
  11 This library is distributed in the hope that it will be useful, but
  12 WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 Library General Public License for more details.
  15
  16 You should have received a copy of the GNU Library General Public
  17 License along with this library; if not, write to the Free Software
  18 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  19
  20 Dan Stanger dan.stanger@internut.com
  21 */
  22
  23 /* this macro was defined as length did not work correctly on arrays.
  24    put it in your startup file as autoload
  25 newArray1(a,i)::=buildq([a,i],(array(a,i),put(a,arrayinfo(a),"arrayinfo")));
  26 alength(x):=block([l:get(x,"arrayinfo")],if listp(l) then first(third(l)) else
  27    error(x," was not declared with newArray1"))$
  28 */
  29
  30 pade2(f,IndVarList, pointlist, norderlist,dorderlist):=
  31 block([atlist,TransIndVarList,transindvar,NumIndVar,norder,
  32         dorder,torder,makeIndexList,makeIndexListAux],
  33    NumIndVar:length(IndVarList),
  34    newArray1(transindvar,NumIndVar),
  35    newArray1(norder,length(norderlist)),
  36    newArray1(dorder,length(dorderlist)),
  37    newArray1(torder,max(length(norderlist),length(dorderlist))),
  38    atlist:maplist(lambda([l,r],l=r),IndVarList,pointlist),
  39    fillarray(norder,norderlist),
  40    fillarray(dorder,dorderlist),
  41    fillarray(torder,norderlist+dorderlist),
  42    TransIndVarList:maplist(lambda([l,r],l-r),IndVarList,pointlist),
  43    fillarray(transindvar,TransIndVarList),
  44
  45    makeIndexListAux:lambda([x,c],
  46       block([l:[]],for i:0 thru c do
  47          block([y:copylist(x)],push(i,y),push(y,l)),l)),
  48    makeIndexList:lambda([order],
  49       block([l],l:makeIndexListAux([],order[0]),
  50          for j:1 thru (alength(order)-1) do
  51             block([l1:[]],
  52                for i in l do
  53                   block(l1:append(l1,makeIndexListAux(i,order[j])),l:l1)),
  54                   sort(l))),
  55    block([nindexlist:makeIndexList(norder),
  56           dindexlist:rest(makeIndexList(dorder)),
  57           tindexlist:makeIndexList(torder),
  58           dcoeff, ncoeff, rcoeff, lcoeff, p],
  59       newArray1(dcoeff,length(dindexlist)),
  60       newArray1(ncoeff,length(nindexlist)),
  61       newArray1(rcoeff,length(tindexlist)),
  62       newArray1(lcoeff,length(tindexlist)),
  63       p:block([den:1,num:0,i:1],
  64          for d in dindexlist do
  65             block(den:den+dcoeff[i]*
  66                 apply("*",maplist("^",TransIndVarList,d)),i:i+1),
  67          i:0,
  68          for n in nindexlist do
  69             block(num:num+ncoeff[i]*
  70                 apply("*",maplist("^",TransIndVarList,n)),i:i+1),
  71         num/den),
  72       block([i:0,lf:[f],lp:[p]],
  73          for r in tindexlist do
  74             block([l:[]],
  75                 maplist(lambda([x,y],l:append([x,y],l)),IndVarList,r),
  76                 rcoeff[i]:at(apply('diff,append(lf,l)),atlist),
  77                 lcoeff[i]:at(apply('diff,append(lp,l)),atlist),
  78                 i:i+1)),
  79       block([l:[],v:[],s],
  80          for i:0 thru alength(rcoeff)-1 do
  81                 l:append([lcoeff[i]=rcoeff[i]],l),
  82          for i:0 thru alength(ncoeff)-1 do
  83                 v:append(v,[ncoeff[i]]),
  84          for i:1 thru alength(dcoeff) do
  85                 v:append(v,[dcoeff[i]]),
  86          s:solve(l,v),
  87          at(p,first(s))))
  88 );
  89
  90 /* here is a example and the result
  91 pade2(exp(x+y),[x,y],[0,0],[1,1],[1,1])
  92
  93 ((((x * y)/4) + (y/2) + (x/2) + 1)/(((x * y)/4) - (y/2) - (x/2) + 1))
  94
  95 */