spde.red 36 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044
  1. module spde; % Determine Lie symmetries of partial differential eqns.
  2. % Author: Fritz Schwarz.
  3. %*******************************************************************$
  4. % $
  5. % This is the REDUCE package SPDE for determining $
  6. % Lie symmetries of partial differential equations $
  7. % Version of November 1986 $
  8. % $
  9. % $
  10. % Fritz Schwarz $
  11. % GMD Institut F1 $
  12. % Postfach 1240 $
  13. % 5205 St. Augustin $
  14. % West Germany $
  15. % $
  16. % Tel. 02241-142782 $
  17. % EARN Id. DBNGMD21.GF1002 $
  18. %*******************************************************************$
  19. create!-package('(spde),'(contrib spde));
  20. algebraic operator x,u,xi,eta,c,xi!*,eta!*$
  21. algebraic operator deq,dx,du,gl,gen,sder,rule$
  22. fluid '(depl!*);
  23. global'(pclass mm nn num!-cgen num!-dgen)$
  24. share pclass,mm,nn$
  25. lisp(pclass:=mm:=nn:=num!-cgen:=num!-dgen:=0)$
  26. lisp(operator simpsys,result,prsys,prsys!*)$
  27. fluid '(!*list kord!*)$
  28. fluid'(uhf dfsub csub czero rdep !*rational)$
  29. fluid'(list!-m list!-deq list!-pq)$
  30. %symbolic procedure prload$
  31. % begin
  32. % if not getd 'solve1 then load solve1,solvetab,quartic;
  33. % if not getd 'depend1 then load depend;
  34. % if not getd 'ratfunpri then load ratprin;
  35. % end$
  36. symbolic procedure prload; nil;
  37. %*******************************************************************$
  38. % Auxiliary RLISP procedures $
  39. %*******************************************************************$
  40. remflag('(ordp ordpa),'lose); % We must use these definitions.
  41. symbolic procedure ordp(u,v)$
  42. %Modified ordering function which orders kernels with CAR parts;
  43. %DF, ETA, XI and C ahead of anything else;
  44. if null u then null v else if null v then t else
  45. if eq(u,'df) or eq(u,'eta) and not eq(v,'df)
  46. or eq(u,'xi) and not(eq(v,'df) or eq(v,'eta))
  47. or eq(u,'c) and not(eq(v,'df) or eq(v,'eta) or eq(v,'xi)) then t else
  48. if eq(u,'eta) and eq(v,'df)
  49. or eq(u,'xi) and (eq(v,'df) or eq(v,'eta))
  50. or eq(u,'c) and (eq(v,'df) or eq(v,'eta) or eq(v,'xi))
  51. or eq(v,'df) or eq(v,'eta) or eq(v,'xi) or eq(v,'c) then nil else
  52. if atom u then if atom v then
  53. if numberp u then numberp v and not(u<v) else
  54. if numberp v then t else orderp(u,v) else nil else
  55. if atom v then t else
  56. if car u=car v then ordp(cdr u,cdr v) else ordp(car u,car v)$
  57. symbolic procedure ordpa(u,v); ordp(u,v);
  58. symbolic procedure makeset u$
  59. if not u then nil else
  60. if member(car u,cdr u) then makeset cdr u else
  61. car u . makeset cdr u$
  62. symbolic procedure lastmem u$
  63. if cdr u then lastmem cdr u else car u$
  64. symbolic procedure xmember(u,v)$ reverse member(u,reverse v)$
  65. symbolic procedure sacar(a,u)$
  66. if atom u then nil else
  67. if eq(a,car u) and cdr u then list u else
  68. append(sacar(a,car u),sacar(a,cdr u))$
  69. symbolic procedure scar(a,u)$
  70. if atom u then nil else if a=car u then u else
  71. scar(a,car u) or scar(a,cdr u)$
  72. symbolic procedure inter(u,v);
  73. if not u then nil else
  74. if member(car u,v) then
  75. (car u) . inter(cdr u,v) else inter(cdr u,v)$
  76. symbolic procedure compl(u,v)$
  77. if not u then nil else if member(car u,v) then
  78. compl(cdr u,v) else car u . compl(cdr u,v)$
  79. symbolic procedure vlist u$
  80. %U is list of items, returns U with all integers omitted;
  81. if not u then nil else
  82. if numberp car u then vlist cdr u else (car u) . vlist cdr u$
  83. symbolic procedure delnil u$
  84. %U is list, returns U with all occurences of nil deleted;
  85. if not u then nil else
  86. if car u then (car u) . delnil cdr u else delnil cdr u$
  87. symbolic procedure prlist u$
  88. %U is list of items, returns list of all pairs in U;
  89. if not u then nil else
  90. if pairp car u then (car u) . prlist cdr u else prlist cdr u$
  91. symbolic procedure appends(u,v,w)$ append(u,append(v,w))$
  92. symbolic procedure propa(fn,u)$
  93. %FN is predicate of a single argument, U a list;
  94. %Returns T if predicate is true for all elements of U;
  95. begin scalar ind;
  96. ind:=t;
  97. while ind and u do <<ind:=apply1(fn,car u); u:=cdr u>>;
  98. return ind;
  99. end$
  100. symbolic procedure sortx(fn,u)$
  101. begin scalar v,w;
  102. while u do<<v:=maxmem(fn,u); u:=delete(v,u); w:=v . w>>;
  103. return w;
  104. end$
  105. symbolic procedure maxmem(fn,u)$
  106. %FN is function of a single argument, U a list;
  107. %Returns element of U for which FN is maximal;
  108. begin scalar v;
  109. v:=car u;
  110. foreach x in cdr u do
  111. if greaterp(apply1(fn,x),apply1(fn,v)) then v:=x;
  112. return v;
  113. end$
  114. symbolic procedure maxl u$
  115. %U is list of integers, returns largest element of U;
  116. if not u then -10000 else max(car u,maxl cdr u)$
  117. symbolic procedure suml u$
  118. %U is list of integers, returns sum of all elements;
  119. if not u then 0 else plus2(car u,suml cdr u)$
  120. symbolic procedure spde!-subsetp(u,v)$
  121. %U and V are list representing sets;
  122. %Returns T if set U is subset of V;
  123. if not u then t else
  124. member(car u,v) and spde!-subsetp(cdr u,v)$
  125. symbolic procedure product!-set2(u,v)$
  126. %U and V are lists representing sets, returns list representing;
  127. %product set of sets represented by U and V;
  128. begin scalar w;
  129. foreach x in u do foreach y in v do w:=list(x,y) . w;
  130. return w;
  131. end$
  132. symbolic procedure leqgrt(l,i,j)$
  133. i leq j and eqn(l,i) or i geq add1 j$
  134. symbolic procedure fidep u$
  135. assoc(u,depl!*) and cdr assoc(u,depl!*)$
  136. symbolic procedure mkdep u$
  137. foreach x in cdr u do depend1(car u,x,t)$
  138. symbolic procedure rmdep u$
  139. <<rmsubs(); foreach x in cdr u do depend1(car u,x,nil)>>$
  140. symbolic procedure blanks l;
  141. begin scalar u;
  142. u := '(!");
  143. for k:=1:l do u:='! . u;
  144. return compress('!" . u)
  145. end$
  146. symbolic procedure terpri2$ <<terpri(); terpri()>>$
  147. %*******************************************************************$
  148. % Auxiliary procedures for manipulating standard forms $
  149. %*******************************************************************$
  150. symbolic procedure lcf u$ not domainp u and lc u$
  151. symbolic procedure minus!-f u$
  152. %U is s.f., returns T if lnc U is negative;
  153. minusf numr simp reval u$
  154. lisp operator minus!-f$
  155. symbolic procedure lengthn u$
  156. if not u then 0 else
  157. if numberp car u then plus(sub1 car u,lengthn cdr u) else
  158. plus(1,lengthn cdr u)$
  159. symbolic procedure lengthf u$
  160. %U is prefix s.f., returns printlength for U;
  161. if not u then 0 else
  162. if atom u then flatsizec u else
  163. if eqcar(u,'plus)
  164. then plus(times(3,sub1 length cdr u),lengthf cdr u) else
  165. if eqcar(u,'times) or eqcar(u,'minus)
  166. then plus(sub1 length cdr u,lengthf cdr u) else
  167. if eqcar(u,'quotient) then
  168. if !*rational then add1 add1 max(flatsizec cadr u,flatsizec caddr u)
  169. else add1 plus(flatsizec cadr u,flatsizec caddr u) else
  170. if eqcar(u,'expt) then add1 flatsizec cadr u else
  171. if eqcar(u,'dx) or eqcar(u,'du) then plus(flatsizec cadr u,4) else
  172. if eqcar(u,'xi) or eqcar(u,'eta) or eqcar(u,'c)
  173. or eqcar(u,'x) or eqcar(u,'u) then times(2,length u) else
  174. if eqcar(u,'df) then plus(4,lengthf cadr u,lengthf cddr u) else
  175. plus(lengthf car u,lengthf cdr u)$
  176. lisp operator lengthf$
  177. symbolic procedure diford u$ lengthn cddr u$
  178. symbolic procedure adiff(u,v)$
  179. %U is kernel with CAR part DF, V is kernel;
  180. %Returns U integrated with respect to V;
  181. if not member(v,u) then u else
  182. if length u=3 and member(v,u) then cadr u else
  183. if not cdr member(v,u)
  184. or not numberp cadr member(v,u) then delete(v,u) else
  185. if cadr member(v,u)=2 then
  186. append(xmember(v,u),cddr member(v,u)) else
  187. append(xmember(v,u),(sub1 cadr member(v,u)) . cddr member(v,u))$
  188. symbolic procedure sub!-int!-df u$
  189. %U is kernel with CAR part INT, returns integrated kernel if CADR;
  190. %part of U is DF and integration variable occurs as argument of DF;
  191. if eqcar(cadr u,'df) and member(lastmem u,cadr u) then
  192. adiff(cadr u,lastmem u) else u$
  193. symbolic procedure subintf u$
  194. %U is s.f., performs all integrations which may be done;
  195. %by cancellation of corresponding differentiation;
  196. begin
  197. foreach x in makeset sacar('int,u) do
  198. u:=subst(sub!-int!-df x,x,u);
  199. return numr simp prepf u;
  200. end$
  201. symbolic procedure monop u$
  202. %Returns T if u is monomial;
  203. domainp u or not red u and monop lc u$
  204. symbolic procedure solvef(u,v)$ car solve0(prepf u,v)$
  205. symbolic procedure comfacn u$ lnc ckrn u$
  206. symbolic procedure remfacn u$ quotf(u,lnc ckrn u)$
  207. %*******************************************************************$
  208. % Procedures for manipulating l.d.f.'s, U is always l.d.f. $
  209. % in this section $
  210. %*******************************************************************$
  211. symbolic procedure ldf!-mvar u$
  212. %Returns function argument of mvar U;
  213. (if eqcar(x,'df) then cadr x else x) where x=mvar u;
  214. symbolic procedure ldf!-fvar u$
  215. %Returns all function arguments occuring in U;
  216. makeset foreach x in u collect ldt!-tvar x$
  217. symbolic procedure ldf!-fvar!-part(u,v)$
  218. %V is function xi(i), eta(alpha) or c(k), returns l.d.f. of those;
  219. %terms in U with ldt-tvar x equal to V, overall factors not removed;
  220. begin scalar w;
  221. foreach x in u do if eq(ldt!-tvar x,v) then w:=x . w;
  222. return reverse w;
  223. end$
  224. symbolic procedure ldf!-dep!-var u$
  225. %Returns all variables x(i) or u(alpha) which occur as;
  226. %arguments of XI, ETA or C;
  227. begin scalar v;
  228. foreach x in u do if assoc(ldt!-tvar x,depl!*) then
  229. v:=append(cdr assoc(ldt!-tvar x,depl!*),v);
  230. return makeset v;
  231. end$
  232. symbolic procedure ldf!-pow!-var u$
  233. %Returns all variables x(i) or u(alpha) which occur as powers;
  234. begin scalar v,z;
  235. foreach x in u do v:=append(v,kernels tc x);
  236. foreach y in prlist makeset v do
  237. if eqcar(y,'x) or eqcar(y,'u) then z:=y . z;
  238. return makeset z;
  239. end$
  240. symbolic procedure ldf!-deg(u,v)$
  241. %V is kernel x(i) or u(alpha), returns degree of U in V;
  242. maxl foreach x in u collect degreef(tc x,v)$
  243. symbolic procedure ldf!-spf!-var u$
  244. %Returns all variables x(i) or u(alpha) which occur as;
  245. %arguments of any other kernel than xi, eta or c;
  246. begin scalar v,z;
  247. foreach x in u do v:=append(v,kernels tc x);
  248. foreach y in prlist makeset v do
  249. if not eqcar(y,'x) and not eqcar(y,'u) then
  250. z:=appends(sacar('x,cdr y),sacar('u,cdr y),z);
  251. return makeset z;
  252. end$
  253. symbolic procedure ldf!-all!-var u$
  254. %Returns all variables x(i) or u(alpha) which occur in U;
  255. makeset appends(ldf!-dep!-var u,ldf!-pow!-var u,ldf!-spf!-var u)$
  256. symbolic procedure ldf!-sep!-var u$
  257. %Returns all variables w.r.t. which U may be separated;
  258. compl(compl(ldf!-pow!-var u,ldf!-dep!-var u),ldf!-spf!-var u)$
  259. symbolic procedure ldf!-int!-var u$
  260. %Returns all variables w.r.t. which U may be integrated;
  261. if eqcar(mvar u,'df) then
  262. begin scalar v;
  263. v:=ldf!-all!-var u;
  264. while v and u do
  265. <<v:=compl(v,compl(ldt!-dep car u,ldt!-dfvar car u)); u:=cdr u>>;
  266. return v;
  267. end$
  268. symbolic procedure ldf!-int u$
  269. %U is l.d.f, returns U with all possible integrations performed;
  270. %or unchanged if integration is not possible;
  271. begin scalar v,w,z,test; integer nfun;
  272. a:
  273. test:=nil;
  274. w:=ldf!-int!-var u;
  275. nfun:=find!-nfun();
  276. foreach x in w do
  277. if not smember('int,z:=caadr algebraic int(lisp prepf u,x))
  278. or not smember('int,z:=subintf z) then
  279. <<v:=!*a2k list('c,nfun:=nfun+1); test:=t;
  280. mkdep(v . delete(x,ldf!-all!-var u));
  281. u:=addf(z,!*k2f v)>>;
  282. if test then go to a;
  283. return u;
  284. end$
  285. symbolic procedure ldf!-df!-diff u$
  286. %Returns list of all df-kernels which may be obtained;
  287. %from U by differentiation or nil;
  288. begin scalar dfvar,dfsub,v,w,z0,z; integer n0,nmax;
  289. v:=compl(ldf!-dep!-var u,ldf!-spf!-var u);
  290. if not v then return;
  291. w:=foreach x in v collect list(x,add1 ldf!-deg(u,x));
  292. nmax:=maxl foreach x in w collect cadr x;
  293. while (n0:=n0+1) leq nmax and not(z0:=nil) do
  294. <<foreach x in w do if cadr x geq n0 then z0:=(car x) . z0;
  295. z:=z0 . z>>;
  296. z:=reverse z;
  297. dfvar:=foreach x in car z collect list x;
  298. foreach x in cdr z do dfvar:=
  299. append(dfvar,foreach y in dfvar collect car product!-set2(x,y));
  300. foreach x in dfvar do
  301. begin scalar p,q;
  302. p:=x; q:=u;
  303. while p and q and red q do
  304. <<q:=ldf!-simp numr difff(q,car p); p:=cdr p>>;
  305. if pairp q and not red q and eqcar(mvar q,'df) then
  306. dfsub:=(mvar q) . dfsub;
  307. end;
  308. return makeset dfsub;
  309. end$
  310. symbolic procedure ldf!-sub!-var u$
  311. %Returns function w.r.t. which U may be resolved;
  312. begin scalar v,w,z;
  313. w:=ldf!-all!-var u;
  314. foreach x in u do if not v and not eqcar(z:=tvar x,'df)
  315. and monop tc x and spde!-subsetp(w,ldt!-dep x)
  316. and not smember(z,delete(x,u)) then v:=z;
  317. return v;
  318. end$
  319. symbolic procedure ldf!-simp u$
  320. %Returns l.d.f. form of U;
  321. if not u then nil else
  322. if not red u then numr simp prepf !*k2f mvar u else
  323. begin scalar v;
  324. v:=numr simp prepf u;
  325. if not domainp v then v := quotf(v,cdr comfac v);
  326. return absf v
  327. end$
  328. symbolic procedure ldf!-sep u$
  329. %Returns list of l.d.f. into which U has been separated;
  330. begin scalar v; integer k;
  331. if not(v:=ldf!-sep!-var u) then return list u;
  332. foreach x in v do u:=subst(list('ux,1,k:=k+1),x,u);
  333. return foreach x in coeff!-all(u,'ux) collect
  334. ldf!-simp numr simp prepf x;
  335. end$
  336. symbolic procedure ldf!-subf0 u$
  337. %Returns U with CZERO substituted;
  338. ldf!-simp delnil foreach x in u collect ldt!-subt0 x$
  339. %*******************************************************************$
  340. % Procedures for manipulating l.d.t.'s, U is always l.d.t. $
  341. % in this section $
  342. %*******************************************************************$
  343. symbolic procedure ldt!-tvar u$
  344. %U is l.d.t., returns function argument of tvar U;
  345. (if eqcar(x,'df) then cadr x else x) where x=tvar u$
  346. symbolic procedure ldt!-dfvar u$
  347. %U is l.d.t., returns variables w.r.t. which tvar u is derived or nil;
  348. (if eqcar(x,'df) then vlist cddr x else nil) where x=tvar u$
  349. symbolic procedure ldt!-dep u$
  350. %U is l.d.t., returns list of variables x or y which occur as;
  351. %arguments LDT-tvar u;
  352. (if x then cdr x else nil) where x=assoc(ldt!-tvar u,depl!*)$
  353. symbolic procedure ldt!-subt0 u$
  354. %U is l.d.t., returns U if LDT-tvar u is not on czero;
  355. if not member(ldt!-tvar u,czero) then u else nil$
  356. %*******************************************************************$
  357. % Procedures for constructing the determining system $
  358. %*******************************************************************$
  359. symbolic procedure cresys u$
  360. begin scalar r,v,w,lgl,lsub,depl!*!*,list!-sder;
  361. remprop('df,'kvalue); remprop('df,'klist);
  362. remprop('c,'kvalue); remprop('c,'klist);
  363. prload();
  364. rmsubs();
  365. depl!*:=nil;
  366. if car u then
  367. list!-deq:=foreach x in u collect assoc(x,get(car x,'kvalue)) else
  368. list!-deq:=get('deq,'kvalue);
  369. if eqn(length list!-deq,1) then
  370. begin scalar p;
  371. p:=maxmem(function length,makeset sacar('u,list!-deq));
  372. p:=mk!*sq !*k2q p;
  373. list!-sder:=list list(list('sder, cadaar list!-deq),p);
  374. end else if car u then
  375. list!-sder:=foreach x in list!-deq collect
  376. assoc(list('sder,cadar x),get('sder,'kvalue)) else
  377. list!-sder:=get('sder,'kvalue);
  378. if not list!-deq then rerror(spde,1,
  379. "Differential equations not defined");
  380. if not list!-sder then rerror(spde,2,
  381. "Substitutions for derivatives not defined");
  382. mm:=find!-m list!-deq; nn:=find!-n list!-deq;
  383. list!-m:=
  384. makeset foreach x in sacar('u,list!-deq) collect cadr x;
  385. for k:=1:nn do<<w:=!*a2k list('xi,k) . w; v:=!*a2k list('x,k) . v>>;
  386. for k:=1:mm do if member(k,list!-m) then
  387. <<w:=!*a2k list('eta,k) . w; v:=!*a2k list('u,k) . v>>;
  388. for k:=1:nn do r:=(!*a2k list('dx,k)) . r;
  389. for k:=1:mm do r:=(!*a2k list('du,k)) . r;
  390. for k:=1:mm do depl!*!*:=(!*a2k list('eta,k) . v) . depl!*!*;
  391. for k:=1:nn do depl!*!*:=(!*a2k list('xi,k) . v) . depl!*!*;
  392. depl!*:=depl!*!*;
  393. kord!*:=reverse r;
  394. foreach x in list!-sder do
  395. lsub:=((mvar caadr cadr x) . prepsq caar solvef(caadr cadr assoc
  396. (list('deq,cadar x),list!-deq),mvar caadr cadr x)) . lsub;
  397. foreach x in list!-deq do
  398. begin scalar s,z,lx,lu;
  399. z:=caadr cadr x;
  400. lx:=makeset sacar('x,z);
  401. lu:=makeset sacar('u,z);
  402. foreach y in lx do s:=addf(s,
  403. multf(!*k2f !*a2k list('xi,cadr y),numr simp prepsq difff(z,y)));
  404. foreach y in lu do if length y=2 then
  405. s:=addf(s,multf
  406. (!*k2f !*a2k list('eta,cadr y),numr simp prepsq difff(z,y))) else
  407. s:=addf(s, multf(numr zeta!* cdr y,numr simp prepsq difff(z,y)));
  408. s:=numr subf(s,lsub);
  409. s:=numr subf(s,lsub);
  410. lgl:=append(coeff!-all(s,'u),lgl);
  411. end;
  412. uhf:=list(makeset lgl,foreach x in reverse w collect !*k2q x);
  413. end$
  414. lisp rlistat'(cresys)$
  415. symbolic procedure totder(u,i)$
  416. begin scalar z,v,w;
  417. v:=car difff(u,!*a2k list('x,i));
  418. z:=makeset sacar('u,u);
  419. for k:=1:mm do if member(k,list!-m) then
  420. z:=(!*a2k list('u,k)) . z;
  421. foreach x in makeset z do w:=addf(w,
  422. multf(!*k2f !*a2k append(x,list i),car difff(u,x)));
  423. return numr simp prepf addf(v,w);
  424. end$
  425. symbolic procedure zeta!* u$
  426. if not get('deq,'kvalue) and (eqn(mm,0) or eqn(nn,0)) then
  427. rerror(spde,3,"Number of variables not defined") else
  428. if length u geq 3 then
  429. begin scalar v,w;
  430. prload();
  431. if eqn(nn,0) then nn:=find!-n list!-deq;
  432. v:=totder(numr zeta!* reverse cdr reverse u,car reverse u);
  433. for s:=1:nn do w:=addf(w,
  434. multf(!*k2f !*a2k('u . append(reverse cdr reverse u,list s)),
  435. totder(!*k2f !*a2k list('xi,s),car reverse u)));
  436. return simp prepsq(addf(v,negf w) ./ 1);
  437. end else
  438. begin scalar v,w;
  439. prload();
  440. if eqn(nn,0) then
  441. <<nn :=find!-n list!-deq; mm:=find!-m list!-deq>> else
  442. begin scalar p,z;
  443. for k:=1:mm do z:=cons(k,z);
  444. for k:=1:nn do p:=(!*a2k list('x,k)) . p;
  445. for k:=1:mm do p:=(!*a2k list('u,k)) . p;
  446. for k:=1:nn do mkdep((!*a2k list('xi,k)) . p);
  447. for k:=1:mm do mkdep((!*a2k list('eta,k)) . p);
  448. list!-m:=z;
  449. end;
  450. v:=totder(!*k2f !*a2k list('eta,car u),cadr u);
  451. for s:=1:nn do w:=addf(w,
  452. multf(!*k2f !*a2k list('u,car u,s),
  453. totder(!*k2f !*a2k list('xi,s),car reverse u)));
  454. return simp prepsq(addf(v,negf w) ./ 1);
  455. end$
  456. symbolic procedure simpu u$
  457. !*p2q mksp(('u . (car u . reverse ordn cdr u)),1)$
  458. put('u,'simpfn,'simpu)$
  459. put('zeta,'simpfn,'zeta!*)$
  460. symbolic procedure coeff!-all(u,v)$
  461. begin scalar z;
  462. list!-pq:=nil;
  463. splitrec(u,v,1,nil);
  464. foreach x in list!-pq do
  465. z:=(ldf!-simp numr simp prepf cdr x) . z;
  466. return makeset z;
  467. end$
  468. symbolic procedure splitrec(u,v,p,q)$
  469. if domainp u then
  470. begin scalar y;
  471. p:=multf(u,p);
  472. if y:=assoc(q,list!-pq) then
  473. rplacd(y,addf(cdr y,p)) else list!-pq:=(q . p) . list!-pq;
  474. end else
  475. begin
  476. if eqcar(mvar u,v) and length mvar u greaterp 2
  477. then splitrec(lc u,v,p,(lpow u) . q)
  478. else splitrec(lc u,v,!*t2f(lpow(u) .* p),q);
  479. if red u then splitrec(red u,v,p,q);
  480. end$
  481. symbolic procedure find!-m u$
  482. maxl makeset foreach x in sacar('u,u) collect cadr x$
  483. symbolic procedure find!-n u$
  484. begin scalar vx,vu,wx,wu;
  485. vx:=makeset sacar('x,u);
  486. vu:=makeset sacar('u,u);
  487. foreach x in vx do wx:=(cadr x) . wx;
  488. foreach x in vu do if length x geq 3 then
  489. wu:=append(cddr x,wu);
  490. return max(maxl wx,maxl wu);
  491. end$
  492. %*******************************************************************$
  493. % Procedures for solving the determining system $
  494. %*******************************************************************$
  495. symbolic procedure rule0$
  496. %Searches for equations of the form C(I)=0 and stores them on CZERO;
  497. if uhf then foreach x in car uhf do
  498. if not red x and not eqcar(mvar x,'df)
  499. then czero:=(mvar x) . czero$
  500. symbolic procedure rule1$
  501. %Searches for equations of the form DF(function,variable)=0;
  502. %and stores it on the list RDEP;
  503. if uhf and car uhf then
  504. begin scalar dfsub;
  505. foreach x in car uhf do
  506. if not red x and eqcar(mvar x,'df) and eqn(diford mvar x,1)
  507. then rdep:=(mvar x) . rdep;
  508. if rdep then return t;
  509. end$
  510. symbolic procedure rule1!-diff$
  511. %Searches for equations of the form DF(function,variable)=0;
  512. %which may be obtained by a single differentiation and stores it on;
  513. %the list RDEP;
  514. if uhf and car uhf then
  515. begin scalar u,v,z;
  516. foreach x in car uhf do if(z:=ldf!-df!-diff x) then
  517. u:=append(z,u);
  518. foreach x in u do if eqn(diford x,1) then v:=x . v;
  519. rdep:=makeset v;
  520. if rdep then return t;
  521. end$
  522. symbolic procedure rulec l$
  523. %Searches for equations of length L which may be solved for a;
  524. %function and stores the corresponding rules on CSUB;
  525. if uhf and car uhf then
  526. begin scalar v;
  527. foreach u in car uhf do if leqgrt(length u,l,4)
  528. and (v:=ldf!-sub!-var u) and not smember(v,csub)
  529. and not inter(foreach x in csub collect car x,ldf!-fvar u)
  530. then csub:=(v . prepsq caar solvef(u,v)) . csub;
  531. if csub then return t;
  532. end$
  533. symbolic procedure ruledf l$
  534. %Searches for equations of the form DF(function,derivative list)=0;
  535. %the derivative beeing of order L and stores the resulting;
  536. %substitution polynomial on CSUB;
  537. if uhf and car uhf then
  538. begin scalar dfsub;
  539. foreach x in car uhf do
  540. if not red x and eqcar(mvar x,'df) and eqn(diford mvar x,l)
  541. and not smember(ldf!-mvar x,dfsub) then dfsub:=(mvar x) . dfsub;
  542. csub:=foreach x in dfsub collect(cadr x) . crepol x;
  543. if csub then return t;
  544. end$
  545. symbolic procedure ruledf!-diff l$
  546. %Searches for all equations of the form;
  547. %DF(function,derivative list)=0 which may be obtained by;
  548. %differentiation, picks out those of order L and stores;
  549. %the corresponding substitution polynomial on CSUB;
  550. if uhf and car uhf then
  551. begin scalar v,dfsub;
  552. foreach u in car uhf do v:=append(v,ldf!-df!-diff u);
  553. if not(v:=makeset v) then return;
  554. foreach x in v do if eqn(diford x,l) then dfsub:=x . dfsub;
  555. if not dfsub then return;
  556. csub:=((cadar dfsub) . crepol car dfsub) . csub;
  557. if csub then return t;
  558. end$
  559. symbolic procedure rule!-int l$
  560. %Searches for an equation of length L which may be solved for a;
  561. %function after beeing integrated and stores the corresponding;
  562. %rule on CSUB;
  563. if uhf and car uhf then
  564. begin scalar v,w;
  565. foreach u in car uhf do if not csub and leqgrt(length u,l,4)
  566. and (v:=ldf!-sub!-var(w:=ldf!-int u))
  567. then csub:=list(v . prepsq caar solvef(w,v));
  568. if csub then return t;
  569. end$
  570. symbolic procedure simpsys0$
  571. %Removes variable which are stored on list CZERO;
  572. begin scalar u,v;
  573. if pclass=2 then<<write"Entering SIMPSYS0"; terpri2()>>;
  574. u:=delnil foreach x in car uhf collect ldf!-subf0 x;
  575. v:=foreach x in cadr uhf collect ldf!-subf0 numr x ./ denr x;
  576. uhf:=list(makeset u,v);
  577. if pclass=1 then
  578. begin
  579. terpri2();
  580. if eqn(length czero,1) then
  581. write"Substitution" else write"Substitutions";
  582. terpri();
  583. foreach x in czero do
  584. algebraic write (lisp aeval x),":=0";
  585. terpri();
  586. end;
  587. if pclass=2 then<<write"CZERO:="; prettyprint czero; terpri()>>;
  588. czero:=nil;
  589. if pclass=2 then<<write"Leaving SIMPSYS0"; terpri2()>>;
  590. end$
  591. symbolic procedure simpsys!-rdep$
  592. %Removes dependencies which are stored on list RDEP;
  593. begin scalar u,v;
  594. if pclass=2 then<<write"Entering SIMPSYS!-RDEP"; terpri2()>>;
  595. foreach x in rdep do rmdep cdr x;
  596. u:=makeset delnil foreach x in car uhf collect ldf!-simp x;
  597. v:=foreach x in cadr uhf collect simp prepsq x;
  598. uhf:=list(u,v);
  599. if pclass=1 then
  600. begin
  601. terpri();
  602. write"Dependencies removed"; terpri2();
  603. foreach x in rdep do
  604. <<maprin cadr x; prin2!*" independent of ";
  605. maprin caddr x; terpri!* t;>>;
  606. terpri();
  607. end;
  608. if pclass=2 then<<write"RDEP:='"; prettyprint rdep; terpri()>>;
  609. if pclass=2 then<<write"Leaving SIMPSYS!-RDEP"; terpri2()>>;
  610. end$
  611. symbolic procedure simpsys!-sep$
  612. %Performs all possible separations;
  613. if uhf and car uhf then
  614. begin scalar u,v,test;
  615. if pclass=2 then<<write"Entering SIMPSYS!-SEP"; terpri2()>>;
  616. foreach x in car uhf do
  617. if eqn(length(v:=ldf!-sep x),1) then u:=x . u else
  618. begin
  619. u:=append(v,u);
  620. if pclass=1 or pclass=2 then
  621. begin scalar z; integer l;
  622. terpri();
  623. l:=length car uhf-length member(x,car uhf)+1;
  624. write"Equation ",l," separated into the terms";
  625. terpri();
  626. if pclass=1 then for k:=1:length v do
  627. begin
  628. z:=prepf nth(v,k);
  629. !*list := lengthf z geq 50;
  630. algebraic write"Term ",k," ",z;
  631. end;
  632. if pclass=2 then foreach y in v do prettyprint y;
  633. end;
  634. test:=t;
  635. end;
  636. !*list := nil;
  637. if test then uhf:=list(reverse makeset u,cadr uhf);
  638. if pclass=2 then<<write"Leaving SIMPSYS!-SEP"; terpri2()>>;
  639. end$
  640. symbolic procedure simpsys!-sub$
  641. %Performs all substitutions which are stored on CSUB;
  642. if uhf and car uhf then
  643. begin scalar u,v;
  644. if pclass=2 then<<write"Entering SIMPSYS!-SUB"; terpri2()>>;
  645. if pclass=1 then prrule csub;
  646. if pclass=2 then<<write"CSUB:='"; prettyprint csub; terpri()>>;
  647. u:=makeset delnil foreach x in car uhf collect
  648. ldf!-simp numr subf(x,csub);
  649. v:=foreach x in cadr uhf collect subsq(x,csub);
  650. uhf:=list(u,v);
  651. csub:=nil;
  652. if pclass=2 then<<write"Leaving SIMPSYS!-SUB"; terpri2()>>;
  653. end$
  654. symbolic procedure simpsys$
  655. if not uhf then
  656. rerror(spde,4,"The determining system is not defined") else
  657. if not car uhf then
  658. rerror(spde,5,"The determining system completely solved") else
  659. begin scalar u,v; integer nfun;
  660. prload();
  661. u:=makeset delnil foreach x in car uhf collect ldf!-simp x;
  662. v:=foreach x in cadr uhf collect simp prepsq x;
  663. uhf:=list(u,v);
  664. mark0:
  665. if pclass=1 then<<prsys!*"Entering main loop">> else
  666. if pclass=2 then prtlist"Entering main loop";
  667. czero:=csub:=rdep:=nil;
  668. simpsys!-sep();
  669. rule0();
  670. if czero then<<simpsys0(); go to mark0>>;
  671. if rule1() or rule1!-diff() then<<simpsys!-rdep(); go to mark0>>;
  672. if ruledf 2 or rulec 2 or rule!-int 2 or ruledf!-diff 2
  673. or ruledf 3 or rulec 3 or rule!-int 3 or ruledf!-diff 3
  674. or ruledf 4 or rulec 4 or rule!-int 4 or ruledf!-diff 4
  675. or ruledf 5 or rulec 5 or rule!-int 5 or ruledf!-diff 5
  676. then <<simpsys!-sub(); go to mark0>>;
  677. if car uhf then
  678. <<write"Determining system is not completely solved";
  679. terpri2(); prsys!*"The remaining equations are";
  680. if not zerop(nfun:=find!-nfun()) then
  681. write"Number of functions is ",nfun>>;
  682. end$
  683. symbolic procedure crepol u$
  684. begin scalar l1,f; integer pow,nfun;
  685. nfun:=find!-nfun();
  686. l1:=cdr assoc(car(u:=cdr u),depl!*);
  687. while (u:=cdr u) do
  688. begin scalar v;
  689. v:=car u;
  690. if length u=1 or not numberp cadr u then pow:=1 else
  691. <<pow:=cadr u; u:=delete(pow,u);>>;
  692. for k:=1:pow do
  693. begin scalar w;
  694. w:=!*a2k list('c,nfun:=nfun+1);
  695. mkdep(w . delete(v,l1));
  696. if k=1 then f:=w . f;
  697. if k=2 then f:=list('times,w,v) . f;
  698. if k geq 3 then
  699. f:=list('times,w,list('expt,v,k-1)) . f;
  700. end;
  701. end;
  702. return append('(plus),f);
  703. end$
  704. %*************************************************************$
  705. % Procedures for analysing the result $
  706. %*************************************************************$
  707. symbolic procedure cpar u$
  708. begin scalar v;
  709. v:=makeset appends(sacar('xi,u),sacar('eta,u),sacar('c,u));
  710. foreach x in v do if not assoc(x,depl!*) then v:=delete(x,v);
  711. return v;
  712. end$
  713. symbolic procedure makeset!-c!-x u$
  714. if not u then nil else
  715. if member!-c!-x(car u,cdr u) then makeset!-c!-x cdr u else
  716. car u . makeset!-c!-x cdr u$
  717. symbolic procedure member!-c!-x(u,v)$
  718. if not v then nil else
  719. if equal!-c!-x(u,car v) then v else member!-c!-x(u,cdr v)$
  720. symbolic procedure equal!-c!-x(u,v)$
  721. begin scalar p,q;
  722. p:=scar('c,u) or scar('xi,u) or scar('eta,u);
  723. q:=scar('c,v) or scar('xi,v) or scar('eta,v);
  724. return equal(subst('cxx,p,u),subst('cxx,q,v));
  725. end$
  726. symbolic procedure numgen$ length get('gen,'kvalue)$
  727. symbolic operator numgen$
  728. symbolic procedure gengen$
  729. begin scalar u,z,cgen,dgen; integer ngen;
  730. remprop('gen,'kvalue); remprop('gen,'klist);
  731. foreach x in cadr uhf do u:=append(ldf!-fvar numr x,u);
  732. foreach x in makeset u do
  733. begin scalar v,w;
  734. w:=nil ./ 1;
  735. if assoc(x,depl!*) then
  736. v:=foreach y in cadr uhf collect
  737. simp prepsq(ldf!-fvar!-part(numr y,x) ./denr y) else
  738. v:=foreach y in cadr uhf collect
  739. simp prepsq((lcf ldf!-fvar!-part(numr y,x)) ./denr y);
  740. for k:=1:nn do if numr nth(v,k) then
  741. w:=addsq(multsq(nth(v,k),!*k2q !*a2k list('dx,k)),w);
  742. for k:=1:mm do if numr nth(v,nn+k) then
  743. w:=addsq(multsq(nth(v,nn+k),!*k2q !*a2k list('du,k)),w);
  744. if assoc(x,depl!*) then
  745. cgen:=(absf remfacn numr simp prepf numr w) . cgen else
  746. dgen:=(absf remfacn numr simp prepf numr w) . dgen;
  747. end;
  748. dgen:=makeset dgen; cgen:=makeset!-c!-x cgen;
  749. num!-dgen:=length dgen; num!-cgen:=length cgen;
  750. for k:=1:nn do if member(z:=!*k2f !*a2k list('dx,k),dgen) then
  751. <<setk(list('gen,ngen:=add1 ngen),prepf z); dgen:=delete(z,dgen)>>;
  752. for k:=1:mm do if member(z:=!*k2f !*a2k list('du,k),dgen) then
  753. <<setk(list('gen,ngen:=add1 ngen),prepf z); dgen:=delete(z,dgen)>>;
  754. dgen:=sortx(function length,dgen);
  755. foreach x in dgen do setk(list('gen,ngen:=add1 ngen),prepf x);
  756. cgen:=sortx(function length,cgen);
  757. foreach x in cgen do setk(list('gen,ngen:=add1 ngen),prepf x);
  758. end$
  759. symbolic operator gengen$
  760. algebraic procedure comm(a,b)$
  761. begin scalar z;
  762. if (lisp length list!-deq)=0 then
  763. <<write"Differential equations not defined"; return nil>>;
  764. z:= (for k:=1:nn sum df(a,dx k)*df(b,x k)-df(b,dx k)*df(a,x k))
  765. +(for k:=1:mm sum df(a,du k)*df(b,u k)-df(b,du k)*df(a,u k))$
  766. return z;
  767. end$
  768. algebraic procedure result$
  769. begin integer l;
  770. if (l:=lisp length list!-deq)=1 then
  771. write"The differential equation" else
  772. write"The differential equations";
  773. for j:=1:l do
  774. begin scalar z; integer i,k;
  775. lisp(z:=car cadadr nth(list!-deq,j));
  776. i:=lisp cadar nth(list!-deq,j);
  777. k:=lisp lengthf prepf z;
  778. symbolic(!*list := k>40);
  779. write"DEQ(",i,"):=",lisp prepf z;
  780. end;
  781. !*list := nil;
  782. if (lisp length car uhf) neq 0 then
  783. prsys!*"The determining system is not completely solved" else
  784. <<lisp gengen(); prgen(); comm!-tab()>>;
  785. end$
  786. %*************************************************************$
  787. % Procedures for displaying the output $
  788. %*************************************************************$
  789. symbolic procedure prsys!* u$
  790. if uhf and car uhf then
  791. <<terpri(); write u; terpri(); prsys(); terpri()>>$
  792. symbolic procedure prsys$
  793. begin scalar v;
  794. terpri();
  795. remprop('gl,'kvalue); remprop('gl,'klist);
  796. for k:=1:length car uhf do
  797. begin scalar z; integer l;
  798. z:=prepf nth(car uhf,k);
  799. l:=lengthf prepf nth(car uhf,k);
  800. !*list := l>50;
  801. algebraic write"GL(",k,"):=",z;
  802. setk(list('gl,k),z);
  803. end;
  804. terpri2();
  805. write"The remaining dependencies";
  806. terpri2();
  807. v:=makeset
  808. appends(sacar('xi,car uhf),sacar('eta,car uhf),sacar('c,car uhf));
  809. foreach x in v do write!-dep x;
  810. !*list := nil;
  811. end$
  812. symbolic procedure prrule u$
  813. begin
  814. terpri2();
  815. if eqn(length u,1) then
  816. write"Substitution" else write"Substitutions";
  817. terpri2();
  818. foreach x in u do
  819. <<maprin car x; prin2!*" = "; maprin cdr x; terpri!* t;>>;
  820. terpri();
  821. foreach x in u do foreach y in sacar('c,cdr x) do write!-dep y;
  822. end$
  823. symbolic procedure prtlist u$
  824. <<write u; terpri2(); write"DEPL!*:='"; prettyprint depl!*;
  825. write"UHF:='"; prettyprint uhf>>$
  826. symbolic procedure write!-df!-sub$
  827. if get('df,'kvalue) then
  828. begin scalar w;
  829. w:=get('df,'kvalue);
  830. remprop('df,'kvalue);
  831. terpri();
  832. if length w=1 then write"Constraint" else write"Constraints";
  833. terpri2();
  834. foreach x in w do
  835. begin scalar u,v;
  836. u:=car x;
  837. v:=cadadr x;
  838. algebraic write lisp u,":=",lisp prepsq v;
  839. terpri2();
  840. end;
  841. put('df,'kvalue,w);
  842. end$
  843. algebraic procedure prgen$
  844. begin scalar lcpar;
  845. for k:=1:nn do <<order dx k; factor dx k>>;
  846. for k:=1:mm do factor du k$
  847. lisp(lcpar:=cpar get('gen,'kvalue));
  848. write"The symmetry generators are";
  849. for k:=1:numgen() do
  850. if (lisp lengthf reval list('gen,k)) leq 60 then
  851. <<symbolic(!*list := nil); write"GEN(",k,"):=",gen k>> else
  852. begin scalar z; integer r,s,nt; operator gen!*;
  853. nt:=lisp length(z:=numr simp reval list('gen,k));
  854. r:=lisp maxl foreach x in z collect abs comfacn list x;
  855. if r=1 then r:=0 else r:=lisp flatsizec r;
  856. for l:=1:nt do gen!* l:=lisp prepf list nth(z,l);
  857. for l:=1:nt do
  858. begin
  859. symbolic(!*list := lengthf prepf tc nth(z,l) geq 56);
  860. s:=lisp abs comfacn list nth(z,l);
  861. if r=0 then s:=0 else
  862. if s=1 then s:=-1 else s:=lisp flatsizec s;
  863. if l=1 then write"GEN(",k,"):=",lisp blanks(r-s+1),gen!* 1 else
  864. if minus!-f gen!* l then
  865. write lisp blanks(r-s+6),gen!* l else
  866. write lisp blanks(r-s+6)," + ",gen!* l;
  867. end;
  868. clear gen!*;
  869. symbolic(!*list := nil);
  870. end;
  871. if (lisp length lcpar) neq 0 then
  872. <<write"The remaining dependencies"; lisp terpri()>>;
  873. for k:=1:(lisp length lcpar) do
  874. <<lisp write!-dep nth(lcpar,k);>>;
  875. if (lisp length lcpar) neq 0 then lisp terpri();
  876. lisp write!-df!-sub();
  877. end$
  878. algebraic procedure comm!-tab$
  879. if (lisp num!-dgen) geq 2 then
  880. begin integer nd; scalar v;
  881. nd:=lisp num!-dgen;
  882. write"The non-vanishing commutators of the finite subgroup";
  883. for i:=1:nd-1 do for j:=(i+1):nd do
  884. if(v:=comm(gen i,gen j)) neq 0 then
  885. if (lisp lengthf reval v) leq 60 then
  886. <<symbolic(!*list := nil); write"COMM(",i,",",j,"):= ",v>> else
  887. begin integer r,s,nt; scalar z; operator gen!*;
  888. nt:=lisp length(z:=numr simp reval v);
  889. r:=lisp maxl foreach x in z collect abs comfacn list x;
  890. if r=1 then r:=0 else r:=lisp flatsizec r;
  891. for i:=1:nt do gen!* i:=lisp prepf list nth(z,i);
  892. for l:=1:nt do
  893. begin
  894. symbolic(!*list := lengthf reval list('gen!*,l) geq 63);
  895. s:=lisp abs comfacn list nth(z,l);
  896. if r=0 then s:=0 else
  897. if s=1 then s:=-1 else s:=lisp flatsizec s;
  898. if l=1 then
  899. write"COMM(",i,",",j,"):=",lisp blanks(r-s+1),gen!* 1 else
  900. if minus!-f gen!* l then
  901. write lisp blanks(r-s+9),gen!* l else
  902. write lisp blanks(r-s+9)," + ",gen!* l;
  903. end;
  904. clear gen!*;
  905. end;
  906. symbolic(!*list := nil);
  907. end$
  908. symbolic procedure write!-dep u$
  909. if assoc(reval u,depl!*) then
  910. begin scalar v;
  911. v:=cdr assoc(u,depl!*);
  912. write car u,"(",cadr u,") depends on ";
  913. write caar v,"(",cadar v,")";
  914. foreach x in cdr v do write",",car x,"(",cadr x,")";
  915. terpri2();
  916. end$
  917. symbolic operator write!-dep$
  918. symbolic procedure find!-nfun$
  919. if not get('c,'klist) then 0 else
  920. maxl makeset foreach x in get('c,'klist) collect cadar x$
  921. endmodule;
  922. end;