1module sstools; 2 3create!-package('(sstools),nil); 4 5%******************************************************************** 6% * 7% The package SSTOOLS for the study of integrability of * 8% super-symmetric PDEs * 9% * 10% Author: Thomas Wolf and Eberhard Schruefer with contributions * 11% of Winfried Neun * 12% * 13%******************************************************************** 14 15global '(!*sstools!-loaded); 16 17!*sstools!-loaded := t; 18 19%----------------------- 20% The following module for non-commutative differentiation and 21% multiplication of fermionic forms is written by Eberhard Schruefer 2007. 22 23module dop$ 24 25fluid '(t_changes_parity s_changes_parity asymplis!* ncmp!*)$ 26 27switch t_changes_parity,s_changes_parity$ 28 29put('t_changes_parity,'simpfg, 30 '((t (rmsubs) (setq t_changes_parity t) (clear_t_changes_parity)) 31 (nil (rmsubs) (setq t_changes_parity nil) (clear_t_changes_parity))))$ 32 33put('s_changes_parity,'simpfg, 34 '((t (rmsubs) (setq s_changes_parity t) (clear_s_changes_parity)) 35 (nil (rmsubs) (setq s_changes_parity nil) (clear_s_changes_parity))))$ 36 37symbolic procedure clear_t_changes_parity()$ 38 begin 39 asymplis!* := clear_t_changes_power_rules asymplis!*; 40 end$ 41 42symbolic procedure clear_s_changes_parity()$ 43 begin 44 asymplis!* := clear_s_changes_power_rules asymplis!*; 45 end$ 46 47symbolic procedure clear_t_changes_power_rules u$ 48 if null u then u 49 else if smember('t,car u) then clear_t_changes_power_rules cdr u 50 else car u . clear_t_changes_power_rules cdr u$ 51 52symbolic procedure clear_s_changes_power_rules u$ 53 if null u then u 54 else if smember('s,car u) then clear_s_changes_power_rules cdr u 55 else car u . clear_s_changes_power_rules cdr u$ 56 57ncmp!* := t$ 58 59rlistat('(fermion))$ 60 61symbolic procedure fermion u$ fermion1 u$ 62 63symbolic procedure fermion1 u$ 64 <<flag(u,'fermionic); 65 for each j in u 66 do <<flag({j},'full); put(j,'simpfn,'simpfermion)>>; 67 for each j in u do noncom1 j>>$ 68 69symbolic procedure simpfermion u$ 70 begin scalar x,y,z; 71 if x := opmtch u then return simp x; 72 u := car u . for each j in cdr u collect aeval j; 73 y := mksq(u,1); 74 if domainp numr y then return y; 75 z := mvar numr y; 76 if atom z or null(car z eq car u) then return y; 77 if null atsoc(z,asymplis!*) 78 then asymplis!* := (z . 2) . asymplis!*; 79 return y 80 end$ 81 82rlistat('(boson))$ 83 84symbolic procedure boson u$ boson1 u$ 85 86symbolic procedure boson1 u$ 87 for each j in u do 88 <<flag({j},'bosonic); mkop j>>$ 89 90put('d,'simpfn,'simpdop)$ 91 92noncom1 'd$ 93 94symbolic procedure simpdop u$ 95 begin scalar x; 96 u := aeval car u . cdr u; 97 if x := opmtch ('d . u) then return simp x; 98% if x := opmtch u then return simp u; 99 x := simp cadr u; 100 if null numr x then return x; 101 if domainp denr x or null fermionicp mvar denr x 102 then return multsq(dopf(car u, numr x), 1 ./ denr x); 103 rederr "encountered non constant den in simpdop ..."; 104 end$ 105 106symbolic procedure dopf(n,u)$ 107 if domainp u then nil ./ 1 108 else addsq(addsq(multpq(lpow u,dopf(n,if fermionicp mvar u 109 then negf lc u 110 else lc u)), 111 multsq(dopp(n, lpow u),lc u ./ 1)),dopf(n, red u))$ 112 113symbolic procedure dopp(n, u)$ 114 % the following test that car u is a constant wrt d should be checked if 115 % it is complete. 116 if atom car u or null(flagp(caar u,'fermionic) or flagp(caar u,'bosonic) or 117 (caar u eq 'd) or (caar u eq 'df)) 118 then nil ./ 1 119 else if cdr u > 1 120 then multsq(cdr u ./ 1,multsq(simpexpt {car u,cdr u-1}, 121 simpdop {n,car u})) 122 else if caar u eq 'd 123 then if n = cadr car u 124 then simpdf {caddr car u,'x} 125 else if n > cadr car u 126 then negsq simpdop({cadr car u,{'d,n,caddr car u}}) 127 else mkdsq({'d,n,car u},1) 128 else mkdsq({'d,n,car u},1)$ 129 130symbolic procedure mkdsq(u,n)$ 131 begin scalar x,z; 132 if x := opmtch u then return simp x; 133 x := mksq(u,n); 134 if null numr x then return x; 135 z := mvar numr x; 136 if fermionicp z and null atsoc(z,asymplis!*) 137 then asymplis!* := (z . 2) . asymplis!*; 138 return x 139 end$ 140 141symbolic procedure dfdp(u,v,n)$ 142 begin scalar x; 143 x := simpdf {caddr u,v,n}; 144 if null numr x then return x; 145 if ((v eq 't) and t_changes_parity) or 146 ((v eq 's) and s_changes_parity) 147 then x := negsq x; 148 if domainp denr x or 149 null fermionicp mvar denr x 150 then return multsq(dopf(cadr u, numr x), 1 ./ denr x) 151 else rederr "fermioic den encountered in dfdp" 152 end$ 153 154put('d,'dfform,'dfdp)$ 155 156symbolic procedure fermionicp u$ 157 null atom u and (flagp(car u,'fermionic) or 158 (car u eq 'd and null fermionicp caddr u) or 159 (car u eq 'df and fermionicp caddr u) or 160 (car u eq 'df and 161 if fermionicp cadr u 162 then if memq('s,cddr u) then not s_changes_parity 163 else if memq('t,cddr u) then not t_changes_parity 164 else t 165 else if memq('s,cddr u) then s_changes_parity 166 else if memq('t,cddr u) then t_changes_parity 167 else nil))$ 168 169 170symbolic procedure sstools!-multfnc(u,v)$ 171 % Returns canonical product of U and V, with both main vars non- 172 % commutative. 173 begin scalar x,y; 174 x := multf(lc u,!*t2f lt v); 175 if null x then nil 176 else if not domainp x and mvar x eq mvar u 177 then x := addf(if null (y := mkspm(mvar u,ldeg u+ldeg x)) 178 then nil 179 else if y = 1 then lc x 180 else !*t2f(y .* lc x), 181 multf(!*p2f lpow u,red x)) 182 else if noncomp mvar u 183 then if ordop(mvar u,mvar x) 184 then x := !*t2f(lpow u .* x) 185 else if null(y := multf(!*p2f lpow u,lc x)) 186 then x := multf(!*p2f lpow u,red x) 187 else x := addf(!*t2f(lpow x .* 188 if fermionicp mvar u and 189 fermionicp mvar x then negf y 190 else y), 191 multf(!*p2f lpow u,red x)) 192 else x := multf(!*p2f lpow u,x) where !*!*processed=t; 193 return addf(x,addf(multf(red u,v),multf(!*t2f lt u,red v))) 194 end$ 195 196symbolic procedure difff(u,v)$ 197 %U is a standard form, V a kernel. 198 %Value is the standard quotient derivative of U wrt V. 199 % Allow for differentiable domains. 200 if atom u then nil ./ 1 201 else if atom car u 202 then (if diff!-fn then apply2(diff!-fn,u,v) else nil ./ 1) 203 where (diff!-fn = get(car u,'domain!-diff!-fn)) 204 else addsq(addsq(multpq(lpow u,difff(if (fermionicp mvar u and 205 (((v eq 's) and s_changes_parity) or 206 ((v eq 't) and t_changes_parity) or 207 fermionicp v 208 ) 209 ) 210 then negf lc u 211 else lc u,v)), 212 multsq(if cdr lpow u = 1 and (((v eq 's) and s_changes_parity) or 213 ((v eq 't) and t_changes_parity) or 214 fermionicp car lpow u) then diffdp(lpow u,v) 215 else diffp(lpow u,v),lc u ./ 1)), 216 difff(red u,v))$ 217 218symbolic procedure diffdp(u,v)$ 219 begin scalar x,y,z,w; 220 u := car u; 221 if null atom u and (x := get(car u,'dfform)) then return apply3(x,u,v,1); 222 if null depends(u,v) then return nil ./ 1; 223 if u eq v then return 1 ./ 1; 224 if eqcar(u,'df) and v eq 's 225 and t_changes_parity and s_changes_parity 226 and uneven_t_p cddr u 227 then z := -1 ./ 1 else z := 1 ./ 1; 228 if car u eq 'df 229 then if (x := find_sub_df(w := cadr u . merge!-ind!-vars(u,v), 230 get('df,'kvalue))) 231 then <<w := simp car x; 232 for each el in cdr x do 233 for i := 1:cdr el do 234 w := diffsq(w,car el); 235 return w>> 236 else u := 'df . w 237 else u := {'df,u,v}; 238 w := 1 ./ 1; 239 if x := opmtch u then return multsq(z,simp x); 240 x := mksq(u,1); 241 if null numr x then return x; 242 y := mvar numr x; 243 if fermionicp y and null atsoc(y,asymplis!*) 244 then asymplis!* := (y . 2) . asymplis!*; 245 return multsq(z,x) 246end$ 247 248symbolic procedure uneven_t_p u$ 249 if car u eq 't 250 then if cdr u then null(fixp cadr u and evenp cadr u) 251 else t 252 else cdr u and uneven_t_p cdr u$ 253 254bothtimes put('is_fermionic,'boolfn,'evalfermionicp)$ 255 256symbolic procedure evalfermionicp u$ 257fermionicfp numr simp u$ 258% begin scalar x; 259% x := numr simp u; 260% return fermionicfp x 261% end; 262 263symbolic procedure fermionicfp u$ 264if domainp u then nil 265%else (fermionicp mvar u and fermionicfp lc u) or fermionicfp red u; 266 else 267if fermionicp mvar u then not fermionicfp lc u 268 else fermionicfp lc u$ 269 270% The following function makes reordering of sstools expressions possible 271% (as needed by coeffn for example). Unfortunately some redefinitions 272% of standard REDUCE functions are necessary. 273 274symbolic procedure reordablep(u,v)$ 275 atom u or atom v or reordablekp(u,v)$ 276 277symbolic procedure reordablekp(u,v)$ 278 sstools_kernelp u or sstools_kernelp v$ 279 280symbolic procedure sstools_kernelp u$ 281 car u eq 'd or flagp(car u,'fermionic) or flagp(car u,'bosonic)$ 282 283symbolic procedure reordop(u,v)$ 284 if reordablep(u,v) then ordop(u,v) 285 else (!*ncmp and noncomp1 u and noncomp1 v) or ordop(u,v)$ 286 287symbolic procedure rmultpf(u,v); !*q2f simp!* prepf (u .* v .+ nil)$ 288 289 290endmodule$ 291 292%----------------------- further declarations and initializations 293 294%algebraic$ 295 296fermion1 '(f th)$ 297boson1 '(b)$ 298 299algebraic$ 300depend {f,b},t,x$ 301 302!#if (memq 'psl lispsystem!*) 303set_bndstk_size 50000$ 304!#endif 305 306lisp <<simplimit!* := 100000$ 307 !*noarg := nil % not done as 'off noarg$' to prevent compiler message 308 >>$ 309on dfprint$ 310 311lisp(!*allfac := nil)$ % instead of OFF ALLFAC in order to avoid difference 312 % between compiled and uncompiled version 313 314symbolic fluid '(x flist blist fblist sysfbl N_ fname_ fname_list nfct_ 315 use_new_crackout print_ print_more collect_sol sol_list 316 proc_list_ max_factor flin_ subst_1 pdelimit_1 old_history 317 homogen_ max_gc_short max_gc_red_len record_hist 318 !*time lin_test_const lines_written ineq_ html_out 319 !*t_changes_parity !*s_changes_parity)$ 320 321lrule1:={D(~n,df(~h,x )) => lin_test_const**3, 322 D(~n,df(~h,x,~m)) => lin_test_const**(2*m+1)}$ 323lrule2:={D(~n,~h) => lin_test_const, 324 df(~h,x) => lin_test_const**2, 325 df(~h,x,~n) => lin_test_const**(2*n)}$ 326 327subli:={}$ % the list of substitutions coming with the call of ssym 328symansatz:={}$ % the ansatz for equations and symmetry 329 330%----------------------- procedures 331 332if lisp(null(getd 'redfront_color) ) then 333symbolic procedure redfront_color(a)$ a$ 334 335%------- 336 337symbolic operator is_const$ 338symbolic procedure is_const(a)$ 339freeof(a,'f) and freeof(a,'b)$ 340 341%------- 342 343symbolic operator is_fermion$ 344symbolic procedure is_fermion(a)$ 345if pairp a then 346if car a = 'D then is_boson caddr a 347 else 348if car a = 'DF then if is_fermion cadr a then 349 if caddr a ='s then not !*s_changes_parity else 350 if caddr a ='t then not !*t_changes_parity else t 351 else 352 if caddr a ='s then !*s_changes_parity else 353 if caddr a ='t then !*t_changes_parity else nil 354 else 355if car a = 'EXPT then 356 if is_boson cadr a then nil else 357 if (caddr a = 1) then t else 358 if fixp caddr a then nil else 359 write"###### ERROR: UNDEFINED PARITY!! ######" 360 else 361if car a = 'TIMES then begin scalar n,h; 362 n:=0; 363 for each h in cdr a do 364 if is_fermion h then n:=add1 n; 365 return not evenp n 366 end 367 else 368if (car a = 'MINUS) or 369 (car a = 'QUOTIENT) then is_fermion cadr a else 370 % we assume, the denominator can not be fermionic 371if (car a = 'PLUS) then begin scalar h; 372 h:=cdr a; 373 while h and zerop car h do h:=cdr h; 374 return if h then is_fermion car h 375 else nil 376 end 377 else 378if not freeof(a,'f) then t else nil 379 else if not freeof(a,'f) then t else nil$ 380 381symbolic operator is_boson$ 382symbolic procedure is_boson(a)$ 383if pairp a then 384if car a = 'D then is_fermion caddr a 385 else 386if car a = 'DF then if is_boson cadr a then 387 if caddr a ='s then not !*s_changes_parity else 388 if caddr a ='t then not !*t_changes_parity else t 389 else 390 if caddr a ='s then !*s_changes_parity else 391 if caddr a ='t then !*t_changes_parity else nil 392 else 393if car a = 'EXPT then 394 if is_boson cadr a then t else 395 if (caddr a = 1) then nil else 396 if fixp caddr a then t else 397 write"###### ERROR: UNDEFINED PARITY!! ######" 398 else 399if (car a = 'MINUS) or 400 (car a = 'QUOTIENT) then is_boson(cadr a) else 401if car a = 'TIMES then begin scalar n,h; 402 n:=0; 403 for each h in cdr a do 404 if is_fermion h then n:=add1 n; 405 return evenp n 406 end 407 else 408if (car a = 'PLUS) then begin scalar h; 409 h:=cdr a; 410 while h and zerop car h do h:=cdr h; 411 return if h then is_boson car h 412 else t 413 end 414 else 415if not freeof(a,'b) then t else nil 416 else if not freeof(a,'b) then t else nil$ 417 418%------- 419 420symbolic operator ssini$ 421symbolic procedure ssini(N,nf,nb)$ 422% initializes N_, fblist and selects printing routines for f(i) and b(j) 423begin scalar j; 424 N_:=N; % assignung the global number of super symm generators N_ 425 % For customized printing:"$ 426 if nf=1 then put('f,'prifn,'myfpri) % to suppress (1) when printing f(1) 427 else put('f,'prifn,nil)$ % to print index of all f(i) 428 if nb=1 then put('b,'prifn,'myfpri) % to suppress (1) when printing b(1) 429 else put('b,'prifn,nil)$ % to print index of all b(i) 430 fblist:=append(for j:=1:nf collect {'f,j}, 431 for j:=1:nb collect {'b,j} )$ 432end$ 433 434%------- 435 436symbolic procedure presimplify(es,fl)$ 437% es,fl are algebraic lists of equations and functions 438begin scalar w,ineq_bak,nopowersbak,eqns,m,k,cpu,p,len_es$ 439 % Simplifying equations and dropping multiple ones 440 len_es:=sub1 length es$ 441 write len_es," equations result"$terpri()$ 442 cpu:=time(); 443 % write((k:=time())-cpu)/1000, 444 % " s: Now simplifying equations and dropping multiple ones"$terpri()$ 445 w:=nil; % list of equations consisting of only a coefficient 446 ineq_bak:=ineq_$ 447 ineq_:=nil$ % for simplifyterm 448 nopowersbak:=!*nopowers; 449 algebraic(off nopowers); 450 eqns:=nil$ 451 fl:=cdr reval fl; 452 m:=nil; % list of new equations consisting of only a coefficient 453 repeat << 454 if m then << 455 es:=algebraic(sub(lisp cons('LIST,for each k in m collect {'EQUAL,k,0}),lisp es))$ 456 % w:=append(m,w)$ 457 w:=nconc(m,w)$ 458 m:=nil; 459 >>$ 460 es:=cdr es$ % to get rid of the leading 'LIST 461 while es do 462 if zerop car es then es:=cdr es else << 463 k:=algebraic(factorize lisp car es); es:=cdr es$ 464 k:=cdr k; 465 p:=nil; 466 while k do << 467 if not freeoflist(cadr car k,fl) then p:=cons(cadr car k,p); 468 k:=cdr k; 469 >>$ 470 if null p then eqns:=cons(1,eqns) else << 471 if cdr p then p:=cons('TIMES,p) 472 else p:=simplifyterm(reval car p,fl); 473 if atom p then m :=union({p}, m) 474 else eqns:=union({p},eqns) 475 >> 476 >>$ 477 es:=cons('LIST,eqns); eqns:=nil 478 >> until null m; 479 480 write length w," coefficients found to be zero: "$listprint w$ terpri()$ 481 eqns:=nconc(es,w)$ 482 write reval (len_es+1-length eqns)," redundant equations have been droped."$ 483 terpri()$ 484 485 ineq_:=ineq_bak$ 486 if nopowersbak then algebraic(on nopowers)$ 487 write(time()-cpu)/1000," s for pre-simplification"$terpri()$ 488 489 return eqns 490end$ 491 492%------- 493 494symbolic procedure input_consistency_test(afwlist,abwlist)$ 495begin 496 if (length afwlist < 2) and 497 (length abwlist < 2) then 498 rederr "flist AND blist are both not assigned."$ 499 500 afwlist:=cdr afwlist; 501 while afwlist and fixp car afwlist and car afwlist>=0 do afwlist:=cdr afwlist; 502 if afwlist then if not fixp car afwlist then 503 rederr "The fermionic weight list contains not only numbers!" else % <0 504 if not fixp reval algebraic max_deg then 505 rederr "If a fermionic weight is < 0 then max_deg must be set!"$ 506 507 abwlist:=cdr abwlist; 508 while abwlist and fixp car abwlist and car abwlist>0 do abwlist:=cdr abwlist; 509 if abwlist then if not fixp car abwlist then 510 rederr "The bosonic weight list contains not only numbers!" else % <=0 511 if not fixp reval algebraic max_deg then 512 rederr "If a bosonic weight is < 1 then max_deg must be set!" 513end$ 514 515%------- 516 517symbolic procedure crossprodu(seta,setb,wgt,endprod)$ 518begin scalar g,k,setbcp; 519 if setb then 520 for each g in seta do << 521 setbcp:=setb; 522 k:=car g + caar setbcp; 523 while k<=wgt do << 524 if k=wgt then endprod:=cons(cons(k,append(cdr g,cdar setbcp)),endprod); 525 setbcp:=cdr setbcp; 526 if setbcp then k:=car g + caar setbcp 527 else k:=10000000 528 >> 529 >>$ 530 return endprod 531end$ 532 533%------- 534 535symbolic procedure rhs_term_list(N,allf,allb,maxwgt,linsub,forbid,verbose)$ 536% generation of a list of all fermionic terms and a list of all 537% bosonic terms, all of weight maxwgt 538begin scalar g,h,j,k,minwgt,maxdeg,maxtermwgt,allv,newv,maxpow,maxwgtm2, 539 fprod,bprod,allprod,newprod,linonly,p,q,r,s$ 540 541 % As weights can be negative, at first we compute the most negative weight 542 % of a partial product. For that we do not need to consider derivatives as 543 % these have all higher weights than the undifferentiated variables. 544 545 minwgt:=0; maxdeg:=algebraic max_deg; allv:=append(allf,allb); 546 for each g in allv do 547 if car g<0 then minwgt:=minwgt+maxdeg*car g; 548 maxtermwgt:=maxwgt-minwgt; % >=maxwgt as minwgt<=0 549 550 % Add to allv for every field its x-derivatives as high as possible 551 % and sort the entries by their weight 552 553 maxwgtm2:=maxtermwgt-2; 554 % newv:=nil; through initialization 555 for each g in allv do << 556 while (car g <= maxwgtm2) and not member(reval {'DF,cadr g,'x},forbid) do << 557 h:=2 + car g; 558 g:={h,{'DF,cadr g,'x}}; % not: reval {'DF,cadr g,'x} 559 newv:=cons(g,newv) 560 >> 561 >>$ 562 % allv:=idx_sort append(allv,newv); 563 allv:=append(allv,newv); 564 565 % Add to allv for each of its elements all possible combinations 566 % of d(i,..) derivatives 567 for j:=1:N do << 568 newv:=nil; 569 for each g in allv do 570 if car g < maxtermwgt then 571 newv:=cons({add1 car g,{'D,j,cadr g}},newv); 572 % allv:=idx_sort append(allv,newv); 573 allv:=append(allv,newv); 574 >>$ 575 576 % generation of all possible products of elements of allv 577 allprod:={{0,nil}}; 578 if verbose then << 579 write"maxwgt=",maxwgt$terpri()$ 580 q:=length allv; 581 p:=0; 582 s:=0; 583 >>$ 584 for each g in allv do << 585 if verbose then << 586 p:=add1 p; 587 write p,"(",q,"): ",g$terpri() 588 >>$ 589 590 if null linsub then linonly:=nil else 591 if smemberl(linsub,cadr g) then linonly:=t 592 else linonly:=nil; 593 maxpow:=if linonly or is_fermion cadr g then 1 else 594 if car g < 1 then maxdeg else 10000000; 595 newprod:=nil; 596 if verbose then r:=0; 597 for each h in allprod do 598 if null linonly or null cadr h then << 599 % i.e. either the new factor cadr g is not linonly and 600 % can occur arbitrarily often in arbitrary products, or 601 % the factors do not contain a linonly factor yet 602 k:=car h + car g; 603 j:=1$ 604 while (j <= maxpow) and (k <= maxtermwgt) do << 605 h:=cons(k,cons(linonly or cadr h,cons(cadr g,cddr h))); 606 newprod:=cons(h,newprod); 607 k:=car h + car g; j:=add1 j 608 >>; 609 if verbose then r:=r+j-1 610 >>$ 611 if verbose then s:=s+r; 612 allprod:=nconc(allprod,newprod) 613 >>$ 614 615 % sorting of the products into fermionic and bosonic products 616 allprod:=cdr allprod; % drop of {{0,nil}} 617 while allprod do << 618 g:=car allprod; allprod:=cdr allprod; 619 if (car g = maxwgt) and (null linsub or cadr g) then 620 if is_fermion cons('TIMES,cdr g) then 621 fprod:=cons(if null cdddr g then caddr g 622 else cons('TIMES,cddr g),fprod) 623 else 624 bprod:=cons(if null cdddr g then caddr g 625 else cons('TIMES,cddr g),bprod) 626 >>$ 627 628 if verbose then <<write length fprod," fermionic and ", 629 length bprod," bosonic monomials"$terpri()>>$ 630 return {fprod,bprod} 631 632end$ 633 634%------- 635 636symbolic procedure gen_eqn(awlist,rhslist,eqn_list,fl,fermionic,tname)$ 637begin scalar h,w,rlcp,rhs,nf,print_bak,svar$ 638 svar:=reval tname$ 639 while awlist do << 640 % an equation df(h,t)=... is generated for each car of element h of awlist 641 h:=car awlist; awlist:=cdr awlist; 642 643 % w is the weight of h 644 w:=car h; h:=cadr h; 645 646 % find the proper rhs for the weight w of h 647 rlcp:=rhslist; 648 while caar rlcp neq w do rlcp:=cdr rlcp; 649 rhs:=if fermionic then cadar rlcp 650 else caddar rlcp; 651 % add a constant factor to each term in rhs 652 print_bak:=print_; print_:=nil; % to avoid printing the new constants 653 rhs:=for each r in rhs collect <<nf :=newfct(fname_,nil,nfct_)$ 654 fl :=cons(nf,fl); 655 nfct_:=add1 nfct_$ 656 {'TIMES,nf,r}>>$ 657 print_:=print_bak; 658 rhs:=if null rhs then 0 659 else if cdr rhs then cons('PLUS,rhs) 660 else car rhs; 661 eqn_list:=cons({'EQUAL,{'DF,h,svar},rhs},eqn_list) 662 >>$ 663 664 return {eqn_list,fl} 665end$ 666 667%------- 668 669symbolic procedure sspol(N,sysfbl,afwlist,abwlist,difforder,tname, 670 linsub,forbid,flip_par,verbose)$ 671% returns {eqn_list,fl} where 672% eqn_list is the algebraic list of equations/symmetries and 673% fl is the algebraic list of undetermined constants 674% tname is the lhs differentiation variable 675begin scalar g,h,i,fwlist,bwlist,allf,allb,w_list,rhs_list,eqn_list,fl,awlist$ 676 677 fwlist:=cdr reval afwlist$ 678 bwlist:=cdr reval abwlist$ 679 680 % Add to each field its weight as a prefix to create allf and allb 681 for i:=1:(length flist) do allf:=cons({nth(fwlist,i),nth(flist,i)},allf); 682 for i:=1:(length blist) do allb:=cons({nth(bwlist,i),nth(blist,i)},allb); 683 684 % Generation of the final form of the equations 685 686 % only those weights for which there is an eqn to generate 687 % making a list of all different weights 688 w_list:=union(fwlist,union(bwlist,nil))$ 689 690 % generating the right-hand-side term-list for the different weights 691 % Because expressions are generated according to their weight and are 692 % partitioned only afterwards into fermionic and bosonic, we would not save 693 % computation when generating only fermionic expressions of specific weight 694 % and afterwards the bosonic expressions with specific weight 695 rhs_list:=for each w in w_list collect 696 cons(w,rhs_term_list(N,allf,allb,w+difforder,linsub,forbid,verbose))$ 697 % new weight(t) based on weight(x)=2 698 699 % generating all bosonic equations 700 for each h in sysfbl do 701 if car h = 'b then << g:=allb; while g and cadar g neq h do g:=cdr g$ 702 awlist:=cons(car g,awlist) 703 >>$ 704 nfct_:=1$ 705 % whether difforder is odd or even has no influence on parity (ferm/bos) 706 h:=gen_eqn(awlist,rhs_list,eqn_list,fl,flip_par,tname); 707 708 % generating all fermionic equations 709 awlist:=nil$ 710 for each h in sysfbl do 711 if car h = 'f then << g:=allf; while g and cadar g neq h do g:=cdr g$ 712 awlist:=cons(car g,awlist) 713 >>$ 714 % whether difforder is odd or even has no influence on parity (ferm/bos) 715 h:=gen_eqn(awlist,rhs_list,car h,cadr h,not flip_par,tname); 716 717 return {cons('LIST,car h),cadr h} 718end$ 719 720%------- 721 722symbolic procedure print_list_of_lists(li)$ 723begin scalar k$ 724 for each k in li do << 725 write"{"$ 726 while k do <<write car k$k:=cdr k$if k then write",">>$ 727 write"}"$terpri()$ 728 >>$ 729end$ 730 731%------- 732 733symbolic operator listine$ 734symbolic procedure listine(N,sys,sym,fl,ineql,non_lin_test,save_lists)$ 735% generation of conditions for the unknown coefficients 736begin scalar eqn,ineql,cnd,h,k,fb,flcp,nfl,evolist,rs$ 737 fl:=cdr fl$ 738 ineql:=for each h in cdr ineql collect cdr h$ 739 740 evolist:=for each h in append(cdr sys,cdr sym) collect caddr h$ 741 k:=gensym()$ 742 743 % The online data base was computed with: 744 % % none of the equations and symmetries should vanish identically 745 % for each eqn in evolist do << 746 % if null (h:=smemberl(fl,eqn)) then h:={0}; 747 % ineql:=cons(h,ineql) 748 % >>$ 749 % %write"ineql1=",ineql$terpri()$ 750 751 % none of the equations should vanish identically 752 for each eqn in cdr sys do << 753 if null (h:=smemberl(fl,caddr eqn)) then h:={0}; 754 ineql:=cons(h,ineql) 755 >>$ 756 757 % at least one of the symmetries must not vanish identically 758 h:=for each h in cdr sym collect caddr h$ 759 if null (h:=smemberl(fl,h)) then h:={0}; 760 ineql:=cons(h,ineql)$ 761 762 % None of the evolution equations should involve only one field, 763 % the system should not be triangular. 764 if cdr fblist then % if more than one field 765 for each eqn in cdr sys do << 766 rs:=caddr eqn; % rhs of eqn 767 fb:=delete(cadr cadr eqn,fblist); % one of the functions fb must occur in rs 768 if pairp rs and (car rs = '!*SQ) 769 then rs:={'!*SQ, subsq(cadr rs, for each h in fb collect (h . 0)), t} 770 else for each h in fb do rs:=subst(0,h,rs); 771 h:=setdiff(smemberl(fl,caddr eqn),smemberl(fl,reval rs)); 772 if null h then ineql:=cons({0},ineql) 773 else ineql:=cons(h,ineql) 774 >>$ 775 776 % At least one evolutionary equation should be non-linear. 777 if non_lin_test then << 778 nfl:=nil; % we drop all constants from fl which are coefficient to 779 % a linear term 780 for each eqn in evolist do << 781 flcp:=smemberl(fl,eqn); 782 for each h in fblist do eqn:=subst({'TIMES,k,h},h,eqn); 783 h:=smemberl(flcp,coeffn(reval eqn,k,1)); 784 flcp:=setdiff(flcp,h); 785 nfl:=append(nfl,flcp) 786 >>$ 787 if null nfl then nfl:={0}; % contradiction issued 788 ineql:=cons(nfl,ineql); 789 >>$ 790 791 % at least one equation should contain a fermion field or a D_i derivative 792 if null flist then << % otherwise no chance to get a purely classical equation 793 flcp:=fl; % we drop all constants from fl which are coefficient to 794 % a term without D(i,..) and without fermion 795 algebraic << 796 drop_fd_rules:={d(~i,~j) => 0, f(~i) => 0}; 797 let drop_fd_rules; 798 >>$ 799 for each eqn in evolist do << 800 h:=smemberl(fl,reval eqn); 801 flcp:=setdiff(flcp,h) 802 >>$ 803 algebraic ( clearrules drop_fd_rules )$ 804 if null flcp then flcp:={0}; % contradiction issued 805 ineql:=cons(flcp,ineql); 806 %write"ineql5=",ineql$terpri()$ 807 >>$ 808 809 % N>1 makes only sense if all D_1,..,D_N occur which is required next 810 if N>1 then for h:=1:N do << 811 flcp:=fl; % we drop all constants from fl which are coefficient to 812 % a term without D(h,..) 813 algebraic <<k:=lisp h;drop_fd_rules:={d(k,~j) => 0}; let drop_fd_rules>>; 814 for each eqn in evolist do << 815 k:=smemberl(fl,reval eqn); 816 flcp:=setdiff(flcp,k) 817 >>$ 818 algebraic ( clearrules drop_fd_rules )$ 819 if null flcp then flcp:={0}; % contradiction issued 820 ineql:=cons(flcp,ineql); 821 >>$ 822 823 % dropping redundand conditions 824 repeat << 825 rs:=ineql; 826 while rs do << 827 k:=delete(car rs,ineql); 828 while k and not_included(car rs,car k) do k:=cdr k; 829 if k then << % i.e. cnd is fully included in car k 830 ineql:=delete(car k,ineql); 831 rs:=ineql 832 >> else rs:=cdr rs 833 >> 834 >> until null rs$ 835 836 if ineql then << 837 terpri()$ 838 write"From each of the following lists at least one element must be non-zero:"$ 839 terpri()$terpri()$ 840 print_list_of_lists(ineql)$ 841 terpri()$ 842 >>$ 843 if save_lists then << 844 out "inelist"$ 845 print_list_of_lists(ineql)$ 846 shut "inelist"$ 847 >>$ 848 849 ineql:=for each cnd in ineql collect 850 if null cdr cnd then car cnd 851 else cons('PLUS,for each h in cnd collect {'TIMES,h,gensym()})$ 852 853 return cons('LIST,ineql) 854 855end$ 856 857%------- 858 859algebraic procedure power_parti(eqns)$ 860begin scalar w,afblist,equa,ff,ex,a,p,parti,maxpow,h$ 861 afblist:=lisp(cons('LIST,fblist))$ 862 w:=lisp lin_test_const$ 863 maxpow:=0$ 864 h:=for each equa in eqns collect << 865 ff:=part(equa,1);ex:=part(equa,2); 866 for each a in afblist do ex:=sub(a=a*w,ex)$ 867 parti:=rest coeff(ex,w); 868 if hipow!*>maxpow then maxpow:=hipow!*; 869 ff=parti 870 >>; 871 872 % Filling up each element of h with zero's so that all have same length 873 if length h > 1 then 874 h:=for each equa in h collect part(equa,1)=<< 875 le:=length part(equa,2); 876 if le=maxpow then part(equa,2) 877 else <<ex:={};for g:=1:(maxpow-le) do ex:=cons(0,ex)$ 878 append(part(equa,2),ex)>> 879 >>$ 880 return {maxpow,h} 881end$ 882 883%------- 884 885symbolic procedure add_df_rules_to_D_rule(subl)$ 886% subl is an algebraic list of rules 887% adds to each substitution rule D(i,A)=>.. all differential consequences 888if null cdr subl then subl else 889begin scalar g,h$ 890 g:=cdr subl; 891 for each h in g do 892 if pairp h and car h='REPLACEBY and pairp cadr h and caadr h='D and 893 pairp caddr cadr h and ((caaddr cadr h='f) or (caaddr cadr h='b)) then << 894 % In algebraic mode rules are not properly printed, i.e. without brackets 895 % around sums to be differentiated, so they may look wrong when being right 896 subl:=cons('LIST,union( 897 {{'REPLACEBY,{'DF,caddr cadr h,'x},reval {'D,cadr cadr h,caddr h}}, 898 {'REPLACEBY,{'DF,caddr cadr h,'x,{'~,'n}}, 899 {'D,cadr cadr h,{'DF,caddr h,'x,{'DIFFERENCE,'n,1}}}}, 900 {'REPLACEBY,{'D,cadr cadr h,{'DF,caddr cadr h,{'~,'x}}}, 901 reval {'DF,caddr h,'x}}, 902 {'REPLACEBY,{'D,cadr cadr h,{'DF,caddr cadr h,{'~,'x},{'~,'n}}}, 903 {'DF,caddr h,'x,'n}}, 904 {'REPLACEBY,{'D,cadr cadr h,{'DF,caddr cadr h,{'~,'t},{'~,'x},{'~,'n}}}, 905 {'DF,caddr h,'t,'x,'n}}}, 906 cdr subl)) 907 >>$ 908 return subl 909end$ 910 911%------- 912 913algebraic procedure sieve(inpu,wgl)$ 914% inpu={sym,cl} 915% sym is an algebraic list of equations df(f(i),s)=.., df(b(j),s=.., 916% cl is an algebraic list of coefficients in sym 917% wgl is a list {sw,afwlist,abwlist} 918% It returns a list (new symmetry, list of coefficients} 919begin scalar sym,cl,sw,fli,bli,su,n,h,newsy,to_drop,sy$ 920 cl:=second inpu; 921 sw:=first wgl; 922 fli:=second wgl; 923 bli:=third wgl; 924 925 % creating a substitution rule 926 su:={}; 927 n:=0$for each h in fli do <<n:=n+1;su:=cons(f(n)=f(n)*lin_test_const**h,su)>>$ 928 n:=0$for each h in bli do <<n:=n+1;su:=cons(b(n)=b(n)*lin_test_const**h,su)>>$ 929 sym:=sub(su,first inpu); 930 let lrule1; sym:=sym; clearrules lrule1; 931 let lrule2; sym:=sym; clearrules lrule2; 932 933 newsy:={}; 934 to_drop:={}$ 935 for each sy in sym do << 936 coeff(lhs sy,lin_test_const); 937 n:=hipow!*+sw; 938 if n<0 then <<h:=lin_test_const**(-n)*rhs sy; 939 newsy:=cons(coeffn(h,lin_test_const,n),newsy)>> 940 else newsy:=cons(coeffn(rhs sy,lin_test_const,n),newsy); 941 >>; 942 su:={}; 943 for each h in cl do if freeof(newsy,h) then << 944 su:=cons(h=0,su); 945 to_drop:=cons(h,to_drop) 946 >>$ 947 948 return {sub(su,first inpu),lisp cons('LIST,setdiff(cdr cl,cdr to_drop))} 949end$ 950 951%------- 952 953symbolic procedure non_t_lhs_dvs(pdes)$ 954% listing lhs's of equations and substitutions which are no t-derivatives 955begin scalar h,forbid$ 956 for each h in pdes do 957 if pairp h and 958 ((car h = 'EQUAL) or (car h = 'REPLACE)) and 959 (pairp cadr h) and 960 ((caadr h = 'DF) or (caadr h = 'D) or (caadr h = 'F) or (caadr h = 'B)) and 961 freeof(cadr h,t) then << 962 forbid:=cons(cadr h,forbid); 963 if caadr h = 'D then forbid:=cons({'DF,caddr cadr h,'x},forbid) else 964 if caadr h = 'df then for g:=1:N_ do forbid:=cons({'D,g,cadr h},forbid) 965 >>$ 966 return forbid 967end$ 968 969%------- 970 971algebraic procedure ssym(N,tw,sw,afwlist,abwlist,eqnlist,fl,inelist,mode)$ 972% uses the global algebraic variable x 973 974% N .. the number of superfields theta_i 975% tw .. 2 times the differential order of the equation = weight(t) 976% sw .. 2 times the differential order of the symmetry = weight(s) 977% afwlist .. (algebraic) list of weights of the fermion fields 978% f(1), f(2), ... 979% abwlist .. (algebraic) list of weights of the boson fields 980% b(1), b(2), ... 981% The number of elements in the last two lists determines the 982% number of fermion and boson fields. 983% eqnlist .. in the nonlinear case a list of extra conditions on the 984% undetermined coefficients where conditions a3=.. are executed 985% instantly, these and any, expressions are added to equations 986% when calling crack 987% in the linear case the system in form of replacements and the 988% linearized system in form of equations 989% fl .. extra unknowns in eqnlist to be determined 990% inelist .. a list, each element of it is a list with at least one of 991% its elements being non-zero 992% mode .. list of flags: 993% init : only initialization of global data 994% plain_com : direct computation of the commutator (for tests) 995% power_split_com : alternative power splitting of commutator, 996% (not if substitutions '='> are present) 997% % const_coeff_in_eqn : have only constant coefficients in eqn. 998% % const_coeff_in_sym : have only constant coefficients in sym. 999% zerocoeff : all coefficients = 0 which are not listed in ine 1000% % verbose : extra output for tracing and debugging 1001% interactive : whether the solution process is interactive 1002% % no_t_in_eqn : no explicit t-derivatives occur in equations 1003% % no_t_in_sym : no explicit t-derivatives occur in symmetries 1004% filter: extra homogeneity weights as given in the list"$ 1005% hom_wei have to be satisfied by the symmetry"$ 1006% lin: find symmetry for all linear fields, 1007% if null spar (bosonic s) then 1008% bosonic lin fields are b(nb+1)..b(2*nb) 1009% fermionic lin fields are f(nf+1)..b(2*nf) 1010% else 1011% bosonic lin fields are b(nb+1)..b(nf+nb) 1012% fermionic lin fields are f(nf+1)..b(nf+nb) 1013% tpar: t/nil, whether time variable t changes parity"$ 1014% spar: t/nil, whether symmetry variable s changes parity"$ 1015% log : t/nil, whether files drvlist, evolist, unolist, inelist 1016% are generated to hold data, for example, for 1017% automatic web page generation 1018% 1019begin scalar g,h,k,cpu,gti,fbno,psys,psym,msysp,msymp,totpow,syspow,sympow, 1020 newcd,afblist,sublist,zerocoeff,fl_e,non_lin_test,do_ine_test, 1021 init,subl,subl2,sys,sym,ls,rs,linsub,filter,forbid,rhssyl,sycon, 1022 nw,w,np,p,nq,q,psycon,cn,lp,verbose,plain_com,power_split_com, 1023 msgbak,interactive,nfr,nfe,nbr,nbe,nf,nb,lhslist,save_lists$ 1024 backup_reduce_flags()$ 1025 lisp << 1026 record_hist:=nil; 1027 if !*time then <<cpu:=time()$ gti:=gctime()>>$ 1028 1029 eqnlist:=cdr eqnlist; 1030 fl :=cdr fl; if fl then homogen_:=nil else homogen_:=t; 1031 mode :=cdr mode; 1032 1033 if member('init,mode) then init:=t else init:=nil; 1034 % if member('const_coeff_in_eqn,mode) then cc_eqn:=t else cc_eqn:=nil; 1035 % if member('const_coeff_in_sym,mode) then cc_sym:=t else cc_sym:=nil; 1036 if member('zerocoeff,mode) then zerocoeff:=t else zerocoeff:=nil; 1037 if member('plain_com,mode) then plain_com:=t else plain_com:=nil; 1038 if member('power_split_com,mode) then power_split_com:=t 1039 else power_split_com:=nil; 1040 % if member('term_split_com,mode) then term_split_com:=t 1041 % else term_split_com:=nil; 1042 if member('verbose,mode) then verbose:=t else verbose:=nil; 1043 if member('interactive,mode) then interactive:=t else interactive:=nil; 1044 % if member('no_t_in_eqn,mode) then no_t_in_eqn:=t else no_t_in_eqn:=nil; 1045 % if member('no_t_in_sym,mode) then no_t_in_sym:=t else no_t_in_sym:=nil; 1046 if member('filter,mode) then filter:=t else filter:=nil; 1047 if member('tpar,mode) then algebraic(on t_changes_parity) 1048 else algebraic(off t_changes_parity); 1049 if member('spar,mode) then algebraic(on s_changes_parity) 1050 else algebraic(off s_changes_parity); 1051 if member('log,mode) then save_lists:=t$ 1052 1053 if member('lin,mode) then << % @#@# 1054 nf:=length afwlist - 1; 1055 nb:=length abwlist - 1; 1056 if member('spar,mode) then 1057 if nf neq nb then <<write"In the case of spar #(boson fields)=#(fermion fields)", 1058 " which is not the case."; forbid:=t>> 1059 else % to do more tests one would have to know/find the numbers 1060 % of b/f fields before linearization 1061 else << % null member('spar,mode) 1062 if null evenp nf then <<write"If no spar then #(fermion fields) needs to be even"$ 1063 forbid:=t>>; 1064 if null evenp nb then <<write"If no spar then #(boson fields) needs to be even"$ 1065 forbid:=t>>; 1066 for each g in eqnlist do if pairp g then lhslist:=cons(cadr g,lhslist); 1067 nf:=nf/2; nb:=nb/2; 1068 for each g in lhslist do 1069 1070if car g = 'df then 1071if caadr g = 'f then if ((cadadr g <= nf) and not member({'df,{'f,cadadr g + nf},caddr g},lhslist)) or 1072 ((cadadr g > nf) and not member({'df,{'f,cadadr g - nf},caddr g},lhslist)) then << 1073 write"The counterpart of ",g," is missing on a left hand side."$ forbid:=t>> else 1074 else 1075if caadr g = 'b then if ((cadadr g <= nb) and not member({'df,{'b,cadadr g + nb},caddr g},lhslist)) or 1076 ((cadadr g > nb) and not member({'df,{'b,cadadr g - nb},caddr g},lhslist)) then << 1077 write"The counterpart of ",g," is missing on a left hand side."$ forbid:=t>> else 1078 else 1079 else 1080if car g = 'd then 1081if caaddr g = 'f then if ((car cdaddr g <= nf) and not member({'d,cadr g,{'f,car cdaddr g + nf}},lhslist)) or 1082 ((car cdaddr g > nf) and not member({'d,cadr g,{'f,car cdaddr g - nf}},lhslist)) then << 1083 write"The counterpart of ",g," is missing on a left hand side."$ forbid:=t>> else 1084 else 1085if caaddr g = 'b then if ((car cdaddr g <= nb) and not member({'d,cadr g,{'b,car cdaddr g + nb}},lhslist)) or 1086 ((car cdaddr g > nb) and not member({'d,cadr g,{'b,car cdaddr g - nb}},lhslist)) then << 1087 write"The counterpart of ",g," is missing on a left hand side."$ forbid:=t>> else 1088 else 1089 else 1090 >>; % of null member('spar,mode) 1091 if forbid then <<forbid:=nil;rederr "">>$ 1092 >>$ 1093 1094 % In the linearized system not half of the equations need to 1095 % be = and half => because there may be nonlocal potential variables 1096 % introduced through a => relation 1097 if nil then 1098 if member('lin,mode) then << % @#@# 1099 % determine nf and nb 1100 nfr:=0; nbr:=0; % number of f(i) and b(i) in replaceby 1101 nfe:=0; nbe:=0; % number of f(i) and b(i) in equal 1102 1103 for each g in eqnlist do 1104 if pairp g then 1105 if (car g='equal) then 1106 if not freeof(cadr g,'f) then nfe:=nfe+1 else 1107 if not freeof(cadr g,'b) then nbe:=nbe+1 else 1108 else 1109 if (car g='replaceby) then 1110 if not freeof(cadr g,'f) then nfr:=nfr+1 else 1111 if not freeof(cadr g,'b) then nbr:=nbr+1 else 1112 else $ 1113 1114%check nfr+nfe = length afwlist - 1 1115%check nbr+nbe = length abwlist - 1 1116 1117write"nfr= ",nfr," nfe=",nfe, 1118 "nbr= ",nbr," nbe=",nbe$ 1119 1120 if ( member('spar,mode) and (nfr=nbe) and (nbr=nfe)) or 1121 (not member('spar,mode) and (nfr=nfe) and (nbr=nbe)) then 1122 else << 1123 if member('spar,mode) then << 1124 if nfr neq nbe then <<write "For spar 1125 the number of '=>' relations with a fermion on the lhs should be equal 1126 the number of '=' relations with a boson on the lhs."$ terpri() 1127 >>$ 1128 if nbr neq nfe then <<write "For spar 1129 the number of '=>' relations with a boson on the lhs should be equal 1130 the number of '=' relations with a fermion on the lhs."$ terpri() 1131 >> 1132 >> else 1133 1134 if not member('spar,mode) then << 1135 if nfr neq nfe then <<write "For missing spar 1136 the number of '=>' relations with a fermion on the lhs should be equal 1137 the number of '=' relations with a fermion on the lhs."$ terpri() 1138 >>$ 1139 if nbr neq nbe then <<write "For missing spar 1140 the number of '=>' relations with a boson on the lhs should be equal 1141 the number of '=' relations with a boson on the lhs."$ terpri() 1142 >> 1143 >>$ 1144 rederr "" 1145 >> 1146 >>; 1147 1148 % Initialization of fermionic and bosonic fields + their printing 1149 input_consistency_test(afwlist,abwlist)$ 1150 N_:=N; % to avoid printing the index of D in crack_out for N=1 1151 flist := for g:=1:(length afwlist-1) collect {'f,g}; 1152 if flist and length flist=1 then put('f,'prifn,'myfpri) 1153 else put('f,'prifn,nil); 1154 blist := for g:=1:(length abwlist-1) collect {'b,g}; 1155 if blist and length blist=1 then put('b,'prifn,'myfpri) 1156 else put('b,'prifn,nil); 1157 fblist:=append(flist,blist) 1158 >>$ 1159 1160 afblist:=lisp(cons('LIST,fblist))$ 1161 for each g in afblist do depend g,x,t$ 1162 1163 lisp << 1164 1165 % listing lhs's of equations and substitutions which are no t-derivatives 1166 forbid:=non_t_lhs_dvs(eqnlist)$ 1167 1168 terpri()$ 1169 write"####### This is the case N",N,"f",length flist,"b", 1170 length blist,"t",tw,"s",sw,"w"$ 1171 for each g in append(cdr afwlist,cdr abwlist) do write g; 1172 write". #######"$terpri()$ 1173 % specifying the names of constants in the equations 1174 if fname_list then <<fname_:=car fname_list; 1175 fname_list:=cdr fname_list>> 1176 else fname_:='p$ 1177 1178 % generating the equations 1179 % Has an ansatz for the system been made? 1180 if eqnlist and 1181 (pairp car eqnlist ) and 1182 (((caar eqnlist='EQUAL ) and 1183 (pairp cadar eqnlist ) and 1184 (caadar eqnlist = 'DF ) ) or 1185 ((caar eqnlist='REPLACEBY ) and 1186 (pairp cadar eqnlist ) and 1187 ((caadar eqnlist = 'DF) or 1188 (caadar eqnlist = 'D ) ) ) ) then << 1189 g:=nil; 1190 for each h in eqnlist do if pairp h then 1191 if pairp h and 1192 car h='EQUAL and 1193 pairp cadr h and 1194 caadr h='DF and 1195 reval caddr cadr h ='t then << 1196 sys:=cons(h,sys); 1197 subl2:=cons({'REPLACEBY,cadr h,caddr h},subl2) % substitutions based on sys 1198 >> else 1199 if car h='REPLACEBY and 1200 pairp cadr h and 1201 ((caadr h='DF) or 1202 (caadr h='D) ) then subl:=cons(h,subl) 1203 else g:=cons(h,g)$ 1204 eqnlist:=g$ 1205 1206 sysfbl:=reverse for each h in sys collect cadadr h; 1207 1208 % simplifying the rhs of substitutions 1209 if verbose then <<write"Simplifying the substitutions."$terpri()>>$ 1210 subl2:=cons('LIST,append(subl,subl2)); 1211 algebraic(let subl2); 1212 subl:=for each h in subl collect {'REPLACEBY,cadr h,algebraic(lisp(caddr h))}}$ 1213 subl:=cons('LIST,reverse subl); 1214 algebraic(clearrules subl2); 1215 subl2:=nil; 1216 1217 sys:=cons('LIST,reverse sys); 1218 fl_e:=fl$ 1219 fl:=nil 1220 >> else << 1221 do_ine_test:=t$ 1222 if null eqnlist then non_lin_test:=t$ 1223 sysfbl:=append(flist,blist); 1224 subl2:=cons('LIST,append(subl,subl2)); 1225 if verbose then <<write"Formulating rhs's of the system"$terpri()>>$ 1226 h:=sspol(N,sysfbl,afwlist,abwlist,tw,algebraic t,linsub, 1227 forbid,!*t_changes_parity,verbose)$ 1228 sys:=car h; 1229 fl_e:=cadr h; 1230 subl:={'LIST} 1231 >>$ 1232 1233 algebraic<<nodepend f,s; nodepend b,s>>$ 1234 for each h in fblist do algebraic(nodepend h,s); 1235 for each h in sysfbl do algebraic( depend h,s); 1236 1237 % add to each substitution rule D(i,A)=>.. all differential consequences 1238 subl:=add_df_rules_to_D_rule subl$ 1239 1240 % specifying the names of constants in the symmetries 1241 if fname_list then <<fname_:=car fname_list; 1242 fname_list:=cdr fname_list>> 1243 else fname_:='q$ 1244 % generating the symmetry 1245 if verbose then <<write"Formulating rhs's of the symmetry"$terpri()>>$ 1246 h:=sspol(N,sysfbl,afwlist,abwlist,sw,algebraic s,linsub, 1247 forbid,!*s_changes_parity,verbose)$ 1248 1249 % applying a homogeneity filter 1250 if filter and pairp algebraic hom_wei 1251 and car (algebraic hom_wei) = 'LIST then << 1252 if verbose then <<write"Applying an extra homogeneity filter"$terpri()>>$ 1253 h:={'LIST,car h,cons('LIST,cadr h)}; 1254 algebraic(for each g in hom_wei do h:=sieve(h,g))$ 1255 h:={cadr h,cdaddr h} % converting to lisp list of two lisp lists 1256 >>$ 1257 sym:=car h$ 1258 flin_:=cadr h$ h:=nil; 1259 1260 % Extracting substitutions from eqnlist to be done instantly 1261 g:=eqnlist; eqnlist:=nil; 1262 for each h in g do 1263 if pairp h and car h = 'EQUAL then << 1264 sublist:=cons(h,sublist)$ 1265 eqnlist:=cons({'DIFFERENCE,cadr h,caddr h},eqnlist) 1266 >> else eqnlist:=cons(h,eqnlist)$ 1267 1268 if zerocoeff then 1269 for each h in append(fl_e,flin_) do 1270 if freeof(inelist,h) then sublist:=cons({'EQUAL,h,0},sublist) 1271 else fl:=cons(h,fl) 1272 else fl:=append(fl_e,flin_); 1273 fl_e:=nil$ 1274 fl:=cons('LIST,fl)$ 1275 sublist:=cons('LIST,sublist)$ 1276 eqnlist:=cons('LIST,eqnlist)$ 1277 >>$ 1278 1279 if sublist then << 1280 sys:=sub(sublist,sys)$ 1281 sym:=sub(sublist,sym) 1282 >>$ 1283 if lisp(save_lists) then << 1284 off nat$ 1285 out "drvlist"$ 1286 for each g in sys do write lhs g," := ",rhs g$ 1287 for each g in sym do write lhs g," := ",rhs g$ 1288 write"end$"$ 1289 shut "drvlist"$ 1290 >>$ 1291 on nat$ 1292 on dfprint$ 1293 off noarg$ 1294 if lisp(save_lists) then << 1295 out "evolist"$ 1296 for each g in sys do write lhs g," := ",rhs g$ 1297 lisp <<terpri()$ 1298 write"</pre> with symmetries"$ terpri()$ 1299 write"<pre>">>$ 1300 for each g in sym do write lhs g," := ",rhs g$ 1301 shut "evolist"$ 1302 out "unolist"$ 1303 lisp <<listprint(reverse cdr fl)$terpri()>>$ 1304 shut "unolist"$ 1305 >>$ 1306 1307 if interactive then <<off batch_mode$ lisp(print_:= 6 )>> 1308 else <<on batch_mode$ lisp(print_:=nil)>>$ 1309 lisp << 1310 use_new_crackout:=t; 1311 %confirm_subst:=t; 1312 collect_sol:=nil; % no collection of solutions, i.e. use of crackout() 1313 subst_1:=11; 1314 pdelimit_1:=400; 1315 max_gc_short:=4; 1316 max_gc_red_len:=4; 1317 1318 % instead of assigning proc_list here in the form 1319 % proc_list_:='( 1320 % to_do 1321 % separation 1322 % subst_level_0 1323 % subst_level_04 1324 % alg_length_reduction 1325 % find_and_use_sub_systems22 1326 % diff_length_reduction 1327 % factorize_to_substitute 1328 % subst_level_35 1329 % decoupling 1330 % factorize_any 1331 % subst_level_4 1332 % alg_solve_single 1333 % stop_batch 1334 % ) 1335 % it is done in the assignment of old_history before the call of 1336 % CRACK, in order to drop the separation step after its initial 1337 % execution 1338 1339 >>$ 1340 1341 if length sys=1 then write"The equation:" else write"The system:"$ 1342 for each g in sys do write g$ 1343 if subl neq {} then << 1344 if length subl=1 then write"Extra condition:" else write"Extra conditions:"$ 1345 off nat; 1346 for each g in subl do write g; 1347 on nat 1348 >>$ 1349 if verbose then << 1350 write"The symmetry:"$ 1351 for each g in sym do write g$ 1352 >>$ 1353 1354 if do_ine_test then << 1355 inelist:=listine(N,sys,sym,fl,inelist,non_lin_test,save_lists); 1356 g:=inelist$ 1357 while g neq {} and first g neq 0 do g:=rest g$ 1358 if g neq {} then return << 1359 out "invalid"$ 1360 write"SYSTEM & SYMMETRY DO NOT SATISFY MINIMAL REQUIREMENTS!"$ 1361 shut "invalid"$ 1362 write"SYSTEM & SYMMETRY DO NOT SATISFY MINIMAL REQUIREMENTS!"$ 1363 >> 1364 >>$ 1365 1366 % Formulating all conditions in stages 1367 1368 if subl neq {} then let subl$% this activates the extra conditions 1369 fbno:=lisp length sysfbl$ 1370 if subl neq {} then power_split_com:=nil$ % because automatic substitutions from 1371 % subs most likely change powers by introducing new factors 1372 if power_split_com then << 1373 lisp (lin_test_const:=gensym())$ 1374 psys:=power_parti(sys)$ msysp:=first psys; psys:=second psys; 1375 psym:=power_parti(sym)$ msymp:=first psym; psym:=second psym; 1376 for totpow:=2:(msysp+msymp) do << % generating all terms with total 1377 % power totpow of all elements of fblist 1378 lisp<<write"Generating all terms of total degree ",totpow$terpri()>>$ 1379 % initialize all partial commutators to zero 1380 newcd:={}$ for g:=1:fbno do newcd:=cons(0,newcd)$ 1381 % generate all pairings giving totpow 1382 for syspow:=1:msysp do 1383 for sympow:=1:msymp do 1384 if (syspow+sympow)=totpow then << 1385 lisp << 1386 write" Pairing now degree ",syspow," terms from the system"$terpri()$ 1387 write" with degree ",sympow," terms from the symmetry"$terpri()>>$ 1388 1389 subl2:=append( 1390 for each g in psys collect <<ls:=lhs g;rs:=part(rhs g,syspow);ls => rs>>, 1391 for each g in psym collect <<ls:=lhs g;rs:=part(rhs g,sympow);ls => rs>>)$ 1392 let subl2; 1393 1394 % generation of partial commutators and adding them to newcd 1395 for g:=1:fbno do % for the g'th evol. equation + related symmetry equation 1396 newcd:=(part(newcd,g):=part(newcd,g) + df(part(rhs part(psys,g),syspow),s) + 1397 (if (lisp t_changes_parity) and 1398 (lisp s_changes_parity) then df(part(rhs part(psym,g),sympow),t) 1399 else (-df(part(rhs part(psym,g),sympow),t)))); 1400 % clearing the t- and s-derivatives of all field variables 1401 clearrules subl2 1402 1403 >>$ 1404 eqnlist:=append(newcd,eqnlist) 1405 >>$ 1406 psys:=psym:=0 1407 >>$ 1408 1409 % Assign t- and s-derivatives for all members of fblist for use in crack_out. 1410 subl2:=append(for each g in sys collect <<ls:=lhs g;rs:=rhs g; ls => rs>>, 1411 for each g in sym collect <<ls:=lhs g;rs:=rhs g; ls => rs>> )$ 1412 let subl2; % This is the ansatz for the system and symmetry to be used 1413 % for the commutator. 1414 if verbose then << 1415 write"sys = ",sys$ 1416 write"sym = ",sym$ 1417 write"df(rhs part(sys,1),s) = ",df(rhs part(sys,1),s)$ 1418 write"df(rhs part(sym,1),t) = ",df(rhs part(sym,1),t)$ 1419 >>$ 1420 1421 if not power_split_com and plain_com then 1422 for g:=1:fbno do 1423 eqnlist:=cons(df(rhs part(sys,g),s) + 1424 (if (lisp t_changes_parity) and 1425 (lisp s_changes_parity) then df(rhs part(sym,g),t) 1426 else (-df(rhs part(sym,g),t))), 1427 eqnlist) 1428 else 1429 if not power_split_com then lisp << 1430 1431 sym:=reval sym; 1432 eqnlist:=cdr eqnlist; % to get rid of 'LIST 1433 1434 % for each symmetry h collect the list of terms of the rhs caddr h 1435 rhssyl:= for each h in cdr sym collect 1436 if pairp caddr h and caaddr h='PLUS then cdaddr h 1437 else {caddr h}; 1438 nw:=0$ 1439 for each w in sysfbl do << 1440 if verbose then write"Integrability conditions are computed for ",w," :"$ terpri()$ 1441 % For each function w of sysfbl generate the full condition 1442 nw:=add1 nw$ % the index of the function w in sysfbl 1443 sycon:=0$ % the condition for this function 1444 np:=0$ % the index of the symmetry of which one term is kept non-zero 1445 1446 % Initialize the s-derivative of all functions to zero 1447 msgbak:=!*msg; !*msg:=nil; 1448 for each h in sysfbl do algebraic <<clear df(lisp h,s);df(lisp h,s):=0>>$ 1449 !*msg:=msgbak; 1450 cn:=0$ % counter of computed partial commutators 1451 lp:=for each p in rhssyl sum length p$ % counter of all partial commutators 1452 1453 psycon:=nil$ % a partial condition for this function 1454 for each p in rhssyl do << % for each rhs of a symmetry 1455 1456 % at first set previous s-deriv again to zero. 1457 msgbak:=!*msg; !*msg:=nil; 1458 if np>0 then <<h:=nth(sysfbl,np)$ 1459 algebraic <<clear df(h,s);df(h,s):=0>> >>$ 1460 !*msg:=msgbak; 1461 1462 % Now we take the next rhs of which one term is not vanishing 1463 np:=add1 np; 1464 h:=nth(sysfbl,np)$ 1465 nq:=0; % the index of the term in p which is kept non-zero 1466 for each q in p do << % compute partial condition and add to psycon 1467 nq:=add1 nq$ 1468 cn:=add1 cn$ 1469 if verbose then algebraic write " ",cn,"(",lp,")/",nw,"(",fbno, 1470 "). Currently non-vanishing: ",q$ 1471 msgbak:=!*msg; !*msg:=nil; 1472 algebraic <<clear df(lisp h,s); 1473 df(lisp h,s):=lisp nth(p,nq)>>$ 1474 !*msg:=msgbak; 1475 psycon:=cons(algebraic << 1476 k:=df(w,t);g:=df(w,s); 1477 df(k,s) + (if (lisp t_changes_parity) and 1478 (lisp s_changes_parity) then df(g,t) 1479 else (- df(g,t))) 1480 >>,psycon) 1481 >> 1482 >>; 1483 algebraic << 1484 sycon:=num lisp cons('PLUS,psycon); 1485 h:=length sycon; 1486 if verbose then 1487 if h=1 then write"The symmetry condition for ",w," has ",h," term." 1488 else write"The symmetry condition for ",w," has ",h," terms." 1489 >>$ 1490 eqnlist:=cons(sycon,eqnlist) 1491 >>$ 1492 eqnlist:=cons('LIST,eqnlist)$ 1493 1494 % cleaning up assignments for s-derivatives 1495 msgbak:=!*msg; !*msg:=nil; 1496 for each h in sysfbl do algebraic clear df(lisp h,s); 1497 !*msg:=msgbak; 1498 1499 >>; 1500 1501 lisp<<msgbak:=!*msg; !*msg:=nil>>; 1502 clearrules subl2$ 1503 clearrules subl$ 1504 lisp (!*msg:=msgbak)$ 1505 1506 subli:=subl$ 1507 symansatz:=subl2$ 1508 1509 sys:=sym:=0$ 1510 1511 if not init then << 1512 lisp(if !*time then <<terpri()$ 1513 write "CPU time so far: ", time() - cpu," ms ", 1514 " GC time so far: ", gctime() - gti," ms"$ 1515 terpri() 1516 >>); 1517 1518 lisp << 1519 1520 % Listing the splitting `variables' 1521 h:='(t s x)$ 1522 eqnlist:=split_simplify({eqnlist,{'LIST},fl,cons('LIST,h),t})$ 1523 if print_ then <<write"Now crack is called."$terpri()>>$ 1524 old_history:='( 1525 1526 cp (!*comma!* 1 3 45 11 8 20 27 30 47 21 38) 1527 % to specify an automatic solving strategy with priorities: 1528 % 1529 % 1 : Hot list of urgent steps 1530 % 3 : Substitution of <=2 terms, only fcts. of less vars., no cases 1531 % 45 : Substitution of <=8 terms in <=1000 terms, 1532 % alg. expressions, no cases 1533 % 11 : Algebraic length reduction of equations 1534 % 53 : Find sub-systems with 2 flin_ functions % taken out, 1535 % % may generate too many equations 1536 % 27 : Length reducing decoupling steps 1537 % 8 : Factorization to subcases leading to substitutions 1538 % 20 : Substitution of <=1000 terms, no cases 1539 % 30 : Do one decoupling step 1540 % 47 : Any factorization 1541 % 21 : Substitution of <=1000 terms 1542 % 34 : Solving an algebraic equation. % taken out as too risky 1543 % 38 : Stop batch mode 1544 s % to display a statistics for the original overdetermined system 1545 % w nil n "w" % to store the system in ASCII for later web-page generation 1546 % This has been disabled for the online demo. 1547 dp % to disable parallellism 1548 % sb "0" % to backup the original system 1549 % This has been disabled for the online demo. 1550 ); 1551 >>; 1552 % on time$ 1553 % eqnlist:=lisp presimplify(eqnlist,fl)$ ############ 1554 crack(eqnlist,inelist,fl,{})$ % ============= CRACK ============== 1555 lisp <<terpri()$ 1556 write if null sol_list then "No" else length sol_list, 1557 if length sol_list=1 then " solution was" else " solutions were", 1558 " found."$ terpri(); 1559 >>$ 1560 >> 1561end$ % of ssym 1562 1563%------- 1564 1565algebraic procedure samelists(a,b)$ 1566if a={} then if b={} then true 1567 else nil 1568 else if b={} then nil 1569 else 1570if first a neq first b then nil 1571 else samelists(rest a,rest b)$ 1572 1573%------- 1574 1575algebraic procedure crack_out(eqns,assigns,freef,ineq,sol_count)$ 1576% to be used by crack called in ssym() 1577% lines marked with %#### have to be commented when ssym is called from script 1578% uses global variablse fblist,t,s,html_out 1579% assumes that t- and s-derivatives for all members of fblist are assigned 1580if lisp use_new_crackout then 1581begin scalar g,h,tm,ratbak,allbak,eqlist,eqcp1,eqcp2,afblist,ff,symlist,frsym, 1582 all,msgbak,sol_count2; 1583 %,natbak 1584 allbak:=lisp !*allfac$ on allfac$ 1585 ratbak:=lisp !*rational$ on rational$ 1586 afblist:=lisp(cons('LIST,sysfbl))$ 1587 let symansatz$ 1588% let ssrules; 1589 html_out:=nil$ % could be made a flag in ssym() 1590 1591 % At first printing the particular solutions correspponding to each 1592 % individual free parameter 1593 1594 eqlist:=for each h in afblist collect sub(assigns,df(h,t))$ 1595 eqlist:=append(eqlist,for each h in afblist collect sub(assigns,df(h,s))); 1596 symlist:=for each h in afblist collect df(h,s); 1597 frsym:={}$ 1598 for each ff in freef do if not freeof(symlist,ff) then frsym:=cons(ff,frsym); 1599 all:={}$ 1600 sol_count2:=0; 1601 1602 write "All symmetries of solution ",sol_count," printed separately:"$ 1603 for each ff in frsym do 1604 if not freeof(eqlist,ff) then << % which really appear %#### 1605 1606 eqcp1:=eqlist; 1607 for each h in frsym do if h neq ff then eqcp1:=sub(h=0,eqcp1); 1608 1609 % For the remaining constant ff we take a value that normalizes the numerical 1610 % coefficient of the first term of the first non-vanishing rhs of this 1611 % symmetry to one. 1612 g:=eqcp1; 1613 while (g neq {}) and freeof(first g,ff) do g:=rest g; 1614 if g neq {} then << 1615 if my_freeof(den first g,x) then h:=den first g else h:=1$ 1616 tm:=first g$ 1617 h:=h/(numcoeff tm)$ 1618 eqcp1:=sub(ff=h,eqcp1) 1619 >>$ 1620 1621 % Has this solution already been found from other constants? 1622 g:=all$ 1623 while (g neq {}) and not samelists(eqcp1,first g) do g:=rest g; 1624 if g={} then << % a new solution 1625 g:={}; 1626 all:=cons(eqcp1,all); 1627 sol_count2:=sol_count2+1; 1628 eqcp2:=eqcp1$ 1629 1630 lisp << %#### 1631 terpri()$ 1632 if html_out then write"</pre>"$ 1633 write "Solution ",sol_count,", symmetry ",sol_count2,". :"$ terpri()$terpri()$ 1634 if length afblist=2 then write"The equation: " 1635 else write"The system: "$ terpri()$ 1636 if html_out then <<write"<pre>"$terpri()>> 1637 1638 >>$ 1639 for each h in afblist do << 1640 g:=df(h,t); 1641 msgbak:=!*msg; !*msg:=nil; 1642 clear df(h,t); 1643 !*msg:=msgbak; 1644 write df(h,t),"=",first eqcp1; 1645 df(h,t):=g; 1646 eqcp1:=rest eqcp1 1647 >>$ 1648 1649 % if (lisp print_more) then 1650 << 1651 lisp << 1652 if html_out then write"</pre>"$ 1653 terpri()$ 1654 write"The symmetry: "$ terpri()$ 1655 if html_out then write"<pre>"$ 1656 >>$ 1657 for each h in afblist do << 1658 g:=df(h,s); 1659 msgbak:=!*msg; !*msg:=nil; 1660 clear df(h,s); 1661 !*msg:=msgbak; 1662 write df(h,s),"=",first eqcp1; 1663 df(h,s):=g; 1664 eqcp1:=rest eqcp1 1665 >> 1666 >>$ 1667 lisp << 1668 if html_out then write"</pre>"$ 1669 terpri()$ 1670 write"And now in machine readable form: "$ terpri()$ 1671 if html_out then write"<p>"$ 1672 if length afblist=2 then write"The equation: " 1673 else write"The system: "$ terpri()$ 1674 if html_out then write"<pre>" 1675 >>$ 1676 off nat$ 1677 for each h in afblist do << 1678 write"df(",h,",t)=",first eqcp2; 1679 eqcp2:=rest eqcp2 1680 >>$ 1681 on nat$ 1682 % if (lisp print_more) then 1683 << 1684 lisp << 1685 if html_out then write"</pre>"$ 1686 terpri()$ 1687 write"The symmetry: "$terpri()$ 1688 if html_out then write"<pre>" 1689 >>$ 1690 off nat$ 1691 for each h in afblist do << 1692 write"df(",h,",s)=",first eqcp2; 1693 eqcp2:=rest eqcp2 1694 >>$ 1695 on nat$ 1696 >> 1697 >> % of `g neq {}' 1698 >>$ % of `for each ff' 1699 1700 % If there are 2 or more symmetries then printing again the solution 1701 % with a linear combination of its symmetries: 1702 1703 if length(sfrsym) > 1 then << 1704 1705 lisp << %#### 1706 terpri()$ 1707 if html_out then write"</pre>"$ 1708 terpri()$ 1709 write "Again solution ",sol_count,", now with ALL its parameters (symmetries):"$ terpri()$ 1710 if length afblist=2 then write"The equation: " 1711 else write"The system: "$ terpri()$ 1712 if html_out then <<write"<pre>"$terpri()>> 1713 >>$ 1714 for each h in afblist do << 1715 g:=df(h,t); 1716 msgbak:=!*msg; !*msg:=nil; 1717 clear df(h,t); 1718 !*msg:=msgbak; 1719 write df(h,t),"=",sub(assigns,g)$ 1720 df(h,t):=g 1721 >>$ 1722 1723 % if (lisp print_more) then 1724 << 1725 lisp << 1726 if html_out then write"</pre>"$ 1727 terpri()$ 1728 write"The symmetry: "$ terpri()$ 1729 if html_out then write"<pre>"$ 1730 >>$ 1731 for each h in afblist do << 1732 g:=df(h,s); 1733 msgbak:=!*msg; !*msg:=nil; 1734 clear df(h,s); 1735 !*msg:=msgbak; 1736 write df(h,s),"=",sub(assigns,g)$ 1737 df(h,s):=g; 1738 >> 1739 >>$ 1740 1741 % And now in machine readable form: 1742 1743 lisp << 1744 if html_out then write"</pre>"$ 1745 terpri()$ 1746 write"And now in machine readable form: "$ terpri()$ 1747 if html_out then write"<p>"$ 1748 if length afblist=2 then write"The equation: " 1749 else write"The system: "$ terpri()$ 1750 if html_out then write"<pre>" 1751 >>$ 1752 off nat$ 1753 for each h in afblist do 1754 write"df(",h,",t)=",sub(assigns,df(h,t)); 1755 on nat$ 1756 1757 % if (lisp print_more) then 1758 << 1759 lisp << 1760 if html_out then write"</pre>"$ 1761 terpri()$ 1762 write"The symmetry: "$terpri()$ 1763 if html_out then write"<pre>" 1764 >>$ 1765 off nat$ 1766 for each h in afblist do 1767 write"df(",h,",s)=",sub(assigns,df(h,s)); 1768 on nat$ 1769 >>$ 1770 >>$ % if more than 1 symmetry 1771 1772 lisp <<write"-----------------------------------"$ terpri()>>$ %#### 1773 clearrules symansatz$ 1774 % ssrules should stay on 1775 if allbak=nil then off allfac$ 1776 if ratbak=nil then off rational$ 1777 1778end$ % of crack_out 1779 1780%------- 1781 1782symbolic procedure bigdpri l$ 1783begin scalar dd,f; 1784 if not !*nat or !*fort then return 'failed; 1785 f:=cadr l; dd:=caddr l; 1786 prin2!* '!D; 1787 ycoord!* := ycoord!*-1; 1788 if ycoord!* < ymin!* then ymin!*:=ycoord!*; 1789 if N_ neq 1 then prin2!* f; 1790 ycoord!* := ycoord!*+1; 1791 if ycoord!* < ymin!* then ymin!*:=ycoord!*; 1792 if dd then 1793 << % prin2!* "("; % to print brackets 1794 maprin dd; 1795 % prin2!* ")"; 1796 >>; 1797 return l; 1798end$ 1799 1800%------- 1801 1802put('!d,'prifn,'bigdpri)$ 1803 1804put('bigdpri,'expt,'inbrackets)$ 1805 1806% put('f,'prifn,'myfpri); % done in ssym() if only one fermion field 1807% put('b,'prifn,'myfpri); % done in ssym() if only one boson field 1808 1809symbolic procedure myfpri(x)$ 1810% if not !*nat or !*fort or null !*no_index then 'failed 1811if not !*nat or !*fort then 'failed 1812 else << prin2!* car x;x>>$ 1813 1814%------- 1815%>>>>>>> The 2001 patch, now in the Reduce development version and in Reduce 3.8 1816 1817% symbolic procedure noncomp u$ 1818% !*ncmp and noncomp1 u$ 1819 1820%------- 1821 1822% symbolic procedure noncomp1 u$ 1823% if null pairp u then nil 1824% else if pairp car u then noncomfp u 1825% else flagp(car u,'noncom) or noncomlistp cdr u$ 1826 1827%------- 1828 1829% symbolic procedure noncomlistp u$ 1830% pairp u and (noncomp1 car u or noncomlistp cdr u)$ 1831 1832%------- 1833 1834% symbolic procedure mchcomb(u,v,op)$ 1835% begin integer n; 1836% n := length u - length v +1; 1837% if n<1 then return nil 1838% else if n=1 then return mchsarg(u,v,op) 1839% else if not smemqlp(frlis!*,v) then return nil; 1840% return for each x in comb(u,n) join 1841% if null ncmp!* then mchsarg((op . x) . setdiff(u,x),v,op) 1842% else (if null y then nil 1843% else mchsarg((op . x) . car y, 1844% if cdr y then reverse v else v,op)) 1845% where y = mchcomb2(x,u,nil,nil,nil) 1846% end$ 1847 1848%------- 1849 1850% symbolic procedure mchcomb2(u,v,w,bool1,bool2)$ 1851% if null u then nconc(reversip w,v) . bool2 1852% else if car u = car v 1853% then if noncomp car u then mchcomb2(cdr u,cdr v,w,t,bool2) 1854% else mchcomb2(cdr u,cdr v,w,bool1,bool2) 1855% else if noncomp car u 1856% then if bool1 then nil 1857% else mchcomb2(u,cdr v,car v . w,t,if bool2 then bool2 else t) 1858% else mchcomb2(u,cdr v,car v . w,bool1,bool2)$ 1859 1860%------- 1861 1862symbolic procedure search_li3(l,carli)$ 1863% Find all sublists of l which start with an element of carli (no nesting) 1864begin scalar b,res$ 1865 1866 if pairp l then << 1867 b:=carli; 1868 while b and null res do if car l = car b then res:=list l 1869 else b:=cdr b 1870 >> else if l and member(l,carli) then res:=list l; 1871 if null res then 1872 while pairp l do << 1873 if b:=search_li3(car l,carli) then res:=union(b,res); 1874 l:=cdr l 1875 >>$ 1876 return res 1877end$ 1878 1879%------- 1880 1881%symbolic procedure itercoeff(sf,splitvar)$ %################# 1882%if not pairp sf then {!*f2a sf} else %################# 1883%if not memq(mvar sf,splitvar) then {!*f2a sf} else %################# 1884%nconc(itercoeff(lc sf,splitvar),itercoeff(red sf,splitvar))$ %################# 1885 1886%------- 1887 1888symbolic procedure add_terms(h)$ 1889% h is a lisp list of terms that each has to get a constant coefficient 1890% and the terms have to be added and returned 1891begin scalar fl,nf$ 1892 h:=for each r in h collect <<nf :=newfct(fname_,nil,nfct_)$ 1893 fl :=cons(nf,fl); 1894 nfct_:=add1 nfct_$ 1895 {'TIMES,nf,r}>>$ 1896 return {if null h then 0 else 1897 if null cdr h then reval car h else reval cons('PLUS,h), 1898 fl} 1899end$ 1900 1901%------- 1902 1903algebraic procedure ssconlaw(N,tw,cw,afwlist,abwlist,pdes,fermi)$ 1904 % N .. the number of superfields theta_i 1905 % tw .. 2 times the differential order of the equation = weight(t) 1906 % cw .. 2 times the differential order of the conservation law 1907 % afwlist .. (algebraic) list of weights of the fermion fields 1908 % f(1), f(2), ... 1909 % abwlist .. (algebraic) list of weights of the boson fields 1910 % b(1), b(2), ... 1911 % The number of elements in the last two lists determines the 1912 % number of fermion and boson fields. 1913 % pdes .. an algebraic list of the pde(s) for which a conservation 1914 % law is to be found, for example 1915 % {df(f(1),t)=df(f(1),x)*b(1)*p9, 1916 % df(b(1),t)=b(1)**3*p3 + d(1,f(1))**2*p2 + 1917 % df(b(1),x)*b(1)*p9 + df(f(1),x)*f(1)*p4}$ 1918 % fermi .. if =t then conserved flow is fermionic, if =nil then bosonic 1919begin 1920 scalar cpu,gti,sw,fwlist,bwlist,allf,allb,g,h,k,l,m,Pt,Px,Pdl,fl,nf,fbno,eqn, 1921 Ptcp,Pxcp,Pdlcp,sol,spezsol,Ql,terms,te,minu,factors,prefac,fa,dli, 1922 dfli,inte,delta,dlirevcp,tr_cl,pdes2,pdes3,fno,bno,syli,deno,allcl, 1923 nontriv,forbid,Pt1,Px1,Pdl1,verbose$ 1924 1925% tr_cl:=t$ 1926 1927 backup_reduce_flags()$ 1928 1929 lisp << 1930 record_hist:=nil;homogen_:=t; 1931 if !*time then <<cpu:=time()$ gti:=gctime()>>$ 1932 input_consistency_test(afwlist,abwlist)$ 1933 N_:=N; % to avoid printing the index of D in crack_out for N=1 1934 fno:=length afwlist-1; 1935 bno:=length abwlist-1; 1936 flist := for g:=1:fno collect {'f,g}; 1937 if flist and length flist=1 and tr_cl=nil then put('f,'prifn,'myfpri) 1938 else put('f,'prifn,nil); 1939 blist := for g:=1:bno collect {'b,g}; 1940 if blist and length blist=1 and tr_cl=nil then put('b,'prifn,'myfpri) 1941 else put('b,'prifn,nil); 1942 fblist:=append(flist,blist) 1943 >>$ 1944 afblist:=lisp(cons('LIST,fblist))$ 1945 for each g in afblist do depend g,x,t$ 1946 1947 lisp << 1948 1949 % creating more substitutions based on pdes 1950 for each g in cdr pdes do << 1951 if is_fermion cadr g then <<fno:=add1 fno; h:={'f,fno}>> 1952 else <<bno:=add1 bno; h:={'b,bno}>>; 1953 algebraic(depend lisp(h),x,t); 1954 pdes2:=cons({'REPLACEBY,cadr g,{'PLUS,caddr g,h}},pdes2); 1955 pdes3:=cons({h,0,{'DIFFERENCE,cadr g,caddr g}},pdes3); 1956 syli:=cons(h,syli) 1957 >>$ 1958 % pde2 is a substitution list introducing a symbol for each equation 1959 % pde3 is an assoc list with an entry for each equation, eg. 1960 % { {symbol_for_eqn, Q(1), df(f(1),t)-...}, 1961 % {symbol_for_eqn, Q(2), df(b(1),t)-...}, ... } 1962 1963 % listing lhs's of equations and substitutions which are no t-derivatives 1964 forbid:=non_t_lhs_dvs(cdr pdes)$ 1965 1966 pdes2:=cons('LIST,pdes2); 1967 if tr_cl then << 1968 algebraic(write"pdes2=",lisp pdes2)$ 1969 write"pdes3="$prettyprint pdes3$ terpri()$ 1970 write"syli=",syli$terpri() 1971 >>$ 1972 1973 % specifying the names of constants in the equations 1974 if fname_list then <<fname_:=car fname_list; 1975 fname_list:=cdr fname_list>> 1976 else fname_:='r$ 1977 1978 fwlist:=cdr reval afwlist$ 1979 bwlist:=cdr reval abwlist$ 1980 1981 % Add to each field its weight as a prefix to create allf and allb 1982 for i:=1:(length flist) do allf:=cons({nth(fwlist,i),nth(flist,i)},allf); 1983 for i:=1:(length blist) do allb:=cons({nth(bwlist,i),nth(blist,i)},allb); 1984 1985 nfct_:=1$ print_:=nil$ 1986 1987 h:=rhs_term_list(N,allf,allb,cw-tw,nil,forbid,verbose); 1988 if fermi then h:= car h else h:=cadr h; 1989 h:=add_terms(h)$ Pt:=car h$ fl:=append(cadr h,fl)$ 1990 eqn:={{'DF,Pt,'t}}$ 1991 1992 if N=0 then << 1993 h:=rhs_term_list(N,allf,allb,cw-2,nil,forbid,verbose); 1994 if fermi then h:=car h else h:=cadr h; 1995 h:=add_terms(h)$ Px:=car h$ fl:=append(cadr h,fl)$ 1996 eqn:=cons({'DF,Px,'x},eqn)$ 1997 Pdl:=nil$ 1998 >> else << 1999 Px:=0; % which is no loss of generality 2000 2001 h:=rhs_term_list(N,allf,allb,cw-1,nil,forbid,verbose); 2002 if !*t_changes_parity then 2003 if fermi then h:= car h else h:=cadr h else 2004 if fermi then h:=cadr h else h:= car h; 2005 Pdl:=nil$ 2006 for g:=N step -1 until 1 do << % step -1 is crucial for backintegration 2007 k:=add_terms(h)$ 2008 Pdl:=cons({'LIST,g,car k},Pdl)$ 2009 fl:=append(cadr k,fl)$ 2010 eqn:=cons({'D,g,car k},eqn)$ 2011 >> 2012 >>$ 2013 2014 eqn:=cons('PLUS,eqn)$ 2015 Pdl:=cons('LIST,Pdl)$ 2016 fl:=cons('LIST,fl) 2017 2018 >>$ 2019 2020 if tr_cl then << 2021 write"Ansatz for Pt: ",Pt$ 2022 write"Ansatz for Px: ",Px$ 2023 write"Ansatz for Pdl: ",Pdl$ 2024 write"list of unknown coefficients: ",fl$ 2025 >>$ 2026 2027 if fl={} then return << 2028 write "No valid ansatz in this case."$ 2029 nil 2030 >>$ 2031 2032 let pdes; % initializing all substitutions based on equations to generate eqn. 2033 eqn:=eqn; % only now after `let pdes' 2034 clearrules pdes; 2035 2036 lisp << % Formulating the conservation law conditions 2037 2038 % Listing the splitting `variables' 2039 h:=search_li3(reval eqn,{'DF,'F,'B,'D,'X,'T}); 2040 if !*time then <<terpri()$ 2041 write "CPU time so far: ", time() - cpu," ms ", 2042 " GC time so far: ", gctime() - gti," ms"$ 2043 terpri()$ 2044 >>$ 2045 2046 % Splitting 2047 k:=setkorder h$ 2048 eqn:=cons('LIST, 2049 for each g in itercoeff(if not pairp eqn then eqn else 2050 if car eqn = '!*sq then reorder numr cadr eqn 2051 else numr simp eqn,h) 2052 collect {'!*SQ,(g . 1),t})$ 2053 setkorder k$ 2054 2055 !!arbint:=0$ 2056 terpri()$ 2057 write length cdr eqn," conditions for ",length cdr fl," unknowns"$ 2058 terpri() 2059 2060 >>$ 2061 2062 % Solving the conservation law conditions 2063 sol:=first solve(eqn,fl)$ 2064 2065 lisp << 2066 % terpri()$ 2067 write if 0=!!arbint then "No" else !!arbint, 2068 if 0 neq !!arbint then " (possibly trivial) " else " ", 2069 if fermi then "fermionic" else "bosonic"," conservation ", 2070 if 1=!!arbint then "law" else "laws"," of weight ",cw, 2071 if 0=!!arbint then "." else ":"$ 2072 terpri()$ 2073 2074 >>$ 2075 2076 nontriv:=0$ 2077 for g:=1:lisp !!arbint do << 2078 spezsol:=sub(arbcomplex(g)=1,sol); 2079 for h:=1:lisp !!arbint do 2080 if h neq g then spezsol:=sub(arbcomplex(h)=0,spezsol); 2081 2082 Ptcp :=sub(spezsol,Pt)$ 2083 Pxcp :=sub(spezsol,Px)$ 2084 Pdlcp:=sub(spezsol,Pdl)$ 2085 2086 % Counting how many P's are nonzero 2087 if Ptcp=0 then k:=0 else k:=1; 2088 if Pxcp neq 0 then k:=k+1; 2089 h:=Pdlcp$ 2090 while h neq {} do << if second first h neq 0 then k:=k+1; 2091 h:=rest h >>$ 2092 2093 % Drop conservation laws with only one non-vanishing component 2094 if k=0 then write "We drop a conservation law with Pt, Px, Pd[i] ", 2095 "all being zero" 2096 else 2097 if k=1 then write "We drop a conservation law with only ", 2098 if Ptcp neq 0 then "Pt" else "one Pd[i]"," nonzero." 2099 else << 2100 2101 Pt1:=Ptcp$ 2102 Px1:=Pxcp$ 2103 Pdl1:=Pdlcp$ 2104 2105 if tr_cl then << 2106 write ">>>>> ",g,". Conservation law: "$ 2107 write"Pt = ",Ptcp$ 2108 write"Px = ",Pxcp$ 2109 h:=Pdlcp$ 2110 while h neq {} do << write"Pd[",first first h,"] = ",second first h; 2111 h:=rest h >> 2112 >>$ 2113 2114 % Determination of the characteristic functions 2115 eqn:=df(Ptcp,t)+df(Pxcp,x)+ 2116 for l:=1:N sum <<m:=first part(Pdlcp,l);D(m,second part(Pdlcp,l))>>$ 2117 if tr_cl then write"rhs = ",eqn$ 2118 2119 if eqn=0 then % write "This conservation law is trivial!" 2120 else lisp << 2121 2122 Ql:=pdes3$ 2123 2124 pdes2:=add_df_rules_to_D_rule pdes2$ 2125 algebraic(let pdes2); 2126 eqn:=reval eqn; 2127 algebraic(clearrules pdes2); 2128 2129 if tr_cl then <<write"rhs in terms of equations = "$ 2130 % by now eqn should not involve any lhs or derivative 2131 % of any lhs of any equation in pdes. 2132 prettyprint eqn; 2133 terpri() 2134 >>$ 2135 2136 % taking care of quotients 2137 deno:=nil; 2138 if pairp eqn and car eqn = 'QUOTIENT then 2139 if is_const caddr eqn then <<deno:=caddr eqn; eqn:=cadr eqn>> 2140 else rederr" Something is wrong!! (0) *****"$ 2141 % converting to a list of terms 2142 if not pairp eqn or car eqn neq 'PLUS then terms:={eqn} else terms:=cdr eqn; 2143 for each te in terms do << 2144 2145 % taking care of the minus sign 2146 minu:=nil; if pairp te and car te = 'MINUS then <<minu:=t; te:=cadr te>>$ 2147 % checking each factor of which only one should have a t-derivative 2148 if not pairp te or car te neq 'TIMES then factors:={te} 2149 else factors:=cdr te; 2150 prefac:=nil; 2151 while factors do << % search until the first factor including an 2152 % element of syli is found 2153 fa:=car factors; factors:=cdr factors; 2154 if freeoflist(fa,syli) then prefac:=cons(fa,prefac) 2155 else << 2156 if tr_cl then << 2157 write" Next term to be partially integrated:";terpri()$ 2158 prettyprint prefac$ terpri()$ write" times "$terpri()$ 2159 prettyprint fa$ terpri()$ write" times "$terpri()$ 2160 prettyprint factors$terpri()$ 2161 write"minu=",minu$ terpri() 2162 >>$ 2163 2164 % In the case of an exponent, a single factor is good enough for us 2165 if pairp fa and car fa='EXPT then << 2166 factors:=cons(if caddr fa=2 then cadr fa 2167 else {'EXPT,cadr fa,sub1 caddr fa},factors); 2168 fa:=cadr fa 2169 >>$ 2170 2171 if null factors then factors:=1 else 2172 if null cdr factors then factors:=car factors 2173 else factors:=cons('TIMES,factors); 2174 2175 if null prefac then prefac:=1 else 2176 if null cdr prefac then prefac:=car prefac 2177 else prefac:=cons('TIMES,reverse prefac); 2178 2179 % Merging prefac and factors, i.e. moving fa to be last factor 2180 if is_fermion fa and is_fermion factors then minu:=not minu; 2181 prefac:={'TIMES,prefac,factors}; 2182 factors:=nil; % i.e. stop of while loop 2183 if minu then prefac:={'MINUS,prefac}; minu:=nil; 2184 2185 dli:=nil; 2186 while pairp fa and car fa='D do <<dli:=cons(cadr fa,dli);fa:=caddr fa>>$ 2187 dli:=reverse dli; 2188 2189 if pairp fa and car fa = 'DF then << 2190 delta:=cadr fa; 2191 dfli:=cddr fa; 2192 h:=nil; 2193 while dfli do << 2194 if fixp car dfli then for k:=2:(car dfli) do h:=cons(car h,h) 2195 else h:=cons(car dfli,h); 2196 dfli:=cdr dfli 2197 >>$ 2198 dfli:=reverse h 2199 >> else <<delta:=fa;dfli:=nil>>; 2200 k:=assoc(delta,Ql); % k has form: {symb for eqn, Q_f(i), df(f(i),t)-rhs} 2201 Ql:=delete(k,Ql); 2202 delta:=caddr k; 2203 if deno then delta:={'QUOTIENT,delta,deno}; 2204 2205 if tr_cl then <<write"delta=",delta," dfli=",dfli$ terpri()>>$ 2206 2207 % At first partial integration of the D-derivatives 2208 while dli do << 2209 h:=car dli; dli:=cdr dli; % h in index of supersym. generator 2210 if is_fermion prefac then prefac:={'MINUS,prefac}; 2211 inte:=delta; 2212 if dfli then inte:=cons('DF,cons(inte,dfli)); 2213 2214 dlirevcp:=reverse dli; 2215 while dlirevcp do << 2216 inte:={'D,car dlirevcp,inte}; 2217 dlirevcp:=cdr dlirevcp 2218 >>$ 2219 algebraic(if tr_cl then write"Pdlcp before = ",Pdlcp)$ 2220 algebraic(Pdlcp:=(part(Pdlcp,h):={first part(Pdlcp,h), 2221 second part(Pdlcp,h) - 2222 lisp({'TIMES,prefac,inte})}))$ 2223 algebraic(if tr_cl then write"Pdlcp after = ",Pdlcp)$ 2224 prefac:=reval {'MINUS,{'D,h,prefac}} 2225 >>$ 2226 % Then partial integration of the x,t-derivatives 2227 while dfli and not zerop prefac do << 2228 if tr_cl then <<write"Partial ",car dfli,"-integration:"$ terpri()>>$ 2229 inte:=delta; 2230 if cdr dfli then inte:=cons('DF,cons(inte,cdr dfli)); 2231 if car dfli='x then << 2232 algebraic(if tr_cl then write"Pxcp before = ",Pxcp)$ 2233 Pxcp:={'DIFFERENCE,Pxcp,{'TIMES,prefac,inte}}; 2234 algebraic(if tr_cl then write"Pxcp after = ",Pxcp)$ 2235 >> else << 2236 algebraic(if tr_cl then write"Ptcp before = ",Ptcp)$ 2237 Ptcp:={'DIFFERENCE,Ptcp,{'TIMES,prefac,inte}}; 2238 algebraic(if tr_cl then write"Ptcp after = ",Ptcp)$ 2239 >>$ 2240 2241 prefac:=reval {'MINUS,{'DF,prefac,car dfli}}; 2242 dfli:=cdr dfli 2243 >>$ 2244 2245 if deno then prefac:={'QUOTIENT,prefac,deno}; 2246 Ql:=cons({car k,{'PLUS,prefac,cadr k},caddr k},Ql); 2247 >> % of not freeoflist(fa,syli) 2248 >> % while factors do 2249 >>; % for each te in terms do 2250 2251 % The integrating factors as well as the conserved densities may now 2252 % involve the symbols, like f(4), b(3) introduced for the equations. 2253 % These symbols have to be replaced by the equations now. 2254 nf:=cons('LIST,for each g in pdes3 collect {'EQUAL,car g,caddr g}) ; 2255 Ql:=for each h in Ql collect {car h,algebraic(sub(nf,lisp cadr h)),caddr h}; 2256 2257 % Is the conservation law trivial? 2258 h:=Ql$ 2259 while h and zerop cadar h do h:=cdr h; 2260 if tr_cl then << 2261 if null h then write"CL is TRIVIAL!" 2262 else write"CL is GENUINE!"$ 2263 terpri()$ 2264 >>$ 2265 2266 if h then algebraic << % CL is not trivial 2267 2268 Ptcp :=sub(nf,Ptcp)$ 2269 Pxcp :=sub(nf,Pxcp)$ 2270 if N>0 then << 2271 % We use D_x(Pxcp) = D(m,D(m,Pxcp)), set Pxcp=0 and add D(m,Pxcp) to 2272 % Pd(m) where m is one index in 1..N (probably 1). 2273 m:=first first Pdlcp; 2274 if tr_cl then write"We add Pxcp=",Pxcp," to Pd[",m,"]."$ 2275 Pdlcp:=cons({m, D(m,Pxcp)+second first Pdlcp}, 2276 rest Pdlcp); 2277 Pxcp:=0$ 2278 Pdlcp:=sub(nf,Pdlcp)$ 2279 >>$ 2280 2281 % Test of correctness 2282 lisp(h:=if length Ql=1 then {'TIMES,cadar Ql,caddar Ql} 2283 else cons('PLUS,for each h in Ql collect 2284 {'TIMES,cadr h,caddr h})); 2285 2286 k:=df(Ptcp,t)+df(Pxcp,x) - (lisp h) + 2287 for g:=1:N sum <<m:=first part(Pdlcp,g);D(m,second part(Pdlcp,g))>>$ 2288 2289 % Making Ql an algebraic list 2290 lisp (Ql:=cons('LIST,for each h in Ql collect {'LIST,caddr h,cadr h})); 2291 2292 % Is conservation law new? --> compare with CL in allcl. 2293 % structure of allcl: {cl1,cl2,...} where each cli is a conservation law 2294 % and has the form {Ql,Pdlcp,Pxcp,Ptcp,no_of_terms_in_all_P} 2295 % h:=allcl$ 2296 % while h neq {} do << 2297 % 2298 % >>$ 2299 %. . . . . . . 2300 2301 if t then << % Print new conservation law 2302 lisp << 2303 nontriv:=add1 nontriv; 2304 terpri()$ 2305 write ">>>>> ",nontriv,". Non-trivial conservation law: "$ 2306 terpri()$terpri()$ 2307 write"At first in the form of a conserved current vanishing mod eqn.s:"$ 2308 terpri()$terpri()$ 2309 >>$ 2310 write"Pt = ",Pt1$ 2311 write"Px = ",Px1$ 2312 % Printing of Pdl1 2313 h:=Pdl1$ 2314 while h neq {} do <<write"Pd[",first first h,"] = ",second first h; 2315 h:=rest h >>$ 2316 2317 write"----- and in machine readable form:"$ 2318 off nat$ 2319 write"Pt = ",Pt1$ 2320 write"Px = ",Px1$ 2321 h:=Pdl1$ 2322 while h neq {} do << write"Pd[",first first h,"] = ",second first h; 2323 h:=rest h >>$ 2324 on nat; 2325 write"Now as conserved current with characteristic functions:"$ 2326 2327 write"Pt = ",Ptcp$ 2328 write"Px = ",Pxcp$ 2329 % Printing of Pdlcp 2330 h:=Pdlcp$ 2331 while h neq {} do <<write"Pd[",first first h,"] = ",second first h; 2332 h:=rest h >>$ 2333 for each h in Ql do write"Q[",first h,"] = ",second h; 2334 write"----- and in machine readable form:"$ 2335 off nat$ 2336 write"Pt = ",Ptcp$ 2337 write"Px = ",Pxcp$ 2338 h:=Pdlcp$ 2339 while h neq {} do << write"Pd[",first first h,"] = ",second first h; 2340 h:=rest h >>$ 2341 for each h in Ql do write"Q[",first h,"] = ",second h; 2342 on nat; 2343 >>; 2344 2345 % We print contradiction message after the values to be able to check 2346 if k neq 0 then << 2347 write"***** ERROR: "$ 2348 write" 0 <> ",k$ 2349 rederr" A test shows a contradiction! *****"$ 2350% write" A test shows a contradiction! *****"$ %@@@@@@@@@@@@@@@@ 2351 >> 2352 2353 >> % there is a non-zero Q after substitution of equations 2354 >> % divergence does not vanish identically 2355 >> % if k>1 then << (more than one non-vanishing component) 2356 >>$ % for g:=1:lisp !!arbint do << 2357 %#1# 2358end$ 2359 2360%------- 2361 2362algebraic procedure ssconl(N,tw,mincw,maxcw,afwlist,abwlist,pdes)$ 2363<<lisp << 2364 terpri()$ 2365 write"##### ssconl: This is the case N",N,"f",length cdr afwlist,"b", 2366 length cdr abwlist,"t",tw,"w"$ 2367 for each g in append(cdr afwlist,cdr abwlist) do write g; 2368 write". #####"$terpri()$ 2369 >>$ 2370 write"The ",if length pdes=1 then "equation" else "system", 2371 " to be investigated:"$ 2372 for each h in pdes do write h$ 2373 pdes:= for each h in pdes collect 2374 lisp {'REPLACEBY,cadr algebraic h,caddr algebraic h}$ 2375 lisp << 2376 terpri()$ 2377 write"Each CL has the form: "$terpri()$ 2378 if N=0 then write" Dt(Pt) + Dx(Px) = Q1*PDE1 + .." 2379 else write" Dt(Pt) + sum_i Di(Pd(i)) = Q1*PDE1 + .."$ 2380 terpri()$ 2381 write"where, e.g. PDE1 is (left hand side) - (right hand side) of PDE 1."$ 2382 terpri() 2383 >>$ 2384 for h:=mincw:maxcw do << 2385 write"NEXT: FERMIONIC CONSERVATION LAWS OF WEIGHT ",h$ 2386 ssconlaw(N,tw,h,afwlist,abwlist,pdes,t)$ 2387 write"NEXT: BOSONIC CONSERVATION LAWS OF WEIGHT ",h$ 2388 ssconlaw(N,tw,h,afwlist,abwlist,pdes,nil) 2389 >> 2390>>$ 2391 2392%------- 2393 2394symbolic procedure lofl(n,mwt)$ 2395% produces lists with n elements, each having a value 0...mwt 2396if n=1 then for i:=0:mwt collect {i,{i}} else 2397for each l in lofl(n-1,mwt) join 2398for j:=0:(mwt-car l) collect {j+car l,cons(j,cadr l)}$ 2399 2400%------- 2401 2402symbolic procedure wgtof(h,ali,wl)$ 2403% subroutine for FindSSWeights below, determines weight of a kernel 2404begin scalar w,wt,m$ 2405 return 2406 if w:=assoc(h,ali) then {cdr w,ali,wl} else 2407 if pairp h and car h='D then << 2408 w:=wgtof(caddr h,ali,wl); 2409 cons({'PLUS,car w,1},cdr w) 2410 >> else 2411 if pairp h and car h='DF then << 2412 w:=wgtof(cadr h,ali,wl); 2413 wt:=car w; 2414 ali:=cadr w; 2415 wl:=caddr w; 2416 h:=cddr h; % h is now a list of differentiations 2417 while h do << 2418 w:=wgtof(car h,ali,wl); 2419 h:=cdr h; 2420 ali:=cadr w; 2421 wl:=caddr w; 2422 if null h or not fixp car h then m:=1 2423 else <<m:=car h; h:=cdr h>>; 2424 wt:={'DIFFERENCE,wt,{'TIMES,car w,m}} 2425 >>$ 2426 {wt,ali,wl} 2427 >> else << 2428 wt:=mkid('w_,if pairp h then gensym() 2429 else h ); 2430 ali:=cons((h . wt), ali); 2431 wl:=cons(wt,wl); 2432 {wt,ali,wl} 2433 >> 2434end$ 2435 2436%------- 2437 2438symbolic operator FindSSWeights$ 2439symbolic procedure FindSSWeights(N,nf,nb,exli,zerowei,verbose)$ 2440% This procedure computes all homogeneities of a list of 2441% expressios exli which can but need not be equations. 2442% zerowei is a list of identifiers which should get zero weight. 2443begin 2444 scalar j,k,s,bsn,w,jmax,h,ex,fwli,bwli,sc,fh,bh,eh,posi,zrop,posh,negh, 2445 bs, % list of all 2^N sequences of 0's and 1's 2446 ali, % association list of (kernel . weight_name) 2447 allali,% list of association lists for all homogeneities 2448 wl, % list of all weight names (needed for solveeval) 2449 sf, % numerator of an equation (minus first terms..) 2450 tf, % (first) term of sf 2451 wtli, % list of conditions resulting from one equation 2452 wt, % weight of one term in one equation 2453 eli, % list of all conditions on weights 2454 ewli, % list of weights of expressions/equations in exli 2455 hlist % list of all homogeneities 2456$ 2457 2458 % Consistency of input: 2459 if not pairp exli or car exli neq 'LIST then return << 2460 write"The input expression is not a list { } ."$terpri()$ 2461 write"Try again."$terpri() 2462 >>$ 2463 2464 % Generation of bs 2465 bs:={nil}; 2466 for j:=1:N do << 2467 for each s in bs do << 2468 bsn:=cons(cons(0,s),bsn); 2469 bsn:=cons(cons(1,s),bsn) 2470 >>; 2471 bs:=bsn; bsn:=nil 2472 >>$ 2473 bs:=for each s in bs collect cons(for each j in s sum j,s)$ 2474 2475 % initialization of ali 2476 ali:={('x . -2)}$ 2477 2478 % identifiers with zero weight 2479 for each h in cdr reval zerowei do 2480 ali:=cons((h . 0), ali); 2481 2482 % the weights of further differentiation variables will be 2483 % added dynamically as they occur 2484 2485 % the weights of all th-variables 2486 for j:=1:N do ali:=cons(({'th,j} . -1), ali); 2487 2488 % the weights of all fermionic fields 2489 for j:=1:nf do << 2490 2491 % adding f(j): 2492 w:=mkid('w_,mkid('f,j)); 2493 wl:=cons(w,wl); 2494 ali:=cons(({'f,j} . w), ali); 2495 2496 % now all coefficients f(j,..), b(j,..) of all terms in the 2497 % th(k) power expansion of f(j): 2498 for each s in bs do 2499 ali:=cons(((if evenp car s then append({'f,j},cdr s) 2500 else append({'b,j},cdr s)) . 2501 (if zerop car s then w 2502 else {'PLUS,w,car s}) ), ali) 2503 >>$ 2504 2505 % the weights of all bosonic fields 2506 for j:=1:nb do << 2507 2508 % adding b(j): 2509 w:=mkid('w_,mkid('b,j)); 2510 wl:=cons(w,wl); 2511 ali:=cons(({'b,j} . w), ali); 2512 2513 % now all coefficients f(j,..), b(j,..) of all terms in the 2514 % th(k) power expansion of b(j): 2515 for each s in bs do 2516 ali:=cons(((if evenp car s then append({'b,j},cdr s) 2517 else append({'f,j},cdr s)) . 2518 (if zerop car s then w 2519 else {'PLUS,w,car s}) ), ali) 2520 >>$ 2521 2522 for each ex in cdr exli do << 2523 2524 sf:=numr (if pairp ex and (car ex = 'EQUAL) then 2525 subtrsq(simp cadr ex,simp caddr ex) else simp ex); 2526 % sf is the numerator of the expression 2527 2528 wtli:=nil; 2529 while sf do << 2530 tf:=first_term_SF sf; sf:=subtrf(sf,tf); % tf is the first term 2531 wt:=nil; 2532 while tf and not domainp tf do << 2533 2534 w:=wgtof(mvar tf,ali,wl); 2535 wt:=cons({'TIMES,ldeg tf,car w},wt); 2536 ali:=cadr w; 2537 wl:=caddr w; 2538 tf:=lc tf 2539 >>; 2540 wt:=if null wt then 0 else 2541 if cdr wt then cons('PLUS,wt) else 2542 car wt$ 2543 wtli:=cons(reval wt,wtli) 2544 >>; 2545 2546 if wtli and cdr wtli then % ie. if there are at least 2 terms 2547 % then formulate a weight condition for each term to have 2548 % the same weight as the first term, eli is the resulting 2549 % list of conditions 2550 for each w in cdr wtli do 2551 eli:=cons(reval {'DIFFERENCE,car wtli,w}, eli)$ 2552 ewli:=cons(car wtli,ewli) 2553 >>; 2554 2555 % solving the system of conditions 2556 !!arbint:=0; 2557 if null eli then 2558 if null wl then s:=nil 2559 else << 2560 s:=nil; 2561 for each w in wl do << 2562 !!arbint:=add1 !!arbint; 2563 s:=cons({'EQUAL,w,{'ARBCOMPLEX,!!arbint}},s) 2564 >>$ 2565 s:=cons('LIST,s) 2566 >> else << 2567 s:=solveeval {cons('LIST,eli),cons('LIST,wl)}; 2568 if s then s:=cdr s; % to get rid of 'LIST 2569 if s then s:=car s; % This was a linear homogeneous problem with 2570 % only one solution which we take (possibly with free parameters). 2571 2572 % If there is just a single unknown then the solution will consist 2573 % only on one {'EQUAL,w_..,..} and not a list --> we make it a list 2574 if pairp s and car s = 'EQUAL then s:={'LIST,s} 2575 2576 >>; 2577 2578 if s then << 2579 2580 fwli:=cons('LIST,for j:=1:nf collect cdr assoc({'f,j},ali)); 2581 bwli:=cons('LIST,for j:=1:nb collect cdr assoc({'b,j},ali)); 2582 ewli:=cons('LIST,ewli); 2583 2584 2585 jmax:=if !!arbint=0 then 1 else !!arbint$ 2586 for j:=1:jmax do << 2587 sc:=s; 2588 for k:=1:!!arbint do 2589 if k neq j then sc:=algebraic(sub(arbcomplex(lisp k)=0,sc)); 2590 2591 % storing all homogeneities 2592 k:=ali; 2593 for each h in cdr sc do 2594 k:=subst(caddr h,reval cadr h,k); 2595 allali:=cons(k,allali); 2596 2597 if j neq 0 then 2598 sc:=algebraic(sub(arbcomplex(lisp j)=1,sc)); 2599 2600 % storing the homogeneity weights of f(1),f(2),.. b(1),b(2),.. 2601 % in sorted lists: 2602 fh:=algebraic(sub(lisp sc,lisp fwli)); 2603 bh:=algebraic(sub(lisp sc,lisp bwli)); 2604 eh:=algebraic(sub(lisp sc,lisp ewli)); 2605 2606 % All homogeneities with strictly positive weights will be 2607 % collected in posh and all others in negh. Then they are 2608 % appended with posh coming first. 2609 posi:=t; zrop:=t; 2610 for each k in cdr fh do << 2611 if k<=0 then posi:=nil; 2612 if k neq 0 then zrop:=nil 2613 >>; 2614 for each k in cdr bh do << 2615 if k<=0 then posi:=nil; 2616 if k neq 0 then zrop:=nil 2617 >>$ 2618 if null zrop then << 2619 if posi then posh:=cons({'LIST,fh,bh,eh},posh) 2620 else negh:=cons({'LIST,fh,bh,eh},negh); 2621 >> 2622 >> 2623 >>; 2624 hlist:=cons('LIST,append(posh,negh)); 2625 2626 % Printing all homogeneities 2627 if verbose then << 2628 if hlist={'LIST} then write"This system is not homogeneous." 2629 else << 2630 if !!arbint=0 then 2631 write"This system has the following homogeneity:" else 2632 write"This system has the following homogeneities:"$ 2633 terpri()$ 2634 for each ali in allali do << 2635 for each w in ali do 2636 algebraic(write"W[",lisp car w,"] = ",lisp reval cdr w)$ 2637 write"================================="$terpri()$ 2638 >>$ 2639 2640 write"The program returns a list of "$ 2641 h:=length hlist - 1; 2642 if h=1 then <<write"one homogeneity"$terpri()$ 2643 write"which is a list of" 2644 >> else <<write h," homogeneities,"$ terpri()$ 2645 write"each homogeneity being a list of"$ 2646 >>$ 2647 terpri()$ 2648 write"a list of all f(1),b(2),.. weights and"$terpri()$ 2649 write"a list of all b(1),b(2),.. weights and"$terpri()$ 2650 write"a list of the weights of equations/expressions"$terpri()$ 2651 write"in the input. "$ 2652 if !!arbint neq 0 then << 2653 write"In this output free parameters"$terpri()$ 2654 write"arbcomplex(i) are replaced by 1:"$ 2655 >>$ 2656 terpri()$ 2657 mathprint hlist$ 2658 >> 2659 >>$ 2660 2661 return hlist 2662 2663end$ 2664 2665%------- 2666 2667algebraic procedure linearize(pdes,nf,nb,tpar,spar)$ 2668% pdes .. a list of equations with a t-derivative on the left hand side 2669% nf .. number of fermion fields f(1) .. f(nf) 2670% nb .. number of boson fields b(1) .. b(nb) 2671% tpar .. whether the variable t is parity changing 2672% spar .. whether the symmetry parameter s of the symmetry that corresponds 2673% to this linearization is parity changing (t) or not (nil) 2674% linearize returns 2675% {list of relations that define the introduced fields like df(f(3),s)=f(6), 2676% list of linearized equations} 2677% Newly introduced fields depend on x,t,s. 2678begin scalar spar_bak,tpar_bak,n,m,p,npdes,rel$ 2679 2680 spar_bak:=lisp s_changes_parity$ 2681 tpar_bak:=lisp t_changes_parity$ 2682 if tpar then on t_changes_parity 2683 else off t_changes_parity$ 2684 2685 % generation of a substitution list f --> F, b --> B for spar=nil 2686 % or f --> B, f --> F got spar=t 2687 n:=length pdes$ 2688 lisp << terpri()$ 2689 write"The following linearization generates ", 2690 if n = 1 then "a condition" else "conditions", 2691 " for the right hand side", 2692 if n = 1 then " " else "s "$terpri()$ 2693 write"of the symmetry: "$ 2694 terpri()$terpri() 2695 >>$ 2696 rel:={}$ 2697 if spar then << 2698 on s_changes_parity$ 2699 for n:=1:nf do <<m:=nb+n;depend b(m),x,t;depend f(n),s; 2700 rel:=cons(df(f(n),s)=b(m),rel);df(f(n),s):=b(m); 2701 lisp write " D_s f(",n,") = b(",m,")">>$ 2702 for n:=1:nb do <<m:=nf+n;depend f(m),x,t;depend b(n),s; 2703 rel:=cons(df(b(n),s)=f(m),rel);df(b(n),s):=f(m); 2704 lisp write " D_s b(",n,") = f(",m,")">> 2705 >> else << 2706 off s_changes_parity$ 2707 for n:=1:nf do <<m:=nf+n;depend b(m),x,t;depend f(n),s; 2708 rel:=cons(df(f(n),s)=f(m),rel);df(f(n),s):=f(m); 2709 lisp write " D_s f(",n,") = f(",m,")">>$ 2710 for n:=1:nb do <<m:=nb+n;depend f(m),x,t;depend b(n),s; 2711 rel:=cons(df(b(n),s)=b(m),rel);df(b(n),s):=b(m); 2712 lisp write " D_s b(",n,") = b(",m,")">> 2713 >>$ 2714 lisp terpri()$ 2715 2716 write"The original system:"$ 2717 for each p in pdes do write p$ 2718 2719 npdes:={}; 2720 for each p in pdes do << 2721 m:=df(lhs p,s)$ 2722 if part(m,0)=minus then n:=t else n:=nil$ 2723 npdes:=sqcons(if tpar and spar then if n then (-m = df(rhs p,s)) 2724 else ( m = -df(rhs p,s)) 2725 else if n then (-m = -df(rhs p,s)) 2726 else ( m = df(rhs p,s)), npdes) 2727 >>$ 2728 2729 npdes:=reverse npdes; 2730 write"The linearized system: "$ 2731 for each p in npdes do write p$ 2732 write"and in machine readable form: "$ off nat$ 2733 for each p in npdes do write p$ 2734 on nat$ 2735 2736 if lisp(spar_bak neq s_changes_parity) then 2737 if spar_bak then on s_changes_parity 2738 else off s_changes_parity$ 2739 2740 if lisp(tpar_bak neq t_changes_parity) then 2741 if tpar_bak then on t_changes_parity 2742 else off t_changes_parity$ 2743 2744 % for n:=1:nf do clear df(f(n),s)$ 2745 % for n:=1:nb do clear df(b(n),s)$ 2746 2747 return {rel,npdes} 2748end$ 2749 2750%------- 2751 2752symbolic operator GenSSPoly$ 2753symbolic procedure GenSSPoly(N,wgtlist,cname,mode)$ 2754begin 2755 scalar linsub,g,h,N,flist,blist,afblist,forbid,non0coeff,fname_bak, 2756 fl,allf,allb,fwlist,bwlist,nf,nb,verbose,afwlist,abwlist,pol, 2757 fonly,bonly,wgt,newf,newb,print_bak$ 2758 2759 wgtlist:=cdr wgtlist; % to drop 'LIST 2760 afwlist:=cadr car wgtlist; 2761 abwlist:=caddr car wgtlist; 2762 wgt :=cadddr car wgtlist; 2763 wgtlist:=cdr wgtlist; 2764 2765 nf:=length afwlist - 1; fwlist:=cdr reval afwlist$ 2766 nb:=length abwlist - 1; bwlist:=cdr reval abwlist$ 2767 2768 % Initialization of flags 2769 mode:=cdr mode; % to drop 'LIST 2770 while mode do << 2771 if car mode = 'lin then << 2772 if not evenp nf then 2773 rederr"The flag `lin' can not be run with odd many fermions"$ 2774 linsub:=for h:=(nf/2+1):nf collect {'f,h}$ 2775 if not evenp nb then 2776 rederr"The flag `lin' can not be run with odd many bosons"$ 2777 linsub:=cons('LIST,append(linsub,for h:=(nb/2+1):nb collect {'b,h})) 2778 >> else 2779 if car mode = 'verbose then verbose:=t else 2780 if car mode = 'fonly then fonly:=t else 2781 if car mode = 'bonly then bonly:=t else 2782 if pairp car mode and cadar mode = 'forbid then forbid :=cddar mode else 2783 if pairp car mode and cadar mode = 'non0coeff then non0coeff:=cddar mode$ 2784 mode:=cdr mode 2785 >>$ 2786 2787 % Initialization of fermionic and bosonic fields + their printing 2788 input_consistency_test(afwlist,abwlist)$ 2789 N_:=N; % to avoid printing the index of D in crack_out for N=1 2790 2791 flist := for g:=1:nf collect {'f,g}; 2792 % if nf=1 then put('f,'prifn,'myfpri) else 2793 % if nf>1 then put('f,'prifn,nil); 2794 % This was commented out to avoid confusion when nf=1 and using th(i) 2795 % expansions and coefficients f(1,1,0,0,1,..). 2796 2797 blist := for g:=1:nb collect {'b,g}; 2798 % if nb=1 then put('b,'prifn,'myfpri) else 2799 % if nb>1 then put('b,'prifn,nil); 2800 2801 % Making f(i) and b(i) x,t-dependent 2802 algebraic << 2803 afblist:=lisp(cons('LIST,append(flist,blist)))$ 2804 for each g in afblist do depend g,x,t$ %<############## t-dependence? 2805 >>$ 2806 2807 % specifying the names of constants in the symmetries 2808 fname_bak:=fname_; fname_:=reval cname; 2809 2810 % Add to each field its weight as a prefix to create allf and allb 2811 for i:=1:nf do allf:=cons({nth(fwlist,i),nth(flist,i)},allf); 2812 for i:=1:nb do allb:=cons({nth(bwlist,i),nth(blist,i)},allb); 2813 2814 % generation of a list of all terms 2815 pol:=rhs_term_list(N,allf,allb,wgt,linsub,forbid,verbose)$ % returns {fprod,bprod} 2816 pol:={'LIST,cons('LIST,car pol),cons('LIST,cadr pol)}$ 2817 2818 %>>>>>> add to avoid all forbid factors also if they are f(i) or b(i) or d(i,..) 2819 2820 % applying the remaining homogeneities 2821 while wgtlist do << 2822 afwlist:=cadr car wgtlist; 2823 abwlist:=caddr car wgtlist; 2824 wgt :=cadddr car wgtlist; 2825 wgtlist:=cdr wgtlist; 2826 2827 algebraic << 2828 % creating a substitution rule 2829 su:={}; 2830 g:=0$for each h in afwlist do <<g:=g+1;su:=cons(f(g)=f(g)*lin_test_const**h,su)>>$ 2831 g:=0$for each h in abwlist do <<g:=g+1;su:=cons(b(g)=b(g)*lin_test_const**h,su)>>$ 2832 pol:=sub(su,pol); 2833 let lrule1; pol:=pol; clearrules lrule1; 2834 let lrule2; pol:=pol; clearrules lrule2; 2835 2836 newf:={}$ 2837 for each h in first pol do << 2838 h:=h/lin_test_const**wgt; 2839 if freeof(h,lin_test_const) then newf:=cons(h,newf) 2840 >>$ 2841 2842 newb:={}$ 2843 for each h in second pol do << 2844 h:=h/lin_test_const**wgt; 2845 if freeof(h,lin_test_const) then newb:=cons(h,newb) 2846 >>$ 2847 2848 pol:={newf,newb} 2849 >> 2850 >>$ 2851 2852 print_bak:=print_; print_:=nil; 2853 newf:=if bonly then 0 2854 else << 2855 h:=for each g in cdadr pol collect <<h :=newfct(fname_,nil,nfct_)$ 2856 fl :=cons(h,fl); 2857 nfct_:=add1 nfct_$ 2858 {'TIMES,h,g}>>$ 2859 if h then if cdr h then cons('PLUS,h) 2860 else car h 2861 else 0 2862 >>$ 2863 newb:=if fonly then 0 2864 else << 2865 h:=for each g in cdaddr pol collect <<h :=newfct(fname_,nil,nfct_)$ 2866 fl :=cons(h,fl); 2867 nfct_:=add1 nfct_$ 2868 {'TIMES,h,g}>>$ 2869 if h then if cdr h then cons('PLUS,h) 2870 else car h 2871 else 0 2872 >>$ 2873 print_:=print_bak; 2874 fname_:=fname_bak; 2875 flin_:=fl$ 2876 2877 return {'LIST,newf,newb,cons('LIST,fl)} 2878end$ 2879 2880%------- 2881 2882algebraic procedure thgen(k)$ 2883% generates a list with 2^N elements, each element generating 2884% a term in the th-power series and has the form: 2885% {product_of_th(i), 2886% number_of_factors_th(i), 2887% {0's and 1's, i.e. a 1 for each th() otherwise zeros} } 2888if k=0 then {{1,0,{}}} 2889 else begin 2890 scalar a,b$ 2891 a:= thgen(k-1)$ 2892 return for each b in a 2893 join {{first b, second b, append (third b,{0})}, 2894 {th(k)*first(b),second(b)+1,append (third b,{1})} } 2895end$ 2896 2897Drule:={D(~i,~a)=>df(a,th(i))+th(i)*df(a,x)}$ % can not be defined within tocoo 2898 2899algebraic procedure tocoo(N,nf,nb,ex)$ 2900% nf: number of fermion fields 2901% nb: number of boson fields 2902% ex: polynomial expression in field form to be transformed into coordinate form 2903begin 2904 scalar thli,subfnb,i,h$ 2905 lisp(N_:=N); 2906 put('f,'prifn,nil); 2907 put('b,'prifn,nil); 2908 thli:=thgen(N)$ 2909 2910 % Generate the coordinate form of f(i) 2911 subfnb:={}; 2912 for i:=1:nf do 2913 subfnb:=cons(f(i)=for each h in thli sum 2914 if evenp second h then lisp(cons('f,cdr algebraic cons(i,third h)))*first h 2915 else lisp(cons('b,cdr algebraic cons(i,third h)))*first h, 2916 subfnb); 2917 2918 % Generate the coordinate form of b(i) 2919 for i:=1:nb do 2920 subfnb:=cons(b(i)=for each h in thli sum 2921 if evenp second h then lisp(cons('b,cdr algebraic cons(i,third h)))*first h 2922 else lisp(cons('f,cdr algebraic cons(i,third h)))*first h, 2923 subfnb); 2924 2925 % Do the substitution 2926 ex:=sub(subfnb,ex)$ 2927 % Replace D(i,..) which has to be done AFTER the above substitution 2928 let Drule; 2929 ex:=ex; 2930 clearrules Drule; 2931 2932 return ex 2933 2934end$ 2935 2936%------- 2937 2938algebraic procedure tofield(N,nf,nb,ex,zerowei)$ 2939% tw the total weight of the input expression ex which is 2940% expected to be either purely fermionic or purely bosonic. 2941 2942% In a call to GenSSPoly() an ansatz (ssp) in field form with the proper 2943% homogeneity weights is generated, then converted to coordinate form (ans) 2944% in a call to tocoo() and compared to the input in coordinate form (ex). 2945% The difference is splitted through a call to split_simp() and the 2946% resulting system of conditions for the undetermined coefficients of 2947% the ansatz has to vanish. 2948 2949if lisp (pairp algebraic ex and car ex = 'LIST) then 2950write"The input should be a polynomial, not an equations!" else 2951begin scalar ssp,ans,fieldans,rlist,indepv,h,eqnli,solu,wflist,wblist,tw, 2952 gs,k,ft,sb,notnum,notpos$ 2953 2954 w:=FindSSWeights(N,nf,nb,{ex},zerowei,t)$ % e.g. w={{{1,3},{},{5}}} 2955 if w={} then return lisp << 2956 write"The input expression does not have a homogeneity"$ terpri()$ 2957 write"which could be used to compute a field form."$ terpri() 2958 >>$ 2959 2960 w:=first w$ % We just use the fist homogeneity because if there are 2961 % several and there is one with strictly positive weights 2962 % then this comes first. 2963 2964 wflist:=first w; 2965 wblist:=second w; 2966 tw:=first third w; 2967 2968 for each h in append(wflist,wblist) do 2969 if not fixp h then notnum:=t else 2970 if h <= 0 then notpos:=t$ 2971 2972 if notnum then lisp << 2973 write"##### Strange error: returned weights returned"$terpri()$ 2974 write" from FindSSWeights() are not integers! #####"$terpri() 2975 >> else 2976 if notpos then lisp << 2977 gs:=gensym()$ 2978 k:=setkorder list gs$ 2979 ft:=append(search_li2(algebraic ex,'f), 2980 search_li2(algebraic ex,'b) )$ % all f(..), b(..) functions 2981 for each h in ft do sb:=cons((h . {'TIMES,h,gs}), sb)$ 2982 h:=numr simp {'!*sq,subsq(simp algebraic ex,sb),nil}$ 2983 algebraic(max_deg):=if mvar h neq gs then 0 else 2984 ldeg numr simp {'!*sq,subsq(simp algebraic ex,sb),nil}$ 2985 setkorder k 2986 >>$ 2987 2988 % generating an ansatz in field form (ssp) 2989 % and convert it back to coordinate form (ans) 2990 % Next line works so far only if not compilied. 2991 if is_fermionic(lisp{'!*SQ,((first_term_SF numr simp algebraic ex) . 1),t}) then << 2992 % if (lisp{'!*SQ,((first_term_SF numr simp algebraic ex) . 1),t})**2=0 then << 2993 ssp:=gensspoly(N,{{wflist,wblist,tw}},r,{fonly}); 2994 fieldans:=first(ssp); 2995 >> else << 2996 ssp:=gensspoly(N,{{wflist,wblist,tw}},r,{bonly}); 2997 fieldans:=second(ssp); 2998 >>; 2999 3000 ans:=tocoo(N,nf,nb,fieldans)$ 3001 3002 %spliting the equation 3003 rlist:=third(ssp); % a list of undetermined coefficients r1,r2, ... 3004 indepv:=append({x,t},for h:=1:N collect th(h))$ % all independent variables 3005 lisp(print_:=nil)$ % to suppress printing about the splitting 3006 eqnli:=split_simp({ans-ex},{},rlist,indepv,t)$ 3007 3008 %applying crack to find solutions 3009 solu:=crack(eqnli,{},rlist,{})$ 3010 3011 return if solu={} then lisp << 3012 write"The expression to be transformed has a homogeneity"$terpri()$ 3013 write"but a field form does not exist."$terpri() 3014 >> else <<solu:=second first (solu); 3015 sub(solu,fieldans)>> 3016 3017end$ 3018 3019%------- 3020 3021symbolic procedure restore_interactive_prompt$ 3022 begin scalar !*redefmsg,!*usermode$ 3023 copyd('update_prompt,'restore_update_prompt) 3024 end$ 3025 3026%------- 3027 3028symbolic procedure change_prompt$ 3029 begin scalar !*usermode$ 3030 if null promptstring!* then promptstring!* := ""; 3031 setpchar promptstring!*; 3032 promptexp!* := promptstring!* 3033 end$ 3034 3035%------- 3036 3037symbolic procedure change_prompt_to u$ 3038 begin scalar oldprompt,!*redefmsg,!*usermode$ 3039 oldprompt := promptstring!*; 3040 promptstring!* := u; 3041 copyd('restore_update_prompt,'update_prompt); 3042 copyd('update_prompt,'change_prompt); 3043 update_prompt(); 3044 restore_interactive_prompt()$ 3045 return oldprompt 3046 end$ 3047 3048%------- 3049 3050symbolic procedure cnt_l_$ 3051if lines_written geq 10 then << 3052 change_prompt_to ""$ 3053 write" Press Enter to continue "$ 3054 restore_interactive_prompt()$ 3055 rds nil; wrs nil$ % Switch I/O to terminal 3056 readch()$ 3057 if ifl!* then rds cadr ifl!*$ % Resets I/O streams 3058 if ofl!* then wrs cdr ofl!*$ 3059 lines_written:=0 3060>>$ 3061 3062%------- 3063 3064symbolic procedure wl1(l)$ 3065<<write l$ terpri()$ cnt_l_()$ lines_written:=lines_written+1>>$ 3066symbolic procedure wl2(l)$ 3067<<write l$ terpri()$ terpri()$ cnt_l_()$ lines_written:=lines_written+2>>$ 3068 3069%------- 3070 3071symbolic procedure sshelp1$ 3072<<wl2"Purpose: "$ 3073 3074 wl1"to determine evolutionary supersymmetric PDEs with higher"$ 3075 wl1"symmetries, to compute symmetries and conservation laws for"$ 3076 wl1"given supersymmetric PDEs, or to do any other algebra or"$ 3077 wl2"differentiations of polynomial supersymmetric expressions."$ 3078>>$ 3079 3080%------- 3081 3082symbolic procedure sshelp2$ 3083<<wl2"Notation: "$ 3084 3085 wl1"Fields f(1),f(2),.. are treated as fermionic, and b(1),.. are"$ 3086 wl1"treated as bosonic. The independent variables are t, x and in"$ 3087 wl1"symmetry computations the symmetry `time' variable is s. "$ 3088 wl2"Further fermionic variables can be defined, for example, theta."$ 3089 3090 wl1"The super derivatives D_i which satisfy (D_i)**2 = d/dx "$ 3091 wl2"are implemented as d(i, .. ) . "$ 3092 3093 wl1"When refering to homogeneity weights, all weights are scaled"$ 3094 wl1"such that the weight of d(i, .. ) is 1 and thus the weight"$ 3095 wl2"of df( .. , x) is two, i.e. twice the usual value. "$ 3096>>$ 3097 3098%------- 3099 3100symbolic procedure sshelp3$ 3101<<wl2"Initializations: "$ 3102 3103 wl1"Before starting super-computations the number N of superfields theta_i"$ 3104 wl1"needs to be known. Also the number nb of boson fields and the nf of"$ 3105 wl1"fermion fields is needed to have compact printing on the screen (to"$ 3106 wl1"avoid printing the index 1 if nf=1 or nb=1). The input of N,nf,nb is"$ 3107 wl1"done either through a call of the procedure ssym (to compute"$ 3108 wl1"symmetries) or of the procedure ssconl (to compute conservation laws)"$ 3109 wl1"or with a call of procedure ssini. This help item describes ssini."$ 3110 wl2"It's format is:"$ 3111 3112 wl2" ssini(N,nf,nb)$"$ 3113 3114 wl2"where"$ 3115 3116 wl1" N .. the number of superfields theta_i"$ 3117 wl1" nf .. number of fermion fields f(1) .. f(nf)"$ 3118 wl2" nb .. number of boson fields b(1) .. b(nb)"$ 3119 3120 wl2"-----------------"$ 3121 3122 wl2"Example: For N=1 and 2 fermion and 2 boson fields we do"$ 3123 3124 wl2"ssini(1,2,2)$ "$ 3125 3126 wl2"and verify easily that fermions anticommute: "$ 3127 3128 wl1"f(1)**2; "$ 3129 wl2"f(1)*f(2)+f(2)*f(1);"$ 3130 3131 wl2"that bosons commute:"$ 3132 3133 wl1"b(1)**3;"$ 3134 wl2"b(1)*b(2)-b(2)*b(1);"$ 3135 3136 wl2"that x,t-differentiations by default do not change parity"$ 3137 3138 wl1"df(b(1),t)**2;"$ 3139 wl2"df(b(1),x)**2;"$ 3140 3141 wl2"unless one explicitly defines d/dt to be parity changing through"$ 3142 3143 wl2"on t_changes_parity$"$ 3144 3145 wl2"demonstrated by"$ 3146 3147 wl2"df(b(1),t)**2;"$ 3148 3149 wl1"which is also possible to define for the symmetry variable s "$ 3150 wl1"(through on s_changes_parity$ , see below)"$ 3151 wl2"but not for x. For the remainder of this demo we change back"$ 3152 3153 wl1"off t_changes_parity$"$ 3154 wl2"df(b(1),t)**2;"$ 3155 3156 wl1"Differentiations d/dt, d/dx, D_i follow the supersymmetric"$ 3157 wl2"product rule:"$ 3158 3159 wl1"df(b(1)*f(2),x);"$ 3160 wl1"d(1,b(1)*f(2));"$ 3161 wl1"df(f(1)*f(2),x);"$ 3162 wl2"d(1,f(1)*f(2));"$ 3163 3164 wl2"furthermore D_i changes parity:"$ 3165 3166 wl1"d(1,f(1))**2;"$ 3167 wl2"d(1,b(1))**2;"$ 3168 3169 wl2"and satisfies D^2 = d/dx:"$ 3170 3171 wl2"d(1,d(1,b(1)));"$ 3172 3173 wl1"To compactify output, in the special case N=1 the index of D_i"$ 3174 wl1"is not printed as seen above. Similarly, indices of f(1), b(1)"$ 3175 wl2"are not printed if we have only one fermion and/or one boson field:"$ 3176 3177 wl1"ssini(2,1,2)$"$ 3178 wl2"f(1)*b(1)*d(1,b(2));"$ 3179 3180 wl1"Any initializations made by ssini are overwritten by a call to the"$ 3181 wl2"procedures ssym or ssconl."$ 3182 3183>>$ 3184 3185%------- 3186 3187symbolic procedure sshelp4$ 3188<<wl2"Coefficients: "$ 3189 3190 wl1"coeffn is a standard REDUCE command although its implementation"$ 3191 wl1"(at least in REDUCE distributions up to version 3.8) does not work"$ 3192 wl1"with non-commuting variables. By loading SsTools some of the REDUCE"$ 3193 wl1"routines are re-defined so that coeffn can be used for supersymmetric"$ 3194 wl2"expressions. The call is unchanged from the standard REDUCE command:"$ 3195 3196 wl2"coeffn(ex,h,n)$"$ 3197 3198 wl1"which computes the coefficient of h**n in the expression ex which"$ 3199 wl1"is supposed to be polynomial in h. The returned result is the"$ 3200 wl1"coefficient AFTER h**n has been moved to be the first factor in"$ 3201 wl2"each term of ex. This matters if h is fermionic and thus n=1."$ 3202 3203 wl2"Example:"$ 3204 3205 wl1"SSIni(2,0,2)$"$ 3206 wl1"ex:=D(2,b(1))*D(1,df(b(2),x,2));"$ 3207 wl2"coeffn(ex,D(1,df(b(2),x,2)),1);"$ 3208 3209 wl1"Because D_2 b(1) and D_1 D_x D_x b(2) are fermionic a minus sign is"$ 3210 wl1"introduced when making D_1 D_x D_x b(2) the first factor in the"$ 3211 wl2"product. This explains the minus in the result of the coeffn call."$ 3212 3213 wl1"NOTE: When using SsTools the REDUCE command noncom has no effect,"$ 3214 wl1"i.e. it is not possible to define additional noncommuting quantities"$ 3215 wl2"besides fermionic ones which are introduced by the command fermionic."$ 3216 3217>>$ 3218 3219%------- 3220 3221symbolic procedure sshelp5$ 3222<<% special use comments: 3223 % in drvlist$ % loads the ansatz of the previous ssym run 3224 % use_new_crackout:=t$ % needed to use the ssym specific crack_out() 3225 3226 wl2"Computing symmetries:"$ 3227 3228 wl1"To compute PDEs together with higher symmetries, or to"$ 3229 wl2"compute higher symmetries for a given PDE(-system) the call is:"$ 3230 3231 wl2" ssym(N,tw,sw,afwlist,abwlist,eqnlist,fl,inelist,mode)$"$ 3232 3233 wl2"where"$ 3234 3235 wl1" N .. the number of superfields theta_i"$ 3236 wl1" tw .. 2 times the differential order of the equation = weight(t)"$ 3237 wl1" sw .. 2 times the differential order of the symmetry = weight(s)"$ 3238 wl1" afwlist .. (algebraic mode) list of weights of the fermion fields"$ 3239 wl1" f(1), f(2), ... . The number of elements of this list"$ 3240 wl1" determines the number of fermion fields. "$ 3241 wl1" abwlist .. (algebraic mode) list of weights of the boson fields"$ 3242 wl1" b(1), b(2), ... . The number of elements of this list"$ 3243 wl1" determines the number of boson fields. "$ 3244 wl1" eqnlist .. - in the nonlinear case a list of extra conditions on the"$ 3245 wl1" undetermined coefficients where conditions a3=.. are executed"$ 3246 wl1" instantly. These and any expressions are added to equations"$ 3247 wl1" when calling crack"$ 3248 wl1" - in the linear case the system in form of replacements => and"$ 3249 wl1" the linearized system in form of equations = ."$ 3250 wl1" fl .. extra unknowns in eqnlist to be determined"$ 3251 wl1" inelist .. a list, each element of it is a list with at least one of"$ 3252 wl1" its elements being non-zero"$ 3253 wl1" mode .. list of flags: "$ 3254 wl1" init: only initialization of global data"$ 3255 wl1" plain_com : direct computation of the commutator (for tests)"$ 3256 wl1" power_split_com : alternatice power splitting of commutator,"$ 3257 wl1" (not if substitutions '=>' are present,"$ 3258 wl1" see below)"$ 3259 wl1" zerocoeff: all coefficients = 0 which do not appear in inelist"$ 3260 wl1" filter: extra homogeneity weights as given in the list"$ 3261 wl1" hom_wei have to be satisfied by the symmetry"$ 3262 wl1" lin: find symmetry that is linear homogeneous in all fields "$ 3263 wl1" with (index)>(maxindex/2), i.e. if lin then the"$ 3264 wl1" number of fermions and boson fields must both be even"$ 3265 wl1" tpar: t/nil, whether time variable t changes parity"$ 3266 wl2" spar: t/nil, whether symmetry variable s changes parity"$ 3267 3268 wl2"-----------------"$ 3269 3270 wl1"The question whether ssym is to be used to compute a PDE(-system) with"$ 3271 wl1"symmetries, or to compute symmetries for a `given' PDE(-system) is"$ 3272 wl1"decided based on the form of the input eqnlist. "$ 3273 wl1"If symmetries of a specific weight for a given PDE(-system) are to"$ 3274 wl1"be computed then eqnlist has to have the form"$ 3275 wl1"{df(f(1),t)=..., df(b(1),t)=...} with as many t-derivatives"$ 3276 wl1"of fermion fields as there are elements in afwlist and as many"$ 3277 wl2"t-derivatives of boson fields as there are elements in abwlist."$ 3278 3279 wl1"If the right hand sides of the PDE(s) in eqnlist involve any "$ 3280 wl1"constants which are to be computed then these constants have to "$ 3281 wl1"be listed in the input list fl, like {p1,p2,...}. If these constant "$ 3282 wl1"coefficients are not given in the list fl then they are treated as "$ 3283 wl1"independent parameters and symmetries are determined which are "$ 3284 wl1"symmetries for generic values of these constant coefficients."$ 3285 wl1"Also, when computing symmetries for a given PDE(-system) then "$ 3286 wl1"dependencies of all the fields on t,x has to be declared"$ 3287 wl1"beforehand in order to be able to input derivatives, like df(f(1),t) "$ 3288 wl2"which otherwise would be automatically evaluated to zero."$ 3289 3290 wl1"If the input list eqnlist does not start with {df(..)=..} then"$ 3291 wl1"an ansatz for the PDE(-system) is generated and this system and"$ 3292 wl1"its symmetries are computed. In this case eqnlist can contain"$ 3293 wl1"substitutions, like p2=0, or expressions which are to be set to"$ 3294 wl1"zero, like p3*r4+p2*r3. Typically one would run the program first "$ 3295 wl1"to see which ansatz for the PDE(-system) and the symmetry is "$ 3296 wl1"generated and then start ssym again with the intended extra"$ 3297 wl2"conditions."$ 3298 3299 wl1"If a boson weight is non-positive or a fermion weight is negative"$ 3300 wl1"then the global (algebraic mode) variable max_deg must have a positive"$ 3301 wl2"integer value which is the highest degree of such a field in any ansatz."$ 3302 3303 wl1"When determining symmetries (either for a given system of equations,"$ 3304 wl1"or when determining equations and symmetries in one run) one has two"$ 3305 wl1"cases: whether the symmetry variable s flips parity, or not, in other"$ 3306 wl1"words, whether s itself is fermionic or bosonic. This has nothing to do"$ 3307 wl1"with the question whether s has an odd or even weight."$ 3308 wl1"Similarly, when determining equations and symmetries in one computation"$ 3309 wl1"one can consider the two cases whether the `time' derivative changes"$ 3310 wl1"parity or not. A parity changing t is specified through the flag `tpar'"$ 3311 wl2"and a parity changing s through the flag `spar'."$ 3312 3313 wl2"-----------------"$ 3314 3315 wl2"Example for determining PDEs + symmetries:"$ 3316 3317 wl2"ssym(1,4,5,{2},{2},{},{},{},{})$ "$ 3318 3319 wl2"Example for determining the symmetries of a given PDE-system:"$ 3320 3321 wl1"ssym(1,4,5,{2},{2},"$ 3322 wl1" {df(f(1),t)=df(f(1),x)*b(1)*p9,"$ 3323 wl1" df(b(1),t)=df(b(1),x)*b(1)*p9 + df(f(1),x)*f(1)*p4},"$ 3324 wl2" {},{},{})$ "$ 3325 3326 wl2"-----------------"$ 3327 3328 wl1"Example for determining symmetries of given PDEs in the"$ 3329 wl1"presence of extra substitution rules which are given"$ 3330 wl1"using `=>' instead of `=' for the equations for which"$ 3331 wl1"matching symmetries are to be found:"$ 3332 3333 wl1"lisp put('f,'prifn,nil)$ % from now on more than one fermion "$ 3334 wl1"lisp put('b,'prifn,nil)$ % from now on more than one boson "$ 3335 wl1"ssym(1,1,4,{1,1},{1,3,1,3}, "$ 3336 wl1" {df(f(1),t)=>-2*f(1)*b(1)*p1, "$ 3337 wl1" df(b(1),t)=>b(1)**2*p1+d(1,f(1))*p2, "$ 3338 wl1" d(1,b(2)) =>-df(b(1),x)*f(1), "$ 3339 wl1" df(b(2),t)=>-2*d(1,b(1))*f(1)*b(1)*p1+d(1,df(b(1),t))*f(1)"$ 3340 wl1" -d(1,f(1))**2*p2/2-d(1,f(1))*b(1)**2*p1"$ 3341 wl1" +d(1,f(1))*df(b(1),t),"$ 3342 wl1" df(f(2),t)=-2*f(2)*b(1)*p1-2*f(1)*b(3)*p1, "$ 3343 wl1" df(b(3),t)=2*b(1)*b(3)*p1+d(1,f(2))*p2, "$ 3344 wl1" d(1,b(4))=>-df(b(3),x)*f(1)-df(b(1),x)*f(2), "$ 3345 wl1" df(b(4),t)=>-2*d(1,b(3))*f(1)*b(1)*p1-2*d(1,b(1))*f(2)*b(1)*p1"$ 3346 wl1" -2*d(1,b(1))*f(1)*b(3)*p1+d(1,df(b(3),t))*f(1) "$ 3347 wl1" +d(1,df(b(1),t))*f(2)-d(1,f(1))*d(1,f(2))*p2 "$ 3348 wl1" -d(1,f(2))*b(1)**2*p1-2*d(1,f(1))*b(3)*b(1)*p1"$ 3349 wl1" +d(1,f(2))*df(b(1),t)+d(1,f(1))*df(b(3),t)},"$ 3350 wl2" {},{},{}); "$ 3351 3352 wl2"-----------------"$ 3353 3354 wl1"The same example again but now with the flag `lin' to generate a"$ 3355 wl1"symmetry that is linear homogeneous in all fields with an index"$ 3356 wl1"greater than half the maximum index for that type of field. When "$ 3357 wl1"setting the flag 'lin' the number of fermions and boson fields must"$ 3358 wl1"both be even. In the following example that means, the symmetry will"$ 3359 wl2"be linear and homogeneous in f(2),b(3),b(4)."$ 3360 3361 wl1"ssym(1,1,4,{1,1},{1,3,1,3}, "$ 3362 wl1" {df(f(1),t)=>-2*f(1)*b(1)*p1, "$ 3363 wl1" df(b(1),t)=>b(1)**2*p1+d(1,f(1))*p2, "$ 3364 wl1" d(1,b(2)) =>-df(b(1),x)*f(1), "$ 3365 wl1" df(b(2),t)=>-2*d(1,b(1))*f(1)*b(1)*p1+d(1,df(b(1),t))*f(1)"$ 3366 wl1" -d(1,f(1))**2*p2/2-d(1,f(1))*b(1)**2*p1"$ 3367 wl1" +d(1,f(1))*df(b(1),t),"$ 3368 wl1" df(f(2),t)=-2*f(2)*b(1)*p1-2*f(1)*b(3)*p1, "$ 3369 wl1" df(b(3),t)=2*b(1)*b(3)*p1+d(1,f(2))*p2, "$ 3370 wl1" d(1,b(4)) =>-df(b(3),x)*f(1)-df(b(1),x)*f(2), "$ 3371 wl1" df(b(4),t)=>-2*d(1,b(3))*f(1)*b(1)*p1-2*d(1,b(1))*f(2)*b(1)*p1"$ 3372 wl1" -2*d(1,b(1))*f(1)*b(3)*p1+d(1,df(b(3),t))*f(1) "$ 3373 wl1" +d(1,df(b(1),t))*f(2)-d(1,f(1))*d(1,f(2))*p2"$ 3374 wl1" -d(1,f(2))*b(1)**2*p1-2*d(1,f(1))*b(3)*b(1)*p1"$ 3375 wl1" +d(1,f(2))*df(b(1),t)+d(1,f(1))*df(b(3),t)},"$ 3376 wl2" {},{},{lin}); "$ 3377 3378 wl2"-----------------"$ 3379 3380 wl1"If one wants to find for a given system of equations df(..,t)=.. a "$ 3381 wl1"symmetry which satisfies in addition to a symmetry weight sw, a list "$ 3382 wl1"of fermion weights afwlist and list of boson weights abwlist other "$ 3383 wl1"homogeneities, then in the last parameter of the ssym call one adds the "$ 3384 wl1"flag `filter' and has a global variable hom_wei which is a list of "$ 3385 wl1"lists {sw,afwlist,abwlist} each being an additional set of homogeneity"$ 3386 wl2"weights."$ 3387 3388 wl2"Example: "$ 3389 3390 wl2"hom_wei:={{10,{3,3},{2,2}}}$"$ 3391 3392 wl1"ssym(1,1,6,{1,1},{1,1},"$ 3393 wl1" {df(f(1),t)=>-2*f(1)*b(1)*p1,"$ 3394 wl1" df(b(1),t)=>b(1)*b(1)*p1+d(1,f(1))*p2,"$ 3395 wl1" df(f(2),t)= - 2*f(2)*b(1)*p1 - 2*f(1)*b(2)*p1,"$ 3396 wl1" df(b(2),t)=2*b(2)*b(1)*p1 + d(1,f(2))*p2},"$ 3397 wl2" {},{},{lin,filter});"$ 3398 3399 wl2"-----------------"$ 3400 3401 wl1"After a run of ssym() the substitution rules, like df(f(1),t)=> .. "$ 3402 wl1"are not active. If one wanted to check symmetries interactively "$ 3403 wl1"one would have to activate these rules which are stored under the"$ 3404 wl1"name `subli' by: let subli$"$ 3405 wl1"and if one would want to de-activate them then do: clearrules subli$ "$ 3406 wl2"For example, to check symmetries after the above run one could do "$ 3407 3408 wl2" let subli$ "$ 3409 3410 wl1" df(b(3),t):=2*b(3)*b(1)*p1 + d(1,f(2))*p2$ "$ 3411 wl2" df(f(2),t):=- 2*f(2)*b(1)*p1 - 2*f(1)*b(3)*p1$ "$ 3412 3413 wl1" b3t:=df(b(3),t)$ "$ 3414 wl2" f2t:=df(f(2),t)$ "$ 3415 3416 wl1" df(b(3),s):=(1/2*d(1,b(3))*f(1)*b(1)**2*p1*p2 + 4/11*d(1,b(1))*f(2) "$ 3417 wl1" *b(1)**2*p1*p2 + d(1,b(1))*f(1)*b(3)*b(1)*p1*p2 + 7/22*d(1,f(1))**2 "$ 3418 wl1" *b(3)*p2**2 + 7/11*d(1,f(1))*b(3)*b(1)**2*p1*p2 + 1/4*d(1,f(1))* "$ 3419 wl1" d(1,b(3))*f(1)*p2**2 + 5/44*d(1,f(1))*d(1,b(1))*f(2)*p2**2 + 1/44* "$ 3420 wl1" df(b(1),x)*f(2)*f(1)*p2**2 + 1/22*df(f(1),x)*f(2)*b(1)*p2**2 "$ 3421 wl1" + 5/22*df(f(1),x)*f(1)*b(3)*p2**2)/(p1*p2)$ "$ 3422 wl1" df(f(2),s):=(1/4*d(1,f(2))*d(1,f(1))*f(1)*p2**2 + 1/2*d(1,f(2)) "$ 3423 wl1" *f(1)*b(1)**2*p1*p2 + 3/44*d(1,f(1))**2*f(2)*p2**2 + 3/22*d(1,f(1)) "$ 3424 wl2" *f(2)*b(1)**2*p1*p2 + 1/44*df(f(1),x)*f(2)*f(1)*p2**2)/(p1*p2)$ "$ 3425 3426 wl1" b3s:=df(b(3),s)$ "$ 3427 wl2" f2s:=df(f(2),s)$ "$ 3428 3429 wl1" write ""ZERO = "",df(b3t,s)-df(b3s,t)$ "$ 3430 wl2" write ""ZERO = "",df(f2t,s)-df(f2s,t)$ "$ 3431 3432 wl2" clearrules subli$ "$ 3433>>$ 3434 3435%------- 3436 3437symbolic procedure sshelp6$ 3438<<wl2"To compute conservation laws the call is:"$ 3439 3440 wl2" ssconl(N,tw,mincw,maxcw,afwlist,abwlist,pdes)$"$ 3441 3442 wl2"where"$ 3443 3444 wl1" N .. the number of superfields theta_i"$ 3445 wl1" tw .. weight(d_t) = 2 times the differential order of the equation"$ 3446 wl1" mincw .. min weight of the conservation law"$ 3447 wl1" maxcw .. max weight of the conservation law"$ 3448 wl1" afwlist .. (algebraic) list of weights of the fermion fields"$ 3449 wl1" f(1), f(2), ... "$ 3450 wl1" abwlist .. (algebraic) list of weights of the boson fields"$ 3451 wl1" b(1), b(2), ..."$ 3452 wl1" pdes .. an algebraic list of the pde(s) for which a conservation"$ 3453 wl1" law is to be found, for example"$ 3454 wl1" {df(f(1),t)=df(f(1),x)*b(1)*p9,"$ 3455 wl1" df(b(1),t)=b(1)**3*p3 + d(1,f(1))**2*p2 + "$ 3456 wl2" df(b(1),x)*b(1)*p9 + df(f(1),x)*f(1)*p4}$"$ 3457 3458 wl1"The number of elements in afwlist and abwlist "$ 3459 wl1"determines the number of fermion and boson fields. All weights are "$ 3460 wl1"based on weight(x)=2."$ 3461 wl1"Important: All fields must be made dependent on x and t before."$ 3462 wl2"Fermion fields are called f(1),f(2),.. and boson fields b(1),.. ."$ 3463 3464 wl1"If a boson weight is non-positive or a fermion weight is negative"$ 3465 wl1"then the global (algebraic mode) variable max_deg must have a positive"$ 3466 wl1"integer value which is the highest degree of such a field or any of"$ 3467 wl2"its derivatives in any ansatz."$ 3468 3469 wl2"-----------------"$ 3470 3471 wl2"Example for determining conservation laws for a given PDE-system:"$ 3472 3473 wl1"ssconl(1,4,10,15,{2},{2},"$ 3474 wl1" {df(f(1),t)=df(f(1),x)*b(1)*p9,"$ 3475 wl1" df(b(1),t)=b(1)**3*p3 + d(1,f(1))**2*p2 + df(b(1),x)*b(1)*p9 "$ 3476 wl1" + df(f(1),x)*f(1)*p4}"$ 3477 wl2" );"$ 3478>>$ 3479 3480%------- 3481 3482symbolic procedure sshelp7$ 3483<<wl1"To compute all possible homogeneities of a list of expressions"$ 3484 wl1"or equations, i.e. to determine for each homogeneity the weight"$ 3485 wl2"of d_t, of fields f(i),b(j) and other constants, the call is:"$ 3486 3487 wl2" FindSSWeights(N,nf,nb,exli,zerowei,verbose)$"$ 3488 3489 wl2"where"$ 3490 3491 wl1" N .. the number of superfields theta_i"$ 3492 wl1" nf .. number of fermion fields f(1) .. f(nf)"$ 3493 wl1" nb .. number of boson fields b(1) .. b(nb)"$ 3494 wl1" exli .. list of equations or expressions"$ 3495 wl1" zerowei .. list of constants or other kernels that"$ 3496 wl1" should have zero weight"$ 3497 wl2" verbose .. =t/nil whether detailed comments shall be made"$ 3498 3499 wl1"Weights are scaled such that weight of d_x is 2, i.e. the weight of"$ 3500 wl2"any D is 1. Input expressions can be in field form or coordinate form."$ 3501 3502 wl1"The program returns a list of homogeneities,"$ 3503 wl1"each homogeneity being a list of"$ 3504 wl1"a list of the weights of f(1),b(2),.. and"$ 3505 wl1"a list of the weights of b(1),b(2),.. and"$ 3506 wl2"a list of the weights of equations/expressions in the input. "$ 3507 3508 3509 wl2"-----------------"$ 3510 3511 wl1"Example for determining all possible weights for a given"$ 3512 wl2"(artificial) system of equations: "$ 3513 3514 wl1"- th(i) denotes superfields theta_i,"$ 3515 wl1"- b's and f's with 3 indices are the coefficients in the"$ 3516 wl1" Taylor expansions"$ 3517 wl1" f(i)=f(i,0,0)+b(i,1,0)*th(1)+b(i,0,1)*th(2)+f(i,1,1)*th(1)*th(2)"$ 3518 wl1"- x,t,th(i) may occur not only as differentiation variables"$ 3519 wl1" but also explicitly polynomially,"$ 3520 wl2"- p3,p5 are supposed to have zero weight."$ 3521 3522 wl1"FindSSWeights(2,2,0,"$ 3523 wl1" {x*df(b(1,1,0),th(1),t)= d(1,d(2,f(2,0,0)))*p5 - "$ 3524 wl1" df(f(1,1,1)*th(1)*th(2),x,2)*p3,"$ 3525 wl1" df(f(2),t)=(df(f(2),x,2)*p3*p5 - df(f(1),x,3)*p3**2)/p5},"$ 3526 wl1" {p3,p5},"$ 3527 wl2" t)$"$ 3528>>$ 3529 3530%------- 3531 3532symbolic procedure sshelp8$ 3533<<wl1"To compute a linearization of a system of evolution equations,"$ 3534 wl2"the call is:"$ 3535 3536 wl2" linearize(pdes,nf,nb,tpar,spar)$"$ 3537 3538 wl2"where"$ 3539 3540 wl1" pdes .. list of equations with a t-derivative on left hand side"$ 3541 wl1" nf .. number of fermion fields f(1) .. f(nf)"$ 3542 wl1" nb .. number of boson fields b(1) .. b(nb)"$ 3543 wl1" tpar .. t/nil, whether time variable t changes parity"$ 3544 wl2" spar .. t/nil, whether symmetry variable s changes parity"$ 3545 3546 wl1" linearize returns {list of relations like df(f(3),s) = f(6),"$ 3547 wl1" list of linearized equations}"$ 3548 wl2"Newly introduced fields depend on x,t,s."$ 3549 3550 wl1"A linearization is associated with a symmetry D_s. The linearized"$ 3551 wl1"equation depends on whether s is parity changing or not and furthermore"$ 3552 wl1"if s is parity changing then the sign of the lhs of the equation (or,"$ 3553 wl2"equally of the rhs) depends on whether t is parity changing or not."$ 3554 3555 wl2"-----------------"$ 3556 3557 wl2"An (artificial) example for a linearization:"$ 3558 3559 wl1"linearize({df(f(1),t)=df(f(2),x)*b(1)**3*p1 "$ 3560 wl1" + d(1,f(2))**5*df(f(1),x,2)*p2,"$ 3561 wl1" df(f(2),t)=2*d(1,df(b(1),x))*df(f(2),x,2)*df(f(1),x)*p3 "$ 3562 wl1" + df(f(1),x,3)*p3**2,"$ 3563 wl2" d(1,d(2,b(1)))=b(1)*d(1,d(2,f(1)))*f(1) },2,1,nil,nil)$"$ 3564>>$ 3565 3566%------- 3567 3568symbolic procedure sshelp9$ 3569<<wl1"To make an ansatz for a polynom in a number of fermionic and/or"$ 3570 wl1"bosonic variables, and their x-derivatives and D(i,..) derivatives"$ 3571 wl2"The call is:"$ 3572 3573 wl2" GenSSPoly(N,wgtlist,cname,mode)$"$ 3574 3575 wl2"where"$ 3576 3577 wl1" N .. the number of superfields theta_i"$ 3578 wl1" wgtlist .. an algebraic mode list of at least one alg. mode list"$ 3579 wl1" {afwlist,abwlist,wgt} where"$ 3580 wl1" afwlist .. (algebraic mode) list of weights of the"$ 3581 wl1" fermion fields f(1), f(2), ..."$ 3582 wl1" abwlist .. (algebraic mode) list of weights of the"$ 3583 wl1" boson fields b(1), b(2), ..."$ 3584 wl1" wgt .. the weight of each term in the generated"$ 3585 wl1" polynomial according to afwlist and abwlist"$ 3586 wl1" The number of elements in afwlist, abwlist determines"$ 3587 wl1" the number of fermion and boson fields."$ 3588 wl1" If the generated polynomial should have more than one"$ 3589 wl1" homogeneity symmetry then each extra one is specified"$ 3590 wl1" by another element {afwlist,abwlist,wgt} in wgtlist."$ 3591 wl1" cname .. name of the coefficients which gets added an index"$ 3592 wl1" mode .. list of flags: "$ 3593 wl1" lin: the generated polynomial has to be linear homogeneous"$ 3594 wl1" in all fields with (index)>(maxindex/2), i.e."$ 3595 wl1" if lin then the number of fermions and boson fields"$ 3596 wl1" must both be even"$ 3597 wl1" fonly: only fermionic terms are generated"$ 3598 wl1" bonly: only bosonic terms are generated"$ 3599 wl1" {forbid,d(1,f(1)),..}: a list of fermionic and bosonic"$ 3600 wl1" variables and their derivatives that may not occur"$ 3601% wl1" {non0coeff,p7,p18,..}: an exclusive list of coefficients"$ 3602% wl2" that may occur"$ 3603 3604 wl2"The result is a list"$ 3605 3606 wl2"{fermonic_polynomial, bosonic_polynomial, {coefficients}}."$ 3607 3608 wl2"Comments:"$ 3609 3610 wl1"- The weight of d/dx is 2, the weight of d(i,..) is 1, no other"$ 3611 wl2" derivatives occur."$ 3612 3613 wl1"- If in the first homogeneity, i.e. in the first element of wgtlist"$ 3614 wl1" a boson weight is non-positive or a fermion weight is negative"$ 3615 wl1" then the global (algebraic mode) variable max_deg must have a"$ 3616 wl1" positive integer value which is the highest degree of such a"$ 3617 wl1" field in any ansatz. Therefore if there is a set of strictly"$ 3618 wl1" positive homogeneity weights, then they should come first, needing"$ 3619 wl2" no extra restriction through max_deg which might restrict generality."$ 3620>>$ 3621 3622%------- 3623 3624symbolic procedure sshelp10$ 3625<<wl1"The N symmetry generators theta_1,..,theta_N are called th(1),..,th(N)"$ 3626 wl1"in SSTools. As they are of odd parity only 2^N many products of their"$ 3627 wl1"powers exist. The name convention used in SSTools for coefficients"$ 3628 wl1"in the Taylor expansion wrt. theta_i becomes apparent in the following"$ 3629 wl2"two expansions (here for N=2):"$ 3630 3631 wl1"f(i)=f(i,0,0)+b(i,1,0)*th(1)+b(i,0,1)*th(2)+f(i,1,1)*th(1)*th(2)"$ 3632 wl2"b(i)=b(i,0,0)+f(i,1,0)*th(1)+f(i,0,1)*th(2)+b(i,1,1)*th(1)*th(2)."$ 3633 3634 wl2"The theta expansion of derivatives D(i,P) is"$ 3635 3636 wl2"D(i,P) = df(P,th(i))+th(i)*df(P,x)} ."$ 3637 3638 wl1"These expansions can be used to convert a polynomial expression in"$ 3639 wl1"terms of f(i), b(j), D(k,..) (what we call `field form') into what"$ 3640 wl2"we call `coordinate form'. The call is:"$ 3641 3642 wl2" tocoo(N,nf,nb,ex)$"$ 3643 3644 wl2"where"$ 3645 3646 wl1" N .. the number of superfields th(1) .. th(N)"$ 3647 wl1" nf .. number of fermion fields f(1) .. f(nf)"$ 3648 wl1" nb .. number of boson fields b(1) .. b(nb)"$ 3649 wl2" ex .. the expression to be converted"$ 3650 3651 wl2"-----------------"$ 3652 3653 wl2"Example: For N=2 and 2 fermion fields the call is"$ 3654 3655 wl1" hh := tocoo(2,2,0,"$ 3656 wl1" f(1)*d(1,f(2))+d(2,d(1,DF(f(2),x)))+"$ 3657 wl2" df(f(2),x,2)*p3*p5 - df(f(1),x,3)*p3**2);"$ 3658>>$ 3659 3660%------- 3661 3662symbolic procedure sshelp11$ 3663<<wl1"The conversion back from coordinate form to field form is done"$ 3664 wl1"by the function tofield. For this function to be applicable "$ 3665 wl1"the expression in coordinate form has to be homogeneous with "$ 3666 wl1"suitable weights and in the current form only involve x- and "$ 3667 wl2"th(i)-derivatives. The call is:"$ 3668 3669 wl2" tofield(N,nf,nb,ex,zerowei)$"$ 3670 3671 wl2"where"$ 3672 3673 wl1" N .. the number of superfields theta_i"$ 3674 wl1" nf .. number of fermion fields f(1) .. f(nf)"$ 3675 wl1" nb .. number of boson fields b(1) .. b(nb)"$ 3676 wl1" ex .. the expression to be converted"$ 3677 wl1" zerowei .. list of constants or other kernels that"$ 3678 wl2" should have zero weight"$ 3679 3680 wl2"-----------------"$ 3681 3682 wl1"Examples:"$ 3683 wl2"The following call inverts the above computation:"$ 3684 3685 wl2" tofield(2,2,0,hh,{p3,p5});"$ 3686 3687 wl1"The next call shows what happens when the input does"$ 3688 wl2"not have a homogeneity:"$ 3689 3690 wl2" tofield(1,1,0,f(1,0)+df(f(1,0),x)+f(1,0)*th(1)*df(f(1,0),x),{});"$ 3691 3692 wl1"In the next call the input has a homogeneity but the resulting"$ 3693 wl1"homogeneity weight for f(1) is negative which does not prevent"$ 3694 wl1"the computation but makes it more expensive but that does not"$ 3695 wl2"matter in the following short example:"$ 3696 3697 wl1" tofield(1,2,0,f(2,0)*b(1,1) + f(2,0)*th(1)*df(f(1,0),x) + "$ 3698 wl2" f(1,0) + th(1)*b(2,1)*b(1,1) + th(1)*b(1,1),{}); "$ 3699 3700 wl1"In the next call the input does have a homogeneity, the resulting"$ 3701 wl2"weight for f(1) is positive but still a field form does not exist:"$ 3702 3703 wl2" tofield(1,1,0,df(f(1,0),x)+f(1,0)*th(1)*df(f(1,0),x),{});"$ 3704>>$ 3705 3706%------- 3707 3708symbolic operator sshelp$ 3709symbolic procedure sshelp$ 3710begin scalar ps,s$ 3711 lines_written:=0$ 3712 rds nil; wrs nil$ % Switch I/O to terminal 3713 ps:=promptstring!*$ promptstring!*:=redfront_color ""$ 3714 repeat << 3715 write"To read about the following topics input the corresponding number:"$ 3716 terpri()$ terpri()$ 3717 write" Purpose (1)"$ terpri()$ 3718 write" Notation (2)"$ terpri()$ 3719 write" Initializations (3)"$ terpri()$ 3720 write" Coefficients (4)"$ terpri()$ 3721 write" Symmetries (5)"$ terpri()$ 3722 write" Conservation laws (6)"$ terpri()$ 3723 write" Homogeneities (7)"$ terpri()$ 3724 write" Linearizations (8)"$ terpri()$ 3725 write" Generating polynomials (9)"$ terpri()$ 3726 write" Taylor expansions (10)"$ terpri()$ 3727 write" Inverting Taylor expansions (11)"$ terpri()$ 3728 write" Exit help (12)"$ terpri()$ 3729 terpri()$ 3730 write"Choice: "$ 3731 s:=read()$ 3732 if ifl!* then rds cadr ifl!*$ % Resets I/O streams 3733 if ofl!* then wrs cdr ofl!*$ 3734 wl2"=========================================================="$ 3735 if s= 1 then sshelp1() else 3736 if s= 2 then sshelp2() else 3737 if s= 3 then sshelp3() else 3738 if s= 4 then sshelp4() else 3739 if s= 5 then sshelp5() else 3740 if s= 6 then sshelp6() else 3741 if s= 7 then sshelp7() else 3742 if s= 8 then sshelp8() else 3743 if s= 9 then sshelp9() else 3744 if s=10 then sshelp10() else 3745 if s=11 then sshelp11() else 3746 if s neq 12 then write"Wrong input, try again."$ 3747 >> until s=12$ 3748 promptstring!*:=ps 3749end$ 3750 3751%------- 3752 3753lisp <<write"For help type: sshelp()$ "$terpri()>>$ 3754 3755endmodule; 3756 3757end$ 3758 3759 3760