Artifact 9f66b62769d68b1bc105156466302c10f745abe554e89580bf00b6f9c3011775:
- File
r34.1/lib/fide1.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 86748) [annotate] [blame] [check-ins using] [more...]
COMMENT ***** NOTE *****; % For fastloading of this package please follow these steps: % faslout fide1; % in "fide1.red"; % faslend; % load "fide1","gentran"; % faslout fide; % in "fide.red"$ % faslend; % For running this package the matrix, gentran, fide1 and fide packages % have to be loaded. However, loading fide is enough to make this % happen, since matrix, gentran and fide1 are automatically loaded. % Version 1.1 of the FIDE package is the result of porting the FIDE % package version 1.0.0 to REDUCE 3.4. Some addition presented in the % LINBAND module. off echo$ load!-package 'matrix; %*********************************************************************** %***** ***** %***** P a c k a g e F I D E Ver. 1.1.0 17/03/1992 ***** %***** ***** %*********************************************************************** %** ** %**FInite difference method for partial Differential Equation systems ** %** ** %*********************************************************************** %** (C) 1991, Richard Liska ** %** Faculty of Nuclear Science and Physical Engineering ** %** Technical University of Prague ** %** Brehova 7 ** %** 115 19 Prague 1 ** %** Czechoslovakia ** %** Email: Richard Liska <tjerl@cspuni12.bitnet> ** %** ** %** This package can be distributed through REDUCE Network Library. ** %*********************************************************************** % The FIDE package consists of the following modules: % % EXPRES for transforming PDES into any orthogonal coordinate system. % IIMET for discretization of PDES by integro-interpolation method. % APPROX for determining the order of approximation of difference % scheme % CHARPOL for calculation of amplification matrix and characteristic % polynomial of difference scheme, which are needed in Fourier % stability analysis. % HURWP for polynomial roots locating necessary in verifying the von % Neumann stability condition. % LINBAND for generating the block of FORTRAN code, which solves a % system of linear algebraic equations with band matrix % appearing quite often in difference schemes. % % These modules can also be used separately. % Changes since version 1.1: % Patches in SIMPINTERPOL and SIMPINTT 13/06/91 % Patch in TDIFPAIR 08/07/91 % Two FEXPR routines F2VAL, FPLUS changed to MACROs 17/03/92 %*********************************************************************** %***** ***** %***** M O D U L E E X P R E S ***** %***** ***** %*********************************************************************** module expres$ % Author: R. Liska % Program EXPRES, Version REDUCE 3.4 05/1991 symbolic$ global '(!*outp)$ % declarations for 3.4 fluid '(!*wrchri orig!* posn!*)$ %global '(!*outp orig!* posn!*)$ % declarations for 3.3 %fluid '(!*wrchri)$ switch wrchri$ global '(olddimension!* dimension!* coordindx!* cyclic1!* cyclic2!* sfprod!* nscal!*)$ flag('(share),'eval)$ % So that SHARE recognized by FASL. share olddimension!*,dimension!*,coordindx!*,cyclic1!*,cyclic2!*, sfprod!*,nscal!*$ nscal!*:=0$ put('tensor,'tag,'tens)$ put('tensor,'fn,'tensfn)$ put('tensor,'evfn,'expres)$ put('tens,'prifn,'tenspri)$ flag('(tensor),'sprifn)$ put('tens,'setprifn,'settenspri)$ put('tensor,'typeletfn,'tenslet)$ symbolic procedure ptensor x$ 'tensor$ symbolic procedure poptensor u$ if null u then 'tensor else nil$ deflist('((tensor ptensor) (tensop poptensor) (df getrtypecar) (diff getrtypecar) (!& ptensor) (!# ptensor) (!? ptensor) (grad ptensor) (div ptensor) (lapl ptensor) (curl ptensor) (vect ptensor) (dyad ptensor) (dirdf ptensor)),'rtypefn)$ put('cons,'rtypefn,'getrtypecons)$ put('rcons,'evfn,'evrcons)$ remprop('cons,'psopfn)$ symbolic procedure getrtypecons u$ if getrtypecar u eq 'tensor then 'tensor else 'rcons$ symbolic procedure evrcons(u,v)$ rcons cdr u$ symbolic procedure tensor u$ for each a in u do <<put(a,'rtype,'tensop)$ put(a,'avalue,list('tensor,mktensor(0 , (0 . 1)))) >>$ deflist('((tensor rlis)),'stat)$ symbolic procedure tenslet(u,v,typu,b,typv)$ if not atom u then lprie list(" Non atom tensor variable ",u) else if b then <<if not eqcar(v,'tensor) then v:= mktensor(0, if eqcar(v,'!*sq) then cadr v else simp!* v)$ rmsubs()$ put(u,'rtype,'tensop)$ put(u,'avalue,list('tensor,v)) >> else <<remprop(u,'rtype)$ remprop(u,'avalue) >>$ %====================================================================== % Data structure for tensor quantities % ==================================== % (tensor nr rnk val car !*sqvar!*) % nr - integer, should be equal to actual nscal!*, otherwise % the quantity has been defined in previous coor. system % number of coordinate system % rnk - integer, 0,1,2 % rank of the tensor % 0 - scalar % 1 - vertor % 2 - dyad (matrix) % val - value % s.q. for rnk = 0 % list of s.q.s for rnk = 1 % list of lists of s.q.s for rnk = 2 % !*sqvar!* used in resimplification routine %====================================================================== % Smacro definitions for access of data structure subparts %====================================================================== smacro procedure tensrnk u$ % determines rank from cddr of datastructure car u$ smacro procedure tensval u$ % determines value from cddr of datastructure cadr u$ symbolic procedure mktensor(rnk,u)$ 'tensor . nscal!* . rnk . u . if !*resubs then !*sqvar!* else list nil$ symbolic procedure settenspri(u,v)$ if not atom u then lprie list(" Non-atom tensor variable ",u) else <<prin2!* u$ prin2!* get('setq,'prtch)$ tenspri v >>$ symbolic procedure tenspri u$ begin scalar rnk$ u:=cddr u$ rnk:=car u$ u:=cadr u$ if rnk=0 then <<pmaprin u$ terpri!* t >> else if rnk=1 then <<tenspri1 u$ orig!*:=0$ terpri!* t >> else if rnk=2 then <<prin2!* " ( "$ tenspri1 car u$ for each i in cdr u do <<prin2!* " , "$ terpri!* t$ tenspri1 i >>$ prin2!* " ) "$ orig!*:=0$ terpri!* t >> else lprie list(" Can't print tensor ",u," with rank ",rnk) end$ symbolic procedure tenspri1 u$ <<prin2!* " ( "$ orig!*:=posn!*$ pmaprin car u$ for each i in cdr u do <<prin2!* " , "$ terpri!* t$ pmaprin i >>$ orig!*:=orig!* - 3$ prin2!* " ) " >>$ symbolic procedure pmaprin u$ maprin(!*outp:=prepsq!* u)$ symbolic procedure updatedimen()$ if olddimension!* = dimension!* then t else <<if dimension!* =2 then <<coordindx!* :='(2 1)$ cyclic1!* :='(1 2)$ cyclic2!* :='(2 1) >> else if dimension!* =3 then <<coordindx!* :='(3 2 1)$ cyclic1!* :='(2 3 1)$ cyclic2!* :='(3 1 2) >> else lprie list(" Can't handle dimension = ",dimension!*)$ olddimension!* := dimension!* >>$ symbolic procedure expres(expn,v)$ express expn$ symbolic procedure resimptens u$ mktensor(caddr u,if caddr u=0 then resimp cadddr u else if caddr u=1 then for each a in cadddr u collect resimp a else if caddr u=2 then for each a in cadddr u collect for each b in a collect resimp b else lprie list("Can't handle tensor ",u, " with rank ",caddr u))$ symbolic procedure express expn$ begin scalar lst,matrx,rnk,opexpress$ if not atom expn then go to op$ if get(expn,'rtype) eq 'tensop and (rnk:=get(expn,'avalue)) and car rnk memq '(tensor tensop) and (rnk:=cadr rnk) then return if cadr rnk=nscal!* then if car cddddr rnk then rnk else resimptens rnk else lprie list(" You must rebind tensor ",expn, " in the new coordinate system")$ return mktensor(0,simp!* expn)$ op:if car expn = 'vect then return mktensor (1,testdim1 mapcar(cdr expn,function simp!*)) else if car expn = 'dyad then return mktensor (2,testdim2 mapcar(cdr expn,function (lambda a$ mapcar(a,function simp!*)))) else if car expn eq 'tensor then return if cadr expn=nscal!* then if car cddddr expn then expn else resimptens expn else lprie list(" You must rebind tensor ",expn, " in the new coordinate system")$ lst:=mapcar(cdr expn,function (lambda a$ cddr express a))$ if (opexpress:=get(car expn,'express)) then <<rnk:=eval(opexpress . list mkquote lst)$ return mktensor(car rnk,cadr rnk)>>$ if get(car expn,'simpfn) then return mktensor(0,simp( car expn . mapcar(lst,function(lambda a$ if car a=0 then list('!*sq,cdr a,nil) else typerr(expn," well formed scalar "))) ))$ lst:=mapcar(lst,function(lambda a$ if car a=0 then prepsq cdr a else typerr(expn," well formed tensor")))$ return mktensor(0,!*k2q(car expn.lst)) end$ procedure testdim1 u$ if length u=dimension!* then u else <<lprie "Bad number of vector components"$u>>$ procedure testdim2 u$ begin scalar x$ if length u = dimension!* then t else go to er$ x:=u$ a:if length car u = dimension!* then t else go to er$ x:=cdr x$ if x then go to a$ return u$ er:lprie "Bad number of dyad components"$ return u end$ %====================================================================== % Procedures in EXPRESS properties of operators are returning % (rnk val), their argument is list of (rnk val) symbolic procedure vectors arglist$ for each i in arglist do <<put(i,'rtype,'tensop)$ put(i,'simpfn,'simpiden)$ put(i,'avalue,list('tensop, mktensor(1,reverse for each a in coordindx!* collect simp list(i,a) )) )>>$ deflist('((vectors rlis)),'stat)$ symbolic procedure dyads arglist$ for each i in arglist do <<put(i,'rtype,'tensop)$ put(i,'simpfn,'simpiden)$ put(i,'avalue,list('tensop, mktensor(2,reverse for each a in coordindx!* collect reverse for each b in coordindx!* collect simp list(i,a,b)))) >>$ deflist('((dyads rlis)),'stat)$ symbolic procedure plusexpress u$ begin scalar z$ z:=car u$ a:u:=cdr u$ if null u then return z$ z:=plus2ex(z,car u)$ go to a end$ put('plus,'express,'plusexpress)$ symbolic procedure plus2ex(x,y)$ begin scalar mtx,mty,slx,sly,rnk,ans,ans1$ rnk:=tensrnk x$ if not(rnk=tensrnk y) then lprie "Tensor mishmash"$ if rnk=0 then return list(rnk,addsq(cadr x,cadr y)) else if rnk=1 then <<slx:=tensval x$ sly:=tensval y$ while slx do <<ans:=addsq(car slx,car sly) . ans$ slx:=cdr slx$ sly:=cdr sly>>$ ans:= list(1,reverse ans) >> else if rnk=2 then <<mtx:=tensval x$ mty:=tensval y$ while mtx do <<slx:=car mtx$ sly:=car mty$ ans1:=nil$ while slx do <<ans1:=addsq(car slx,car sly) . ans1$ slx:=cdr slx$ sly:=cdr sly>>$ ans:=reverse ans1 . ans$ mtx:=cdr mtx$ mty:=cdr mty>>$ ans:=list(2,reverse ans) >>$ return ans end$ symbolic procedure timesexpress u$ begin scalar z$ z:=car u$ a:u:=cdr u$ if null u then return z$ z:=times2ex(z,car u)$ go to a end$ put('times,'express,'timesexpress)$ symbolic procedure times2ex(x,y)$ begin scalar rnkx,rnky$ rnkx:=tensrnk x$ rnky:=tensrnk y$ return if rnkx=0 then list(rnky,times0ex(tensval x,tensval y,rnky)) else if rnky=0 then list(rnkx,times0ex(tensval y,tensval x,rnkx)) else lprie " Tensor mishmash " end$ symbolic procedure times0ex(x,y,rnk)$ if rnk=0 then multsq(x,y) else if rnk=1 then for each a in y collect multsq(x,a) else if rnk=2 then for each a in y collect for each b in a collect multsq(x,b) else lprie " Tensor mishmash "$ symbolic procedure minusexpress expn$ timesexpress list(list(0,cons(-1,1)),car expn)$ put('minus,'express,'minusexpress)$ symbolic procedure differenceexpress expn$ plusexpress list(car expn,minusexpress list cadr expn)$ put('difference,'express,'differenceexpress)$ symbolic procedure quotientexpress expn$ if tensrnk cadr expn = 0 then times2ex(list(0,simp!* list('recip,prepsq tensval cadr expn)), car expn) else lprie " Tensor mishmash "$ put('quotient,'express,'quotientexpress)$ symbolic procedure exptexpress expn$ if tensrnk car expn=0 and tensrnk cadr expn = 0 then list(0,simp!* list('expt, prepsq tensval car expn, prepsq tensval cadr expn)) else lprie " Tensor mishmash "$ put('expt,'express,'exptexpress)$ symbolic procedure recipexpress expn$ if tensrnk car expn = 0 then list(0,simp!* list('recip, prepsq tensval car expn)) else lprie " Tensor mishmash "$ put('recip,'express,'recipexpress)$ symbolic procedure inprodexpress expn$ begin scalar arg1,arg2,rnk1,rnk2$ arg1:=tensval car expn$ arg2:=tensval cadr expn$ rnk1:=tensrnk car expn$ rnk2:=tensrnk cadr expn$ return if rnk1=1 then inprod1ex(arg1,arg2,rnk2) else if rnk1=2 then inprod2ex(arg1,arg2,rnk2) else lprie " Tensor mishmash " end$ put('cons,'express,'inprodexpress)$ symbolic procedure inprod1ex(x,y,rnk)$ begin scalar lstx,lsty,mty,z,zz$ lstx:=x$ lsty:=y$ if rnk=1 then <<z:=nil . 1$ while lstx do <<z:=addsq(z,multsq(car lstx,car lsty))$ lstx:=cdr lstx$ lsty:=cdr lsty>>$ z:=list(0,z)>> else if rnk=2 then <<z:=nil$ lsty:=mty:=copy2(y,t)$ while car mty do <<lstx:=x$ lsty:=mty$ zz:=nil . 1$ while lstx do <<zz:=addsq(zz,multsq(car lstx,caar lsty))$ rplaca(lsty,cdar lsty)$ lsty:=cdr lsty$ lstx:=cdr lstx>>$ z:=nconc(z,list zz) >>$ z:=list(1,z)>> else lprie " Tensor mishmash "$ return z end$ symbolic procedure inprod2ex(x,y,rnk)$ begin scalar lstx,lsty,mtx,z,zz$ mtx:=x$ if rnk=1 then while mtx do <<z:=nconc(z,cdr inprod1ex(car mtx,y,1))$ mtx:=cdr mtx>> else if rnk=2 then while mtx do <<z:=nconc(z,cdr inprod1ex(car mtx,copy2(y,t),2))$ mtx:=cdr mtx>> else lprie " Tensor mishmash "$ return list(rnk,z) end$ symbolic procedure outexpress expn$ begin scalar x,y,z$ x:=tensval car expn$ y:=tensval cadr expn$ if tensrnk car expn=1 and tensrnk cadr expn=1 and null cddr expn then for each i in x do z:=mapcar(y,function(lambda a$ multsq(a,i))) . z else lprie list(" Outer product of ",expn)$ return list(2,reverse z) end$ put('!&,'express,'outexpress)$ flag('(!&),'tensfn)$ symbolic procedure copy2(x,p)$ if null x then nil else if p then copy2(car x,nil) . copy2(cdr x,t) else car x . copy2(cdr x,nil)$ symbolic procedure listar(arg,j)$ if j=1 then car arg else if j=2 then cadr arg else if j=3 then caddr arg else lprie list(" LISTAR ",arg,j)$ symbolic procedure listarsq(arg,j)$ prepsq listar(arg,j)$ symbolic procedure dinprod expn$ begin scalar x,y,z,xx,yy$ x:=tensval car expn$ y:=copy2(tensval cadr expn,t)$ z:=nil . 1$ if not(tensrnk car expn=2 and tensrnk cadr expn=2 and null cddr expn) then lprie list(" D-scalar product of ",expn)$ a:if null x and null y then go to d else if null x or null y then go to er$ xx:=car x$ yy:=car y$ b:if null xx and null yy then go to c else if null xx or null yy then go to er$ z:=addsq(z,multsq(car xx,car yy))$ xx:=cdr xx$ yy:=cdr yy$ go to b$ c:x:=cdr x$ y:=cdr y$ go to a$ d:return list(0,z)$ er:lprie list(" EXPRESS error ",expn," D-S dyads of dif. size") end$ put('!#,'express,'dinprod)$ put('hash,'express,'dinprod)$ put('hash,'rtypefn,'ptensor)$ symbolic procedure antisymsum(u,v)$ if dimension!* = 2 then difmul(car u,cadr u,cadr v,car v) else if dimension!* = 3 then list (difmul(cadr u,caddr u,caddr v,cadr v), difmul(caddr u,car u,car v,caddr v), difmul(car u,cadr u,cadr v,car v)) else lprie list(" ANTISYMSUM ",u,v)$ symbolic procedure difmul(a,b,c,d)$ % A*C-B*D$ addsq(multsq(a,c),negsq multsq(b,d))$ symbolic procedure vectprod expn$ begin scalar x,y,rnx,rny$ x:=tensval car expn$ y:=tensval cadr expn$ rnx:=tensrnk car expn$ rny:=tensrnk cadr expn$ if rnx=1 and rny=1 then return list(dimension!* - 2,antisymsum(x,y)) else if rnx=2 and rny=1 then return list(dimension!* - 1,mapcar(x,function(lambda a$ antisymsum(a,y)))) else if rnx=1 and rny=2 then return list(dimension!* - 1, if dimension!*=3 then tp1 copy2(mapcar(tp1 copy2(y,t), function(lambda a$ antisymsum(x,a))),t) else mapcar(tp1 copy2(y,t), function(lambda a$ antisymsum(x,a)) )) else lprie list(" VECTPROD of ",expn) end$ put('!?,'express,'vectprod)$ algebraic operator diff$ symbolic procedure gradexpress expn$ begin scalar arg,vt,ans,row,z$ arg:=tensval car expn$ vt:=tensrnk car expn$ if vt=0 then for each i in coordindx!* do ans:=simp!* list('quotient, list('diff, list('!*sq,arg,nil), getel list('coordinats,i)), getel list('sf,i)) . ans else if vt=1 then for each i in coordindx!* do <<row:=nil$ for each j in coordindx!* do <<z:=list('diff, listarsq(arg,j), getel list('coordinats,i))$ z:=list list('quotient, z, getel list('sf,i))$ for each k in coordindx!* do z:=list('times, getel list('christoffel,i,j,k), listarsq(arg,k)) . z$ row:=simp!*('plus.z) . row>>$ ans:=row . ans>> else lprie list(" GRAD of ",expn)$ return list(vt+1,ans) end$ put('grad,'express,'gradexpress)$ symbolic procedure divexpress expn$ begin scalar arg,vt,ans,z$ arg:=tensval car expn$ vt:=tensrnk car expn$ if vt=1 then <<for each i in coordindx!* do z:=list('diff, list('times, sfprod!*, listarsq(arg,i), list('recip, getel list('sf,i))), getel list('coordinats,i)) . z$ z:='plus . z$ z:=list('quotient, z, sfprod!*)$ ans:=simp!* z>> else if vt=2 then for each i in coordindx!* do <<z:=nil$ for j:=1:dimension!* do <<z:=list('quotient, list('diff, list('times, listarsq(listar(arg,j),i), sfprod!*, list('recip, getel list('sf,j))), getel list('coordinats,j)), sfprod!*) . z$ for l:=1:dimension!* do z:=list('times, listarsq(listar(arg,j),l), getel list('christoffel,j,i,l)) . z>>$ ans:=simp!*('plus.z) . ans>> else lprie list(" DIV of ",expn)$ return list(vt-1,ans) end$ put('div,'express,'divexpress)$ symbolic procedure laplexpress expn$ begin scalar arg,vt,ans$ arg:=tensval car expn$ vt:=tensrnk car expn$ if vt=0 then <<for i:=1:dimension!* do ans:=list('diff, list('times, sfprod!*, list('expt, getel list('sf,i), -2), list('diff, list('!*sq,arg,nil), getel list('coordinats,i))), getel list('coordinats,i)).ans$ ans:=list(0,simp!* list('quotient, 'plus . ans, sfprod!*))>> else if vt=1 then ans:=divexpress list gradexpress expn else lprie list(" LAPLACIAN of ",expn)$ return ans end$ put('lapl,'express,'laplexpress)$ symbolic procedure curlexpress expn$ begin scalar arg,vt,ans,ic1,ic2$ arg:=tensval car expn$ vt:=tensrnk car expn$ if vt=1 then for each i in (if dimension!* = 3 then coordindx!* else '(1) ) do <<ic1:=listar(cyclic1!*,i)$ ic2:=listar(cyclic2!*,i)$ ans:=simp!* list('times, getel list('sf,i), list('recip,sfprod!*), list('difference, list('diff, list('times, getel list('sf,ic2), listarsq(arg,ic2)), getel list('coordinats,ic1)), list('diff, list('times, getel list('sf,ic1), listarsq(arg,ic1)), getel list('coordinats,ic2)))) . ans>> else lprie list(" CURL of ",expn)$ return (if dimension!* = 3 then list(1,ans) else list(0,car ans)) end$ put('curl,'express,'curlexpress)$ flag('(cons grad div lapl curl tens vect dyad dirdf !& !# !?) ,'tensfn)$ symbolic procedure exscalval u$ begin scalar fce,args$ fce:=car u$ args:=mapcar(cdr u,function(lambda a$cddr express a))$ fce:=eval(get(fce,'express) . list mkquote args)$ if car fce=0 then return cadr fce else typerr(u," is not scalar ")$ return (nil . 1) end$ algebraic$ infix #,?,&$ precedence .,**$ precedence #,.$ precedence ?,#$ precedence &,?$ symbolic flag('(cons !# !? div lapl curl dirdf),'full)$ symbolic for each a in '(cons !# !? div lapl curl dirdf) do put(a,'simpfn,'exscalval)$ symbolic procedure scalefactors transf$ begin scalar var$ dimension!*:=car transf$ transf:=cdr transf$ if dimension!*=2 then <<var:=cddr transf$ transf:=list( car transf,cadr transf)>> else if dimension!*:=3 then <<var:=cdddr transf$ transf:=list(car transf,cadr transf,caddr transf)>> else lprie list(" Can't handle dimension = ",dimension!*)$ if dimension!*=length var then t else lprie list(" Transformation ",transf,var)$ for i:=1:dimension!* do setel(list('coordinats,i),listar(var,i))$ for row:=1:dimension!* do for col:=1:dimension!* do setel(list('jacobian,row,col), aeval list('df, listar(transf,col), getel list('coordinats ,row)))$ updatedimen()$ rscale() end$ deflist('((scalefactors rlis)),'stat)$ array jacobian(3,3),coordinats (3),sf(3),christoffel(3,3,3)$ procedure rscale()$ begin sfprod!*:=1$ nscal!*:=nscal!* + 1$ for row:=1:dimension!* do <<for col:=1:(row-1) do <<sf(row):=gcov(row,col)$ if sf(row)=0 then nil else write " WARNING: Coordinate system is nonorthogonal" ," unless following simplifies to zero: ", sf(row) >>$ write sf(row):=sqrt gcov(row,row)$ sfprod!*:=sfprod!* *sf(row)>>$ on nero$ for i:=1:dimension!* do for j:=1:dimension!* do for k:=1:dimension!* do begin christoffel(i,j,k):= ((if i=j then df(sf(j),coordinats (k)) else 0) -(if i=k then df(sf(k),coordinats (j)) else 0)) /(sf(j)*sf(k))$ if wrchri(a)=0 then write christoffel(i,j,k):= christoffel(i,j,k) end$ off nero end$ procedure gcov(j,k)$ for l:=1:dimension!* sum jacobian(j,l)*jacobian(k,l)$ symbolic$ symbolic procedure simpwrchri u$ if !*wrchri then nil . 1 else 1 . 1$ put('wrchri,'simpfn,'simpwrchri)$ symbolic procedure rmat$ 'dyad . cdr matstat()$ symbolic procedure formdyad(u,v,m)$ 'list . mkquote 'dyad . cddr formmat(u,v,m)$ put('dyad,'stat,'rmat)$ put('dyad,'formfn,'formdyad)$ symbolic procedure dirdfexpress expn$ begin scalar arg,vt,direc,ans,z,dj,di,argj,sfj,sfi,cooj$ arg:=cadr expn$ vt:=tensrnk arg$ direc:=car expn$ if not (tensrnk direc=1) then return lprie list (" Direction in DIRDF is not a vector ",expn)$ if vt=0 then return inprodexpress list (direc, gradexpress list arg)$ arg:=tensval arg$ direc:=tensval direc$ if not vt=1 then return lprie list (" Argument of DIRDF is dyadic ",expn)$ for each i in coordindx!* do <<z:=nil$ di:=listarsq(direc,i)$ sfi:=getel list('sf,i)$ for j:=1:dimension!* do <<dj:=listarsq(direc,j)$ argj:=listarsq(arg,j)$ sfj:=getel list('sf,j)$ cooj:=getel list('coordinats,j)$ z:=list('times, dj, list('recip,sfj), list('diff, listarsq(arg,i), cooj)) . z$ z:=list('times, di, argj, list('recip,sfi), list('recip,sfj), list('df,sfi,cooj)) . z$ z:=list('minus, list('times, dj, argj, list('recip,sfi), list('recip,sfj), list('df, sfj, getel list('coordinats,i)))) . z>>$ z:='plus . z$ z:=simp!* z$ ans:=z . ans >>$ return list(1,ans) end$ put('dirdf,'express,'dirdfexpress)$ symbolic procedure dfexpress expn$ begin scalar arg,vt,rest$ arg:=tensval car expn$ vt:=tensrnk car expn$ rest:=cdr expn$ rest:=mapcar(rest,function(lambda a$ if tensrnk a=0 then if atom tensval a then tensval a else if cdr tensval a=1 and numberp car tensval a then car tensval a else !*q2k tensval a else lprie list(" Bad arg of DF ", expn)))$ if vt=0 then return list(0,simpdf(list('!*sq,arg,t) . rest)) else if vt=1 then return list(1,mapcar(arg,function( lambda a$simpdf(list('!*sq,a,t) . rest))) ) else if vt=2 then return list(2,mapcar(arg,function( lambda a$mapcar(a,function( lambda b$simpdf(list('!*sq,b,t) . rest))) ))) else lprie list(" Bad tensor in DF ",expn) end$ put('df,'express,'dfexpress)$ symbolic procedure diffexpress expn$ begin scalar arg,vt,rest$ arg:=tensval car expn$ vt:=tensrnk car expn$ rest:=cdr expn$ rest:=mapcar(rest,function(lambda a$ if tensrnk a=0 then if atom tensval a then tensval a else if cdr tensval a=1 and numberp car tensval a then car tensval a else !*q2k tensval a else lprie list(" Bad arg of DIFF ", expn)))$ if vt=0 then return list(0,simp('diff . (prepsq arg . rest))) else if vt=1 then return list(1,mapcar(arg,function( lambda a$simp('diff . (prepsq a . rest))) )) else if vt=2 then return list(2,mapcar(arg,function( lambda a$mapcar(a,function( lambda b$simp('diff . (prepsq b . rest))) )))) else lprie list(" Bad tensor in DIFF ",expn) end$ put('diff,'express,'diffexpress)$ algebraic$ scalefactors 3,x,y,z,x,y,z$ endmodule$ %*********************************************************************** %***** ***** %***** M O D U L E I I M E T ***** %***** ***** %*********************************************************************** module iimet$ % Author: R. Liska % Program IIMET, Version REDUCE 3.4 05/1991$ symbolic$ global '(cursym!* !*val dimension!*)$ fluid '(!*exp alglist!*)$ symbolic procedure array u$ begin scalar msg,erfg$ msg:=!*msg$ !*msg:=nil$ erfg:=erfg!*$ erfg!*:=nil$ arrayfn('symbolic, mapcar (u,function(lambda a$car a . sub1lis cdr a)))$ erfg!*:=erfg$ !*msg:=msg end$ symbolic procedure sub1lis u$ if null u then nil else ((car u - 1) . sub1lis cdr u)$ sfprod!*:=1$ global'(date!*!*)$ date!*!*:= "IIMET Ver 1.1 17/05/91"$ put('version,'stat,'rlis)$ put('diff,'simpfn,'simpiden)$ global '(coords!* icoords!* dvars!* grids!* given!* same!* difml!* iobjs!* !*twogrid !*eqfu !*fulleq !*centergrid)$ switch twogrid,eqfu,fulleq,centergrid$ !*twogrid:=t$ % Given functions can be on both grids. !*eqfu:=nil$ % During pattern matching the given and % looked for functions are different. !*fulleq:=t$ % Optimalization is performed on both sides of PDE. !*centergrid:=t$ % Centers of grid cells are in points I % (otherwise in I+1/2). icoords!*:='(i j k l m n i1 j1 k1 l1 m1 n1)$ % Indices which are given implicit. procedure coordfn$ % Stat procedure of the COORDINATES statement, which defines indexes % of coordinates. begin scalar cor,icor$ flag('(into),'delim)$ cor:=remcomma xread nil$ remflag('(into),'delim)$ if cursym!* eq 'into then icor:=remcomma xread nil else if cursym!* eq '!*semicol!* then icor:=icoords!* else return symerr('coordfn,t)$ return list('putcor, mkquote cor, mkquote icor) end$ put('coordinates,'stat,'coordfn)$ flag('(putcor),'nochange)$ procedure putcor(u,v)$ begin scalar j$ j:=1$ coords!*:=u$ while u do <<if not idp car u then msgpri (" Coordinate ",car u," must be identifier",nil,'hold)$ if not(idp car v or fixp car v) then msgpri (" Index ",car v," must be identifier or integer",nil,'hold)$ put(car u,'index,list car v)$ put(car v,'coord,list car u)$ put(car u,'ngrid,j)$ j:=j+1$ put(car u,'simpfn,'simpiden)$ u:=cdr u$ v:=cdr v >> end$ procedure tcar u$ if pairp u then car u else u$ procedure grid u$ % Procedure definning the statement GRID. eval list(get(car u,'grid), mkquote cdr u)$ put('grid,'stat,'rlis)$ put('uniform,'grid,'gridunif)$ procedure gridunif u$ flag(u,'uniform)$ procedure dependence u$ % Procedure definning the statemnt DEPENDENCE. begin scalar x,y,z,gg,l,te,yy,y1,yl$ if null coords!* then rederr " Coordinates have not been defined yet"$ gg:=explode '!*grid$ l:=list(length coords!* + 1)$ a:x:=car u$ y:=car x$ if idp y then if not y memq dvars!* then dvars!*:=y . dvars!* else nil else return msgpri(" Variable ",y," must be identifier",nil, 'hold)$ z:=cdr x$ x:=car z$ if not numberp x then go to b$ if x=1 then apply('vectors,list y) else if x=2 then apply('dyads,list y) else if x=0 then t else return errpri2(car u,'hold)$ z:=cdr z$ b:yl:=nil$ yy:=explode y$ te:=aeval y$ if eqcar(te,'tensor) then te:=caddr te else te:=nil$ if te=1 then for i:=1:dimension!* do <<y1:=intern compress append(yy,explode i)$ yl:=y1 . yl$ setk1(list(y,i),y1,t)$ x:=intern compress append(explode y1,gg)$ % The name of an array filled with grids of Y(I) put(y1,'grid,x)$ eval list('array,mkquote list(x . l)) >> else if te=2 then for i:=1:dimension!* do for j:=1:dimension!* do <<y1:=intern compress append(yy,append(explode i, explode j))$ yl:=y1 . yl$ setk1(list(y,i,j),y1,t)$ x:=intern compress append(explode y1,gg)$ % The name of an array filled with grids of Y(I) put(y1,'grid,x)$ eval list('array,mkquote list(x . l)) >> else <<yl:=y . yl$ x:=intern compress append(yy,gg)$ put(y,'grid,x)$ eval list('array,mkquote list(x . l)) >>$ for each a in yl do put(a,'simpfn, 'simpiden)$ put(y,'names,reverse yl)$ if te member '(1 2) then <<put(y,'kkvalue,get(y,'kvalue))$ remprop(y,'kvalue) >>$ for each v in z do if v memq coords!* then for each w in yl do depend1(w,v,t) else msgpri(" Identifier ",v," is not coordinate",nil,'hold)$ u:=cdr u$ if u then go to a$ return nil end$ put('dependence,'stat,'rlis)$ procedure given u$ begin scalar x,xnam$ a:x:=car u$ xnam:=get(x,'names)$ if not idp x then msgpri (" Variable ",x," must be identifier",nil,'hold) else if xnam then given!* := union(xnam,given!*) else msgpri (" Identifier ",x," is not variable",nil,'hold)$ u:=cdr u$ if u then go to a$ return nil end$ put('given,'stat,'rlis)$ procedure cleargiven$ <<remflag(given!*,'noeq)$ given!*:=nil >>$ put('cleargiven,'stat,'endstat)$ flag('(cleargiven),'eval)$ newtok'(( !. !. ) isgr)$ algebraic infix ..$ grids!* := '(one half)$ procedure trlis$ % Stat procedure of the statement ISGRID. begin scalar x$ put('!*,'newnam,'tims)$ x:=rlis()$ remprop('!*,'newnam)$ return x end$ procedure formtr(u,vars,mode)$ list('isgrid,mkquote cdr u)$ procedure isgrid u$ % Proceudre definning the statement ISGRID. begin scalar x,y,z,z1,te,gd,lz,lz1$ a:x:=car u$ y:=car x$ x:=cdr x$ if not(y memq dvars!*) then return msgpri (" Identifier ",y," is not variable",nil,'hold)$ if null x then go to er$ te:=aeval y$ te:=if eqcar(te,'tensor) then caddr te else nil$ if (te=1 and null atom x and indexp car x and gridp cdr x) or (te=2 and null atom x and indexp car x and null atom cdr x and indexp cadr x and gridp cddr x) or ((te=0 or null te) and null atom x and gridp x) then t else go to er$ if te=1 then <<z:=car x$ x:=cdr x$ lz:=nil$ if numberp z then lz:=z . lz else for i:=1:dimension!* do lz:=i . lz$ for each a in lz do mapc(x,function(lambda b$ setel(list(get(cadr assoc(list(y,a), get(y,'kkvalue)), 'grid), car b), cadr b . nil))) >> else if te=2 then <<z:=car x$ z1:=cadr x$ x:=cddr x$ lz:=lz1:=nil$ if numberp z then lz:=z . lz else for i:=1:dimension!* do lz:=i . lz$ if numberp z1 then lz1:=z1 . lz1 else for i:=1:dimension!* do lz1:=i . lz1$ for each a in lz do for each b in lz1 do mapc(x,function(lambda c$ setel(list(get(cadr assoc(list(y,a,b), get(y,'kkvalue)), 'grid), car c), cadr c . nil))) >> else mapc(x,function(lambda c$ setel(list(get(y,'grid),car c),cadr c . nil)))$ u:=cdr u$ if u then go to a$ return nil$ er:errpri2(car u,'hold) end$ put('isgrid,'stat,'trlis)$ put('isgrid,'formfn,'formtr)$ procedure indexp u$ u eq 'tims or (numberp u and 0<u and u<dimension!*)$ procedure gridp u$ null atom u and begin scalar x$ a:x:=car u$ if null atom x and car x eq 'isgr and null atom cdr x and cadr x memq coords!* and null atom cddr x and caddr x memq grids!* then rplaca(u,cdr x) else return nil$ x:=car u$ rplaca(x,get(car x,'ngrid))$ u:=cdr u$ if u then go to a$ return t end$ procedure grideq u$ begin scalar x,y,z$ x:=u$ a:y:=car x$ if not car y memq dvars!* then return msgpri( "Identifier ",car y," is not variable",nil,'hold)$ z:=cdr y$ b:if not atom car z and caar z eq 'isgr and cadar z memq coords!* and caddar z memq grids!* then put(car y,cadar z,caddar z) else return errpri2(u,'hold)$ z:=cdr z$ if z then go to b$ x:=cdr x$ if x then go to a$ return nil end$ put('grideq,'stat,'rlis)$ fluid '(coord unvars sunvars interpolp novars ncor nvar intp icor gvar hi hip1 hip2 him1 him2 lhs rhs lsun lun xxgrid resar)$ % GVAR - grid on which is actual equation integrated % NVAR - number of actual equation in IIM21 % NCOR - number of actual coordinate % UNVARS - list of looked for functions % NOVARS - list of functions controlled by the SAME statement % SUNVARS - list of looked for functions, which are not controlled by % SAME, and of given functions, which are not controlled by % SAME and which can be only on one grid (if OFF TWOGRID) % in case ON TWOGRID given functions are not in SUNVARS % distinguishing of given functions defined in ISGRID is done % in procedures TWOGRIDP and GETGRID % LSUN - length of SUNVARS-1 % INTERPOLP-flag for MATCHLINEAR (for OFF EXP), has value T for calling % from INTERPOL and NIL for calling from SIMPINTT - for % NGETVAR in SIMPINTT % RESAR - the name of an array which will be filled by the result procedure zero u$ nil . 1$ procedure ngetvar(u,v)$ if atom u then get(u,v) else if atom car u then get(car u,v) else nil$ procedure ngrid u$ if u eq 'one then 'none else if u eq 'half then 'nhalf else nil$ procedure gnot u$ % Procedure inverts denotation of half-integer and integer grid if u='one then 'half else if u='half then 'one else nil$ procedure same u$ begin scalar x,y,z,bo$ if null same!* then <<same!*:=list u$ return nil >>$ x:=same!*$ a:y:=car x$ z:=u$ bo:=nil$ while z and not bo do <<if car z memq y then bo:=t$ z:=cdr z >>$ if bo then go to b$ x:=cdr x$ if x then go to a$ same!*:= u . same!*$ return nil$ b:rplaca(x,union(y,u))$ return nil end$ put('same,'stat,'rlis)$ procedure clearsame$ same!*:=nil$ put('clearsame,'stat,'endstat)$ flag('(clearsame),'eval)$ procedure mksame$ begin scalar x,y,z,yy,bo$ x:=expndsame()$ a:y:=car x$ yy:=y$ while yy and not(car yy memq unvars) do yy:=cdr yy$ if null yy then <<msgpri (" Same declaration ",nil,list(y, " doesn't contain unknown variable"),nil,'hold)$ go to b >>$ if y neq yy then <<z:=car y$ % The looked for function on the first place rplaca(y,car yy)$ rplaca(yy,z) >>$ z:=car y$ yy:=cdr y$ put(z,'sames,yy)$ novars:=union(novars,yy)$ for each a in cdr y do % Testing if A has appeared in the statement DEPENDENCE if not get(a,'grid) then msgpri (" Identifier ",a," is not variable",nil,'hold) else put(a,'same,z)$ for i:=1:length coords!* do <<yy:=y$ bo:=nil$ while yy and not bo do % Test on ISGRID <<bo:=getel list(get(car yy,'grid),i)$ yy:=cdr yy >>$ if bo then filgrid(y,bo,i) >>$ b:x:=cdr x$ if x then go to a$ sunvars:=setdiff(unvars,novars)$ return unvars end$ procedure filgrid(y,bo,i)$ % Filling up after finding ISGRID according to SAME begin scalar yy,bg$ yy:=y$ while yy do <<bg:=getel list(get(car yy,'grid),i)$ if null bg then setel(list(get(car yy,'grid),i),bo) else if bg eq bo then t else msgpri (" Same declaration ",nil,list(y, " doesn't correspond to isgrid declarations"),nil,'hold)$ yy:=cdr yy>> end$ procedure expndsame$ % Extending SAME!* by new identifiers for vectors and tensors begin scalar x,y,sam$ x:=same!*$ a:y:=mapcan(car x,function(lambda a$copy1 get(a,'names)))$ sam:=y . sam$ x:=cdr x$ if x then go to a$ return sam end$ procedure copy1 u$ if null u then nil else if atom u then u else car u . copy1 cdr u$ procedure nrsame$ % Changing the numbering of variables according to SAME for each a in sunvars do begin scalar x,nx$ x:=get(a,'sames)$ if x then <<nx:=get(a,'nrvar)$ for each b in x do put(b,'nrvar,nx) >>$ return nil end$ procedure iim u$ % Procedure defines the statement IIM begin scalar xx,xxx,be,beb1,val,twogr$ iim1 u$ iobjs!*:=append(unvars,append(coords!*,given!*))$ val:=!*val$ !*val:=nil$ novars:=sunvars:=nil$ if same!* then mksame() else sunvars:=unvars$ twogr:=!*twogrid$ xxx:=setdiff(given!*,novars)$ if !*twogrid then if null xxx then !*twogrid:=nil else flag(xxx,'twogrid) else sunvars:=union(sunvars,xxx)$ flag(given!*,'noeq)$ xxx:=0$ % Numbering of variables and equation for each a in sunvars do <<put(a,'nrvar,xxx)$ xxx:=xxx+1 >>$ if same!* then nrsame()$ xxx:=0$ for each a in unvars do <<put(a,'nreq,xxx)$ xxx:=xxx+1 >>$ lun:=length unvars-1$ lsun:=length sunvars-1$ eval list('array,mkquote list('!*f2 . add1lis list(lun,lsun,1)))$ xxx:=coords!*$ d:coord:=car xxx$ icor:=tcar get(coord,'index)$ difml!*:=nil$ for i:=0:10 do difml!*:=append(difml!*, for each a in getel list('difm!*,i) collect if (xx:=atsoc(coord,cdr a)) then car a . cdr xx else if (xx:=atsoc('all,cdr a)) then car a . cdr xx else nil )$ difml!*:=mapcon(difml!*,function(lambda a$if null car a then nil else list car a ))$ if !*twogrid then difml!*:= for each a in difml!* collect if (xx:=caadr a) and (!*eqfu or memq(caar xx,'(f g))) then (car a . extdif(cdr a,nil)) else a$ be:=iim2 ()$ iim21 be$ if car be then beb1:=iim22 be else beb1:=list(car be,cadr be,car be)$ if not fixp intp then msgpri(" INTP after heuristic search ", nil,list("is not a number, INTP=",intp),nil,nil)$ if not(intp=0) then iim3 beb1$ iim4 ()$ xxx:=cdr xxx$ if xxx then go to d$ iim5 ()$ for each a in '(rtype avalue dimension) do remprop('!*f2,a)$ !*val:=val$ !*twogrid:=twogr$ return nil end$ procedure extdif(x,lg)$ % Performs corrections of diff. schemes for given functions on % two grids - everytime chooses the scheme with minimal penalty. % LG - list of all terms from (U V W F G), which has been in X % already changed and choosen. begin scalar olds,news,y,gy,xx,lgrid,gg,g1$ lgrid:=get('difm!*,'grids)$ gy:=caar x$ gg:=gy$ for each a in lg do gg:=delete(atsoc(a,gg),gg)$ if gg then gg:=caar gg else return x$ x:=mapcar(x,function(lambda a$a))$ a:xx:=x$ y:=car x$ gy:=car y$ g1:=atsoc(gg,gy)$ gy:=(gg . gnot cdr g1) . delete(g1,gy)$ gy:=acmemb(gy,lgrid)$ while cdr xx and not eqcar(cadr xx,gy) do xx:=cdr xx$ if cdr xx then <<olds:=y . (cadr xx . olds)$ y:=if cadr y > cadadr xx or (cadr y=cadadr xx and sublength caddr y > sublength car cddadr xx) then cadr xx else y$ rplacd(xx,cddr xx) >> else olds:=y . olds$ gy:=car y$ g1:=atsoc(gg,gy)$ gy:=delete(g1,gy)$ if null gy then t else if (xx:=acmemb(gy,lgrid)) then gy:=xx else nconc(lgrid, list gy)$ y:=gy . cdr y$ news:=y . news$ x:=cdr x$ if x then go to a$ if(xx:=caar news) and (!*eqfu or memq(caar xx,'(f g))) then <<olds:=extdif(olds,gg . lg)$ news:=extdif(news,lg) >>$ return nconc(olds,news) end$ procedure sublength u$ if atom u then 0 else length u + sublengthca u$ procedure sublengthca u$ if null u then 0 else sublength car u + sublengthca cdr u$ procedure iim1 u$ % Checks the syntax of the IIM statement, calculates scalar PDEs, % vector and tensor equations are expanded to scalar components. begin scalar x,xx,e,te,exp$ terpri()$ prin2t"*****************************"$ prin2 "***** Program ***** "$ prin2t date!*!*$ prin2t"*****************************"$ exp:=!*exp$ !*exp:=t$ rhs:=lhs:=unvars:=nil$ if null coords!* then return lprie " Coordinates defined not yet"$ if null dvars!* then return lprie " Variables defined not yet"$ for each v in dvars!* do if eqcar((te:=aeval v),'tensor) and caddr te member '(1 2) then put(v,'kvalue,get(v,'kkvalue))$ if atom u or not idp car u then return errpri2(u,'hold)$ resar:=car u$ u:=cdr u$ a:if atom u or atom cdr u then return errpri2(u,'hold)$ x:=car u$ if not idp x then return msgpri (" Parameter ",x," must be identifier",nil,'hold) else if not(x memq dvars!*) then return msgpri (" Identifier ",x," is not variable",nil,'hold) else if x memq unvars then return msgpri (" Variable ",x," has second appearance",nil,'hold) else if x memq given!* then return msgpri (" Variable ",x," cannot be declared given",nil,'hold) else unvars:=x . unvars$ e:=cadr u$ if not eqexpr e then return msgpri (" Parameter ",e," must be equation",nil,'hold) else e:=aeval list('times, list('difference,cadr e,caddr e), sfprod!*)$ if atom e then return msgpri (" Equation ",e," must be P.D.E.",nil,'hold)$ te:=aeval x$ te:=if eqcar(te,'tensor) then caddr te else nil$ if(te=1 and car e eq 'tensor and caddr e=1) or (te=2 and car e eq 'tensor and caddr e=2) or (null te and car e eq 'tensor and caddr e=0) then e:=cadddr e % Necessary to carrect after change in EXPRESS else if null te and car e eq '!*sq then e:=cadr e else return msgpri (" Tensor order of",x," does not correspond to order of ",e, 'hold)$ lhs:=e . lhs$ u:=cddr u$ if u then go to a$ for each v in dvars!* do if eqcar((te:=aeval v),'tensor) and caddr te member '(1 2) then remprop(v,'kvalue)$ b:x:=car unvars$ e:=car lhs$ % Transformation of vectors and tensor into components te:=aeval x$ te:=if eqcar(te,'tensor) then caddr te else nil$ if te=1 then rhs:=append(e,rhs) else if te=2 then for each a in reverse e do rhs:=append(a,rhs) else rhs:=e . rhs$ xx:=append(get(x,'names),xx)$ % Add the checking if given equation unvars:=cdr unvars$ % solves given variable (tensor var.) lhs:=cdr lhs$ if unvars then go to b$ unvars:=xx$ lhs:=rhs$ put('diff,'simpfn,'zero)$ % Splitting left and right hand side alglist!*:=nil . nil$ % All derivatives go to left h.s. rhs:=mapcar(rhs,function resimp)$ put('diff,'simpfn,'simpiden)$ alglist!*:=nil . nil$ x:=lhs$ xx:=rhs$ terpri()$ prin2t " Partial Differential Equations"$ prin2t " =============================="$ terpri()$ c:rplaca(xx,negsq car xx)$ rplaca(x,prepsq!* addsq(car x,car xx))$ rplaca(xx,prepsq!* car xx)$ maprin car x$ prin2!* " = "$ maprin car xx$ terpri!* t$ rplaca(x,prepsq simp car x)$ rplaca(xx,prepsq simp car xx)$ x:=cdr x$ xx:=cdr xx$ if x then go to c$ terpri()$ x:=length lhs-1$ if x=0 then eval list('array, mkquote list(resar . add1lis list(1))) else eval list('array,mkquote list(resar . add1lis list(x,1)))$ !*exp:=exp$ return nil end$ procedure iim2$ % Defines the steps of the grid, splits variables to free and predefined % grid in actual coordinate. begin scalar b,e,xx,dihalf,dione,dihalfc$ e:=append(explode 'h,explode coord)$ e:=intern compress e$ if flagp(coord,'uniform) then hi:=hip1:=him1:=him2:=e else <<put(e,'simpfn,'simpiden)$ xx:=tcar get(coord,'index)$ him1:=list(e,list('difference,xx,1))$ him2:=list(e,list('difference,xx,2))$ hi:=list(e,xx)$ hip2:=list(e,list('plus,xx,2))$ hip1:=list(e,list('plus,xx,1)) >>$ dihalf:=list( 'di . list('quotient, list('plus,him1,hi), 2), 'dim1 . him1, 'dip1 . hi, 'dim2 . list('quotient, list('plus,him2,him1), 2), 'dip2 . list('quotient, list('plus,hi,hip1), 2))$ dihalfc:=list( 'di . list('quotient, list('plus,hip1,hi), 2), 'dim1 . hi, 'dip1 . hip1, 'dim2 . list('quotient, list('plus,hi,him1), 2), 'dip2 . list('quotient, list('plus,hip2,hip1), 2))$ dione:=list( 'di . hi, 'dim1 . list('quotient, list('plus,him1,hi), 2), 'dip1 . list('quotient, list('plus,hi,hip1), 2), 'dim2 . him1, 'dip2 . hip1)$ put('steps,'one, ('i . icor) . (if !*centergrid then dione else dihalf))$ put('steps,'half, ('i . list('plus, icor, '(quotient 1 2))) . (if !*centergrid then dihalfc else dione))$ ncor:=get(coord,'ngrid)$ % Number of the COODR coordinate e:=nil$ for each a in sunvars do % Splitting of variables with predefined % grid. if (xx:=getel list(get(a,'grid),ncor)) and car xx then e:=a . e$ b:=setdiff(sunvars,e)$ return list(b,e) end$ procedure filfree(var,vgrid,freelst,pgr,peq)$ begin scalar x,nx,grn,nv,ng,ngrn,g1,g2,saml,bsam,asam,egrid$ x:=ngetvar (var,'nrvar)$ c:put(var,pgr,vgrid)$ egrid:=vgrid$ if flagp(var,'noeq) then go to d$ nx:=ngetvar (var,'nreq)$ % calulating in which point will be the euation for VAR integrated if egrid:=get(var,coord) then go to a else egrid:=vgrid$ put('f2val,'free,'f2vzero)$ if (g1:=f2eval(nx,x,0)) > (g2:=f2eval(nx,x,1)) then egrid:=gnot vgrid$ if not g1=g2 then go to a$ put('f2val,'free,'f2vmin)$ if(g1:=f2eval(nx,x,0)) > (g2:=f2eval(nx,x,1)) then egrid:=gnot vgrid$ if not g1=g2 then go to a$ put('f2val,'free,'f2vmax)$ if (g1:=f2eval(nx,x,0)) > (g2:=f2eval(nx,x,1)) then egrid:=gnot vgrid$ a:put(var,peq,egrid)$ % Penalties for free variables in the equation for VAR grn:=gnot egrid$ ng:=ngrid egrid$ ngrn:=ngrid grn$ for each a in freelst do <<nv:=ngetvar (a,'nrvar)$ asam:=a . get(a,'sames)$ for each aa in asam do put(aa,pgr,egrid)$ put(a,ng,cfplus2(get(a,ng),getel list('!*f2,nx,nv,0)))$ for each aa in asam do put(aa,pgr,grn)$ put(a,ngrn,cfplus2(get(a,ngrn),getel list('!*f2,nx,nv,1)))$ for each aa in asam do remprop(aa,pgr) >>$ if bsam then go to d$ saml:=get(var,'sames)$ bsam:=t$ d:if null saml then go to b$ var:=car saml$ saml:=cdr saml$ go to c$ b:return egrid end$ procedure f2eval(i,j,k)$ eval getel list('!*f2,i,j,k)$ procedure f2plus(i,j,k,l)$ % Procedure fills F2(I,J,K) with the number F2(I,J,K)+L$ begin scalar ma,x,y$ if pairp l then if length car l=2 and cadr l=caddr l then l:=cadr l else if length l=3 and cadr l=caddr l and cadr l=cadddr l and cadr l=car cddddr l then l:=cadr l$ ma:=list('!*f2,i,j,k)$ x:=getel ma$ if numberp l then if numberp x then setel(ma,x+l) else rplaca(x,car x+l) else if numberp x then setel(ma,list(x,l)) else if (y:=assoc(car l,cdr x)) then <<rplaca(cdr y,cadr y+cadr l)$ rplaca(cddr y,caddr y+caddr l)$ if cddar l then <<rplaca(cdddr y,cadddr y + cadddr l)$ rplaca(cddddr y,car cddddr y+car cddddr l)>>>> else rplacd(x,(l . cdr x)) end$ procedure f2var u$ % Forms the elements of array !*F2 into the form % (FPLUS <CISLO> {(F2VAL U V N1 N2)}) if numberp u then u else ('fplus . car u . mapcar(cdr u,function(lambda a$ list('f2val,car a,cdr a))))$ macro procedure f2val x$ % Evaluates the expression (F2VAL ...) begin scalar us,ns,u,v,w,n1,n2,n3,n4,gu,gv,gw$ x:=cdr x; us:=car x$ ns:=cadr x$ u:=car us$ v:=cadr us$ n1:=car ns$ n2:= cadr ns$ gu:=get(u,xxgrid)$ gv:=get(v,xxgrid)$ if cddr us then <<w:=caddr us$ gw:=get(w,xxgrid)$ n3:=caddr ns$ n4:=cadddr ns>>$ return mkquote if w then if gu and gv and gw then if gu eq gv and gu eq gw then n1 else if gu eq gv then n2 else if gu eq gw then n3 else if gv eq gw then n4 else apply(get('f2val,'free),list(us,ns)) else if gu and gv then if gu eq gv then aplf2val(u,w,n1,n2) else aplf2val(u,w,n3,n4) else if gu and gw then if gu eq gw then aplf2val(u,v,n1,n3) else aplf2val(u,v,n2,n4) else if gv and gw then if gv eq gw then aplf2val(u,v,n1,n4) else aplf2val(u,v,n2,n3) else apply(get('f2val,'free),list(us,ns)) else if gu and gv then if gu eq gv then n1 else n2 else apply(get('f2val,'free),list(us,ns)) end$ procedure aplf2val(u,v,n1,n2)$ apply(get('f2val,'free),list(list(u,v),list(n1,n2)))$ macro procedure fplus u$ % Evaluates the expression (FPLUS ...) begin scalar x,y,z$ u:=cdr u; y:=car u$ a:u:=cdr u$ z:=eval car u$ if numberp z then y:=y+z else x:=z . x$ if cdr u then go to a$ return mkquote if x then ('fplus . y . x) else y end$ procedure cfplus2(u,v)$ % Adds the expressions of type (FPLUS ...). % Destroys U, does not destroy V. begin scalar f2v$ f2v:=get('f2val,'free)$ put('f2val,'free,'f2vunchange)$ if not fixp u then u:=eval u$ if not fixp v then v:=eval v$ put('f2val,'free,f2v)$ return if fixp u then if fixp v then (u + v) else ('fplus . (cadr v+u) . cddr v) else if numberp v then <<rplaca(cdr u,cadr u+v)$ u>> else <<rplaca(cdr u,cadr u+cadr v)$ nconc(u,cddr v) >> end$ procedure f2vunchange(us,ns)$ list('f2val,us,ns)$ procedure f2vzero(us,ns)$ 0$ procedure f2vplus(us,ns)$ eval('plus . ns)$ procedure f2vmax(us,ns)$ eval('max . ns)$ procedure f2vmin(us,ns)$ eval('min . ns)$ put('f2val,'fselect,'f2vplus)$ put('f2val,'fgrid,'f2vmin)$ procedure iim21 u$ % Fills the array !*F2 according to the system of PDE and penalties % given. % Fills the properties NONE,NHALF (FONE,FHALF) of free variables. % According to predefined variables filles the properties XGRID and EQ % of predefined variables. begin scalar b,e,lh,lhe,xx,rh,bdef$ b:=car u$ % Free vars e:=cadr u$ % Predefined vars for i:=0:lun do for j:=0:lsun do % Filling the array F2 <<setel(list('!*f2,i,j,0),0)$ setel(list('!*f2,i,j,1),0) >>$ nvar:=0$ % Number of actual variable lh:=lhs$ rh:=rhs$ interpolp:=nil$ put('intt,'simpfn,'simpiden)$ a:lhe:=car lh$ % Actual equation lhe:=formlnr list('intt,lhe,coord)$ rplaca(lh,lhe)$ bdef:=t$ for each var in sunvars do if get(var,coord) then t else bdef:=nil$ if null b and bdef then go to c$ % If there are no free variables it is not necessary to fill % the array F2 - no optimalization is necessary -> To use this % statement, we have to test if we know in which point (over % which interval) will all equation be integrated (discretized). put('intt,'simpfn,'simpintt)$ alglist!*:=nil . nil$ simp lhe$ put('intt,'simpfn,'simpiden)$ if !*fulleq then % Optimalizatioon is performed for both sides of <<lhe:=car rh$ % equations, otherwise only for left hand side. lhe:=formlnr list('intt,lhe,coord)$ put('intt,'simpfn,'simpintt)$ simp lhe$ alglist!*:=nil . nil$ put('intt,'simpfn,'simpiden)$ rh:=cdr rh >>$ c:nvar:=nvar+1$ lh:=cdr lh$ if lh then go to a$ for i:=0:lun do for j:=0:lsun do for k:=0:1 do setel(list('!*f2,i,j,k),f2var getel list('!*f2,i,j,k))$ xxgrid:='xgrid$ for each a in b do <<put(a,'none,0)$ put(a,'nhalf,0) >>$ for each a in e do % Predefined variables filfree(a,car getel list(get(a,'grid),ncor),b,'xgrid,'eq)$ % Predefined penalties intp:=0$ for each a in e do if a memq unvars then <<xx:=ngetvar(a,'nreq)$ for each c in e do if get(c,'xgrid) eq get(a,'eq) then intpfplus(xx,c,0) else intpfplus(xx,c,1) >>$ for each a in b do <<put(a,'fone,get(a,'none))$ put(a,'fhalf,get(a,'nhalf)) >>$ return nil end$ procedure iim22 u$ begin scalar b,e,bb,b1,b2,x,xx,x1,nv,g1,g2$ b:=car u$ e:=cadr u$ bb:=b$ % Heuristic determination of grids for % variables from B f:x:=car bb$ % Chose the next variable X put('f2val,'free,get('f2val,'fselect))$ xx:=abs(eval get(x,'none)-eval get(x,'nhalf))$ b2:=cdr bb$ while b2 do if xx<(x1:=abs(eval get(car b2,'none)-eval get(car b2,'nhalf))) then <<x:=car b2$ xx:=x1$ b2:=cdr b2 >> else b2:=cdr b2$ b1:=x . b1$ % List of variables subsequently choosen from B bb:=delete(x,bb)$ % List of variables remaining in B put('f2val,'free,get('f2val,'fgrid))$ put(x,'xgrid,'one)$ g1:=eval get(x,'none)$ put(x,'xgrid,'half)$ g2:=eval get(x,'nhalf)$ if g1>g2 then xx:='half else xx:='one$ filfree(x,xx,bb,'xgrid,'eq)$ intpgplus(x,xx)$ for each ax in (x . get(x,'sames)) do if ax memq unvars then <<x1:=ngetvar(ax,'nreq)$ for each a in append(b1,e) do if get(a,'xgrid)=get(ax,'eq) then intpfplus(x1,a,0) else intpfplus(x1,a,1) >>$ if bb then go to f$ return list(b,e,b1) end$ procedure intpfplus(nx1,a,n)$ intp:=cfplus2(intp,getel list('!*f2,nx1,ngetvar(a,'nrvar),n))$ procedure intpgplus(a,ga)$ intp:=cfplus2(intp,get(a,ngrid ga))$ procedure iim3 u$ begin scalar b,e,b1,bb$ prin2t" Backtracking needed in grid optimalization"$ b:=car u$ % Free vars e:=cadr u$ % Predefined vars b1:=caddr u$ for each a in b do % Full search - bactracking <<put(a,'none,get(a,'fone))$ put(a,'nhalf,get(a,'fhalf)) >>$ xxgrid:='bxgrid$ nbxgrid(e,'bxgrid,'beq,'xgrid,'eq)$ put('f2val,'free,'f2vunchange)$ varyback(b1,nil)$ for each a in union(unvars,given!*) do <<remprop(a,'bxgrid)$ remprop(a,'beq) >>$ return nil end$ procedure nbxgrid(u,ng,ne,og,oe)$ for each a in u do for each b in (a . get(a,'sames)) do <<put(b,ng,get(b,og))$ put(b,ne,get(b,oe)) >>$ procedure varyback(bb,b1)$ % Performs full search of BB. B1 is B-BB. N is the number of % interpolations performed up to now. if null bb then begin scalar none,nhalf,n,eqg,i,j$ n:=0$ for each a in unvars do <<none:=nhalf:=0$ i:=get(a,'nreq)$ for each b in sunvars do <<j:=get(b,'nrvar)$ none:=none + if get(b,'bxgrid) eq 'one then f2eval(i,j,0) else f2eval(i,j,1)$ nhalf:=nhalf + if get(b,'bxgrid) eq 'half then f2eval(i,j,0) else f2eval(i,j,1) >>$ put(a,'beq,if (eqg:=get(a,coord)) then eqg else if none<=nhalf then 'one else 'half)$ n:=n + if eqg then if eqg eq 'one then none else if eqg eq 'half then nhalf else <<msgpri("VARYBACK ",eqg,nil,nil,'hold)$ 0>> else if none<=nhalf then none else nhalf >>$ if n<intp then <<nbxgrid(b1,'xgrid,'eq,'bxgrid,'beq)$ for each a in unvars do put(a,'eq,get(a,'beq))$ intp:=n >> end else if intp=0 then t else <<varb(bb,b1,'one)$ varb(bb,b1,'half) >>$ procedure varb(bb,b1,xx)$ % Subprocedure of VARYBACK procedure % In BB are temporary free variables % In B1 are temporary predefined variables (over BXGRID property) begin scalar x$ x:=car bb$ for each a in (x . get(x,'sames)) do put(a,'bxgrid,xx)$ return varyback(cdr bb,x . b1) end$ procedure iim4$ begin scalar lh,rh,x,lhe,var$ intp:=intp/6$ prin2 intp$ prin2 " interpolations are needed in "$ prin2 coord$ prin2t " coordinate"$ for each a in unvars do <<prin2" Equation for "$ prin2 a$ prin2" variable is integrated in "$ prin2 get(a,'eq)$ prin2t" grid point" >>$ interpolp:=t$ put('intt,'simpfn,'simpinterpol)$ lh:=lhs$ rh:=rhs$ x:=unvars$ j:var:=car x$ gvar:=get(var,'eq)$ lhe:=car lh$ alglist!*:=nil . nil$ lhe:=prepsq simp lhe$ rplaca(lh,lhe)$ lhe:=car rh$ lhe:=formlnr list('intt,lhe,coord)$ lhe:=prepsq simp lhe$ rplaca(rh,lhe)$ x:=cdr x$ lh:=cdr lh$ rh:=cdr rh$ if x then go to j$ put('intt,'simpfn,'simpiden)$ return lhs end$ procedure iim5$ begin scalar lh,rh,val,nreq,ar$ val:=!*val$ !*val:=nil$ for each a in union(union(unvars,sunvars),given!*) do <<remprop(a,'xgrid)$ remprop(a,'eq)$ remprop(a,'nreq)$ remprop(a,'nrvar)$ remprop(a,'sames)$ remprop(a,'none)$ remprop(a,'nhalf)$ remprop(a,'fone)$ remprop(a,'fhalf) >>$ remflag(given!*,'twogrid)$ remflag(given!*,'noeq)$ terpri()$ prin2t " Equations after Discretization Using IIM :"$ prin2t " =========================================="$ terpri()$ lh:=lhs$ rh:=rhs$ nreq:=0$ k:rplaca(lh,prepsq!* simp!* car lh)$ maprin car lh$ prin2!* " = "$ rplaca(rh,prepsq!* simp!* car rh)$ maprin car rh$ terpri!* t$ terpri()$ ar:=if null cdr lhs then list(resar,0) else list(resar,nreq,0)$ setel(ar,car lh)$ ar:=if null cdr lhs then list(resar,1) else list(resar,nreq,1)$ setel(ar,car rh)$ lh:=cdr lh$ rh:=cdr rh$ nreq:=nreq+1$ if lh then go to k$ !*val:=val$ return nil end$ put('iim,'stat,'rlis)$ array difm!*(10)$ procedure iscomposedof(x,objs,ops)$ if null x then nil else if atom x then if idp x then memq(x,objs) else if fixp x then t else nil else if idp car x and car x memq ops and cdr x then iscompos(cdr x,objs,ops) else nil$ procedure iscompos(x,objs,ops)$ if null x then t else if idp car x then car x memq objs and iscompos(cdr x,objs, ops) else if numberp car x then iscompos(cdr x,objs,ops) else if atom car x then nil else if idp caar x then caar x memq ops and cdar x and iscompos(cdar x,objs,ops) and iscompos(cdr x,objs,ops) else nil$ global'(difconst!* diffuncs!*)$ difconst!* := '(i n di dim1 dip1 dim2 dip2)$ diffuncs!*:=nil$ procedure difconst u$ difconst!* := append(u,difconst!*)$ put('difconst,'stat,'rlis)$ procedure diffunc u$ <<diffuncs!*:=append(u,diffuncs!*)$ for each a in u do put(a,'matcheq,'matchdfunc) >>$ put('diffunc,'stat,'rlis)$ procedure matchdfunc(u,v)$ begin scalar x,y$ return if null u and null v then list t else if null u or null v then nil else if (x:=matcheq(car u,car v)) and (y:=matchdfunc(cdr u,cdr v)) then union(x,y) else nil end$ procedure difmatch u$ begin scalar l,gds,gdsf,pl,x,dx,y,z,coor$ coor:=car u$ if not atom coor then go to er$ u:=cdr u$ x:=car u$ if not iscomposedof(x,'(u f x n v w g), append(diffuncs!*,'(diff times expt quotient recip)))then go to er$ x:=prepsq simp x$ l:=if atom x then 0 else length x$ x:=x . nil$ if null(y:=getel list('difm!*,l)) then setel(list('difm!*,l),list x) else if (z:=assoc(car x,y)) then x:=z else nconc(y,list x)$ y:=cdr u$ a:gds:=nil$ gdsf:=nil$ if not eqexpr car y then go to b$ a1:if not(cadar y memq '(u v w f g) and caddar y memq grids!*) then go to er$ if cadar y memq '(f g) then gdsf:=(cadar y . caddar y) . gdsf else gds:=(cadar y . caddar y) . gds$ y:=cdr y$ if null y then go to er$ if eqexpr car y then go to a1$ b:if not fixp car y then go to er$ pl:=car y$ y:=cdr y$ if null y then go to er$ if not iscomposedof(car y,difconst!*,append(diffuncs!*,'(u x f v w g plus minus difference times quotient recip expt)))then go to er$ dx:=car y$ y:=cdr y$ gds:=nconc(gdsf,gds)$ defdfmatch(x,gds,pl,list dx,coor)$ if y then go to a$ return nil$ er:errpri2(y,'hold) end$ procedure defdfmatch(x,gds,pl,dx,coor)$ begin scalar y,z,yy$ y:=get('difm!*,'grids)$ if null y then put('difm!*,'grids,list gds) else if null gds then t else if (z:=acmemb(gds,y)) then gds:=z else nconc(y,list gds)$ y:=cdr x$ if y then if (yy:=atsoc(coor,y)) then if (z:=assoc(gds,cdr yy)) then <<rplacd(z,pl . dx)$ msgpri(" Difference scheme for ",car x," redefined ", nil,nil) >> else nconc(cdr yy,list(gds . (pl . dx))) else nconc(y,list(coor . list(gds . (pl . dx)))) else rplacd(x,list(coor . list(gds . (pl . dx))))$ return y end$ deflist('((difmatch rlis) (cleardifmatch endstat)),'stat)$ procedure cleardifmatch$ for i:=0:10 do difm!*(i):=nil$ flag('(cleardifmatch),'eval)$ procedure acmemb(u,v)$ if null v then nil else if aceq(u,car v) then car v else acmemb(u,cdr v)$ procedure aceq(u,v)$ if null u then null v else if null v then nil else if car u member v then aceq(cdr u,delete(car u,v)) else nil$ procedure matcheq(u,v)$ if null u or null v then nil else if numberp u then if u=v then list t else nil else if atom u then begin scalar x$ x:=eval list(get(u,'matcheq),mkquote u,mkquote (if atom v then list v else v))$ return if x then x else if null !*exp and pairp v and car v memq '(plus difference) then matchlinear(u,v) else nil end else if atom v then nil else if atom car u and car u eq car v then eval list(get(car u,'matcheq),mkquote cdr u,mkquote cdr v) else if null !*exp and car v memq'(plus difference) and car u eq 'diff then matchlinear(u,v) else nil$ algebraic operator matchplus$ fluid'(uu vv)$ procedure matchlinear(u,v)$ % Construction for OFF EXP and second and next coordinates begin scalar x,uu,vv,alg$ if not atom u then return if car u eq 'diff then matchlindf(u,v) else if car u eq 'times then matchlintimes(u,v) else nil$ uu:=u$ vv:='first$ x:=formlnr list('matchplus,v,coord)$ put('matchplus,'simpfn,'matchp)$ alg:=alglist!*$ alglist!*:=nil . nil$ simp x$ alglist!*:=alg$ put('matchplus,'simpfn,'simpiden)$ return if vv then list(u . (if interpolp then v else vv)) else nil end$ procedure matchp y$ begin scalar x$ if null vv then return(nil . 1)$ x:=matcheq(uu,car y)$ if null x then return begin vv:=nil$ return(nil . 1) end$ if vv eq 'first then return begin vv:=cdar x$ return (nil . 1) end$ if mainvareq(vv,cdar x) then return (nil . 1)$ vv:=nil$ return(nil . 1) end$ unfluid '(uu vv)$ procedure mainvareq(x,y)$ if atom x then eq(x,y) else if car x memq iobjs!* then eq(car x,car y) else if car x memq '(diff expt) then (car y eq car x and mainvareq(cadr x,cadr y) and cddr x=cddr y) else nil$ procedure tlist x$ if atom x then list x else x$ procedure matchlindf(u,v)$ begin scalar x,y,b$ x:=mapcar(cdr v,function fsamedf)$ y:=cdar x$ if null y then return nil$ x:=for each a in x collect if y=cdr a then car a else b:=t$ if b then return nil$ x:=(car v . x) . y$ return matchdf(cdr u,x) end$ procedure fsamedf u$ begin scalar x$ return if atom u then nil . nil else if car u eq 'minus then <<x:=fsamedf cadr u$ list('minus,car x) . cdr x >> else if car u eq 'diff then cadr u . cddr u else if car u eq 'times then begin scalar y,z$ x:=cdr u$ a:if null x or y=t then go to b$ if numberp car x then z:=car x . z else if eqcar(car x,'diff) then <<if y then y:=t else y:=cddar x$ z:=cadar x . z >> else if depends(car x,coord) then y:=t else z:=car x . z$ x:=cdr x$ go to a$ b:return if y=t then nil . nil else ('times . z) . y end else nil . nil end$ procedure matchlintimes(u,v)$ begin scalar x,y,z$ y:=cadr v$ if eqcar(y,'times) then y:=cdr y else if eqcar(y,'minus) and eqcar(cadr y,'times) then y:= (-1) . cdadr y else return nil$ x:=for each a in cdr v collect if eqcar(a,'times) then <<y:=intersect(y,cdr a)$ a>> else if eqcar(a,'minus) and eqcar(cadr a,'times) then <<y:=intersect(y,cdadr a)$ 'times . ((-1) . cdadr a) >> else y:=nil$ if null y then return nil$ x:=for each a in x collect <<z:=setdiff(cdr a,y)$ if null cdr z then car z else 'times . z >>$ x:=car v . x$ return matchtimes(cdr u,x . y) end$ procedure intersect(u,v)$ if null u or null v then nil else if member(car u,v) then car u . intersect(cdr u,v) else intersect(cdr u,v)$ procedure matchu(u,v)$ if car v memq unvars or (!*eqfu and car v memq given!*) then list(u . v) else if car v eq 'diff and not coord memq cddr v and matcheq(u,tlist cadr v) then list(u . (car v . (tlist cadr v . cddr v))) else if car v eq 'times then % Product can be inside brackets or in DIFF begin scalar x,b1,vv$ x:=for each a in cdr v collect a$ % To allow RPLACA vv:=car v . x$ b1:=0$ while x and b1<2 do <<if depends(car x,coord) then <<b1:=b1+1$ if atom car x then rplaca(x,list car x) >>$ x:=cdr x >>$ return if b1=0 or b1>1 then nil else (u . vv) end else nil$ put('u,'matcheq,'matchu)$ put('v,'matcheq,'matchu)$ put('w,'matcheq,'matchu)$ procedure matchf(u,v)$ if car v memq given!* then list(u . v) else if car v eq 'diff and not coord memq cddr v and matchf(u,tlist cadr v) then list(u . (car v . (tlist cadr v . cddr v))) else nil$ put('f,'matcheq,'matchf)$ put('g,'matcheq,'matchf)$ procedure matchx(u,v)$ if car v eq coord then list t else nil$ put('x,'matcheq,'matchx)$ procedure matchtimes(u,v)$ begin scalar bool,bo,x,y,y1,asl$ x:=u$ a:y:=t . v$ d:bool:=nil$ while not bool and cdr y do <<y:=cdr y$ bool:=matcheq(car x,car y)>>$ if null bool then go to b$ bo:=bool$ c: if not atom bo and not atom car bo then y1:=atsoc(caar bo,asl) else y1 := nil$ if y1 and not y1=car bo then go to d$ bo:=cdr bo$ if bo then go to c$ v:=delete(car y,v)$ x:=cdr x$ asl:=union(bool,asl)$ if x then go to a$ if v then return nil$ return asl$ b:return if null cdr v and cdr x then if y:=matcheq('times . x,car v) then union(asl,y) else nil else nil end$ put('times,'matcheq,'matchtimes)$ procedure matchexpt(u,v)$ if fixp cadr u then if cadr u=cadr v then matcheq(car u,car v) else nil else if cadr u='n then begin scalar x$ x:=matcheq(car u,car v)$ return if x then (('n . cadr v) . x) else nil end else nil$ put('expt,'matcheq,'matchexpt)$ procedure matchquot(u,v)$ begin scalar man,mad$ return if(man:=matcheq(car u,car v)) and (mad:=matcheq(cadr u,cadr v)) then union(man,mad) else nil end$ put('quotient,'matcheq,'matchquot)$ procedure matchrecip(u,v)$ matcheq(car u,car v)$ put('recip,'matcheq,'matchrecip)$ procedure matchdf(u,v)$ begin scalar x,asl,y$ asl:=matcheq(car u,car v)$ if null asl then return nil$ y:=x:=append(cdr v,nil)$ while x and car x neq coord do x:=cdr x$ if null x then return nil else if null cddr u then if null cdr x or idp cadr x then go to df1 else return nil else if cdr x and caddr u=cadr x then t else return nil$ rplacd(x,cddr x)$ df1:y:=delete(coord,y)$ if null y or null interpolp then return asl % !!! Be aware !!! in mixed derivations of product else return list(car u . ('diff . (cdar asl . y))) end$ put('diff,'matcheq,'matchdf)$ procedure finddifm u$ begin scalar x,v,asl,eqfu,b,bfntwo,bftwo1$ eqfu:=!*eqfu$ if eqfu then !*eqfu:=nil$ a:x:=t . difml!*$ bftwo1:=bfntwo$ bfntwo:=nil$ if !*eqfu then b:=t$ while cdr x and not asl do <<x:=cdr x$ asl:=matcheq(caar x,u)$ % Test for given variables of type F, which have to be on one grid, % if choosed substitution from DIFMATCH contains definition of F % on grids. if asl then for each a in asl do if not atom a and car a memq'(f g) and (cadr a memq novars or (null !*twogrid and cadr a memq sunvars)) and null atsoc(car a,caadar x) then bfntwo:=t$ if bfntwo then asl:=nil >>$ !*eqfu:=eqfu$ if null asl then if null b and eqfu then go to a else go to nm$ return list(('x . coord) . delete(t,asl),cdar x)$ nm:if eqcar(u,'times) and null !*exp then <<!*exp:=t$ alglist!*:=nil . nil$ v:=prepsq simp u$ v:=formlnr list('intt,v,coord)$ !*exp:=nil$ if null atom v and car v memq'(plus difference) then return ('special . simp v) >>$ msgpri(" Matching of ",u," term not find ",nil,'hold)$ if bfntwo or bftwo1 then lprie(" Variable of type F not defined on grids in DIFMATCH")$ return nil end$ procedure tdifpair x$ % From CDR ATSOC(.,ASL) makes an atom - free variable <<if eqcar(x,'diff) then if eqcar(cadr x,'times) then <<x:=cdadr x$ while x and not depends(car x,coord) do x:=cdr x$ x:=tdifpair car x>> % patch else x:=cadr x$ if pairp x then x:=car x$ x >>$ procedure simpintt u$ begin scalar asl,agdsl,l,x,nv,y,x1,y1,nv1,n1,n2,nn1,nn2, x2,y2,nv2,n3,n4,n5,n6,lgrids,gds$ u:=prepsq simp car u$ if u=1 then go to r$ asl:=finddifm u$ if null asl or eqcar(asl,'special) then go to r$ agdsl:=cadr asl$ % List from DIFML!* asl:=car asl$ % ASOC. list of assignments gds:=caar agdsl$ l:=length gds$ if l=0 then go to r$ a:y:=caar gds$ x:=atsoc(y,asl)$ if null x then go to er1$ x:=tdifpair cdr x$ if !*twogrid and flagp(x,'twogrid) then if l=1 then go to r else <<gds:=cdr gds$ l:=l-1$ go to a>>$ nv:=ngetvar(x,'nrvar)$ if l=1 then go to l1 else go to l2$ l1:x:=assoc(list(y . 'one),agdsl)$ if null x then go to er2$ f2plus(nvar,nv,0,6*cadr x)$ x:=assoc(list(y . 'half),agdsl)$ if null x then go to er2$ f2plus(nvar,nv,1,6*cadr x)$ go to r$ l2:y1:=caadr gds$ x1:=atsoc(y1,asl)$ if null x1 then go to er1$ x1:=tdifpair cdr x1$ if !*twogrid and flagp(x1,'twogrid) then if l=2 then go to l1 else <<gds:=cdr gds$ l:=l-1$ go to l2 >>$ nv1:=ngetvar(x1,'nrvar)$ lgrids:=get('difm!*,'grids)$ if l=3 then go to l3 else if l>3 then go to er$ l20:n1:=atsoc(acmemb(list(y . 'one,y1 . 'one),lgrids),agdsl)$ n2:=atsoc(acmemb(list(y . 'one,y1 . 'half),lgrids),agdsl)$ nn1:=atsoc(acmemb(list(y . 'half,y1 . 'half),lgrids),agdsl)$ nn2:=atsoc(acmemb(list(y . 'half,y1 . 'one),lgrids),agdsl)$ if n1 and n2 and nn1 and nn2 then t else go to er2$ n1:=cadr n1$ n2:=cadr n2$ nn1:=cadr nn1$ nn2:=cadr nn2$ l21:add2sint(nv,nv1,x,x1,n1,n2,nn1,nn2)$ go to r$ l3:y2:=caaddr gds$ x2:=atsoc(y2,asl)$ if null x2 then go to er1$ x2:=tdifpair cdr x2$ if !*twogrid and flagp(x2,'twogrid) then go to l20$ nv2:=ngetvar(x2,'nrvar)$ n1:=atsoc(acmemb(list(y . 'one,y1 . 'one,y2 . 'one),lgrids),agdsl)$ n2:=atsoc(acmemb(list(y . 'half,y1 . 'half,y2 . 'half),lgrids),agdsl)$ nn1:=atsoc(acmemb(list(y . 'one,y1 . 'one,y2 . 'half),lgrids),agdsl)$ nn2:=atsoc(acmemb(list(y . 'half,y1 . 'half,y2 . 'one),lgrids),agdsl)$ n3:=atsoc(acmemb(list(y . 'one,y1 . 'half,y2 . 'one),lgrids),agdsl)$ n4:=atsoc(acmemb(list(y . 'half,y1 . 'one,y2 . 'half),lgrids),agdsl)$ n5:=atsoc(acmemb(list(y . 'one,y1 . 'half,y2 . 'half),lgrids),agdsl)$ n6:=atsoc(acmemb(list(y . 'half,y1 . 'one,y2 . 'one),lgrids),agdsl)$ if n1 and n2 and nn1 and nn2 and n3 and n4 and n5 and n6 then t else go to er2$ n1:=cadr n1$ n2:= cadr n2$ nn1:=cadr nn1$ nn2:=cadr nn2$ n3:=cadr n3$ n4:=cadr n4$ n5:=cadr n5$ n6:=cadr n6$ if n1=nn1 and n2=nn2 and n3=n5 and n4=n6 then <<n2:=n3$ nn1:=nn2$ nn2:=n4$ go to l21 >> else if n1=n3 and n2=n4 and nn1=n5 and nn2=n6 then <<n2:=nn1$ nn1:=n4$ x1:=x2$ nv1:=nv2$ go to l21 >> else if n1=n6 and n2=n5 and nn1=n4 and nn2=n3 then <<n2:=nn2$ nn1:=n5: nn2:=n4$ x:=x2$ nv:=nv2$ go to l21 >>$ add3sint(nv,nv1,nv2,x,x1,x2,n1,n2,n3,n4,n5,n6,nn1,nn2)$ r:return (nil . 1)$ er:msgpri(nil,l," Free vars not yet implemented ",nil,'hold)$ go to r$ er1:msgpri(" Failed matching of variables in ", u,list(asl,agdsl),nil,'hold)$ go to r$ er2:msgpri(" All grids not given for term ",u,list(asl,agdsl), nil,'hold)$ go to r end$ procedure add2sint(nv,nv1,x,x1,n1,n2,nn1,nn2)$ begin % Enhansment for symmetries, when only one variable influence if n1=n2 and nn1=nn2 then <<f2plus(nvar,nv,0,6*n1)$ f2plus(nvar,nv,1,6*nn1)$ return >> else if n1=nn2 and n2=nn1 then <<f2plus(nvar,nv1,0,6*n1)$ f2plus(nvar,nv1,1,6*n2)$ return >>$ n1:=3*n1$ n2:=3*n2$ nn1:=3*nn1$ nn2:=3*nn2$ x:=list(x,x1)$ f2plus(nvar,nv,0,list(x,n1,n2))$ f2plus(nvar,nv,1,list(x,nn1,nn2))$ x:=reverse x$ f2plus(nvar,nv1,0,list(x,n1,nn2))$ f2plus(nvar,nv1,1,list(x,nn1,n2))$ return end$ procedure add3sint(nv,nv1,nv2,x,x1,x2,n1,n2,n3,n4,n5,n6,nn1,nn2)$ begin n1:=2*n1$ n2:=2*n2$ nn1:=2*nn1$ nn2:=2*nn2$ n3:=2*n3$ n4:=2*n4$ n5:=2*n5$ n6:=2*n6$ x:=list(x,x1,x2)$ f2plus(nvar,nv,0,list(x,n1,nn1,n3,n5))$ f2plus(nvar,nv,1,list(x,n2,nn2,n4,n6))$ f2plus(nvar,nv1,0,list(x,n1,nn1,n4,n6))$ f2plus(nvar,nv1,1,list(x,n2,nn2,n3,n5))$ f2plus(nvar,nv2,0,list(x,n1,nn2,n3,n6))$ f2plus(nvar,nv2,1,list(x,n2,nn1,n4,n5))$ return end$ procedure simpinterpol u$ begin scalar asl,agdsl,gds,x,y,xx,a$ u:=prepsq simp car u$ if eqcar(u,'diff) and not coord memq cddr u then % !!!! Be aware !!!! could not work for mixed derivatives return <<asl:=!*exp$ !*exp:=t$ alglist!*:=nil . nil$ % added 13/06/91 u:= simp formlnr('diff . (prepsq simp formlnr list('intt,cadr u,coord) . cddr u))$ !*exp:=asl$ u>>$ asl:=finddifm u$ if null asl then return (nil . 1) else if eqcar(asl,'special) then return cdr asl$ agdsl:=cadr asl$ % Actual list from DIFML!*, contains definition % of grid, penalty and diff. scheme asl:=car asl$ % Assoc. list of assignments of variables X,U,V,W % to actual variables if not gvar memq grids!* then go to erg$ asl:=append(asl,get('steps,gvar))$ % Adding DIM1, DIP1 ... to assoc. % list if null caar agdsl then return simp sublap(asl,caddar agdsl)$ % For a:=caar agdsl$ % DIFMATCH without def. grids b:if null a then go to c$ y:=caar a$ x:=atsoc(y,asl)$ if null x then go to er1$ % GDS is assoc. list of actual xx:=cdr x$ % assignments of grids to x:=getgrid xx$ % variables U, V if gvar eq 'half then x:=gnot x$ if !*twogrid and twogridp xx then t else gds:=(y . x) . gds$ a:=cdr a$ go to b$ c:if null gds then go to a$ % For given functions which can be on y:=get('difm!*,'grids)$ % both grids x:=acmemb(gds,y)$ % Unique GDS if null x then go to er1$ gds:=x$ a:x:=assoc(gds,agdsl)$ if null x then go to erg$ x:=caddr x$ % Actual difference scheme return simp sublap(asl,x)$ er1:msgpri(" Failed matching of ",u,list(asl,agdsl,gds),nil, 'hold)$ return (nil . 1)$ erg:msgpri(" Bad grids ",u,list(asl,agdsl,gds),nil,'hold)$ return (nil . 1) end$ procedure twogridp u$ % Checks if prefix form U can be on both grids begin scalar x$ return if atom u then if flagp(u,'twogrid) then if !*twogrid and u memq given!* and getel list(get(u,'grid),ncor) then nil else t else nil else if flagp(car u,'twogrid) then if !*twogrid and car u memq given!* and getel list(get(car u,'grid),ncor) then nil else t else if car u memq '(diff plus difference) then twogridp cadr u else if car u eq 'times then twogridpti cdr u else nil end$ procedure twogridpti u$ begin scalar x$ a:x:=twogridp car u$ if x then return x$ u:=cdr u$ if u then go to a$ return nil end$ procedure getgrid u$ begin scalar x$ return if atom u then if x:=get(u,'xgrid) then x else if !*twogrid and u memq given!* and (x:=getel list(get(u,'grid),ncor)) then car x else nil else if (x:=get(car u,'xgrid)) then x else if !*twogrid and car u memq given!* and (x:=getel list(get(car u,'grid),ncor)) then car x else if car u eq 'diff then if atom cadr u then getgrid cadr u %else if caadr u eq 'times then % probably can % if (x:=getgrid cadadr u) then x % be deleted % else getgrid caddr cadr u % !!!!!"!!! else getgrid cadr u else if car u memq '(plus difference) then getgrid cadr u else if car u eq 'times then getgti cdr u else nil end$ procedure getgti u$ begin scalar x$ a:x:=getgrid car u$ if x then return x$ u:=cdr u$ if u then go to a$ return nil end$ procedure sublap(u,v)$ % U is assoc. list, V is pattern diff. scheme % Performs substitution of assod. list into the diff. scheme begin scalar x$ return if null u or null v then v else if atom v then if numberp v then v else if x:=atsoc(v,u) then cdr x else v else if flagp(car v,'app) then sublap1(u,v) else (sublap(u,car v) . sublap(u,cdr v)) end$ flag('(u f v w x g),'app)$ procedure sublap1(u,v)$ begin scalar x,y$ x:=atsoc(car v,u)$ if null x then return msgpri(" Substitution for ",v," not find", nil,'hold)$ x:=cdr x$ y:=mapcar(cdr v,function(lambda a$irev sublap(u,a)))$ return if eqcar(x,'diff) then ('diff . (subappend(cadr x,y) . cddr x)) else subappend(x,y) end$ procedure subappend(x,y)$ if atom x then if x memq iobjs!* and depends(x,coord) then (x . y) else x else if car x memq iobjs!* and depends(car x,coord) then append(x,y) else (car x . mapcar(cdr x,function(lambda a$ subappend(a,y))))$ procedure irev u$ begin u:=simp u$ return if cdaaar u=1 and cdaar u=cdr u and fixp cdar u then if cdr u=1 then if cdar u<0 then list('difference, caaaar u, -cdar u) else list('plus, caaaar u, cdar u) else if cdar u<0 then list('difference, caaaar u, list('quotient, -cdar u, cdr u)) else list('plus, caaaar u, list('quotient, cdar u, cdr u)) else prepsq u end$ unfluid '(coord unvars sunvars interpolp novars ncor nvar intp icor gvar hi hip1 hip2 him1 him2 lhs rhs lsun lun xxgrid resar)$ procedure gentempst$ list('gentemp,xread t)$ put('gentemp,'stat,'gentempst)$ put('gentemp,'formfn,'formgentran)$ put('outtemp,'stat,'endstat)$ flag('(outtemp),'eval)$ algebraic$ endmodule$ end$