1%****************************************************************************** 2% * 3% The program LIEPDE for computing point-, contact- and higher * 4% order symmetries of individual ODEs/PDEs or systems of ODEs/PDEs * 5% * 6% Author: Thomas Wolf * 7% Date: 20. July 1996, 26. Nov 2006 * 8% * 9% For details of how to use LIEPDE see the file LIEPDE.TXT or the * 10% header of the procedure LIEPDE below. * 11% * 12% BSDlicense: * 13% * 14% Redistribution and use in source and binary forms, with or without * 15% modification, are permitted provided that the following conditions are met: * 16% * 17% * Redistributions of source code must retain the relevant copyright * 18% notice, this list of conditions and the following disclaimer. * 19% * Redistributions in binary form must reproduce the above copyright * 20% notice, this list of conditions and the following disclaimer in the * 21% documentation and/or other materials provided with the distribution. * 22% * 23% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * 24% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * 25% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * 26% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR CONTRIBUTORS BE * 27% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR * 28% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF * 29% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS * 30% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN * 31% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * 32% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * 33% POSSIBILITY OF SUCH DAMAGE. * 34%****************************************************************************** 35 36create!-package('(liepde), nil); 37 38symbolic fluid '(print_ logoprint_ adjust_fnc proc_list_ 39 prelim_ individual_ prolong_order !*batch_mode 40 collect_sol flin_ lin_problem done_trafo 41 etamn_al max_gc_fac batch_mode_sub 42 inverse_trafo_list_incomplete)$ 43 44lisp << !*batch_mode:=t$ prelim_:=nil$ 45 individual_:=nil$ prolong_order:=0;etamn_al:=nil >>$ 46 47%---------------------------- 48 49algebraic procedure equ_to_expr(a)$ 50% converts an equation into an expression 51begin scalar lde; 52 return 53 if a=nil then a else 54 <<lisp(lde:=reval algebraic a); 55 if lisp(atom lde) then a else num 56 if lisp(car lde = 'equal) then lhs a - rhs a 57 else a 58 >> 59end$ % of equ_to_expr 60 61 62%******************************************************************** 63module pdesymmetry$ 64%******************************************************************** 65% Routines for finding Symmetries of single or systems of ODEs/PDEs 66% Author: Thomas Wolf 67% July 1996 68 69%--------------------- 70 71% Bei jedem totdf2-Aufruf pruefen, ob evtl. kuerzere dylist reicht 72% evtl die combidiff-Kette und combi nicht mit in dylist, sond. erst in 73% prolong jedesmal frisch generieren. 74 75%--------------------- 76 77symbolic procedure comparedif3(du1,u2,u2list)$ 78% checks whether u2 differentiated wrt. u2list 79% is a derivative of du1 80begin 81 scalar u1l; 82 u1l:=combidif(du1)$ % u1l=(u1, 1, 1, ..) 83 if car u1l neq u2 then return nil else 84 return comparedif1(cdr u1l, u2list) 85end$ % of comparedif3 86 87%--------------------- 88 89%symbolic procedure orderofderiv(du)$ 90%if atom du then (length combidif(du))-1 91% else nil$ 92 93%--------------------- 94 95symbolic procedure mergedepli(li1,li2)$ 96% li1,li2 are lists of lists 97% mergedepli merges the sublists to make one list of lists 98begin scalar newdep; 99 while li1 and li2 do << 100 newdep:=union(car li1, car li2) . newdep; 101 li1:=cdr li1; li2:=cdr li2 102 >>; 103 return if li1 then nconc(reversip newdep,li1) else 104 if li2 then nconc(reversip newdep,li2) else reversip newdep 105end$ 106 107%--------------------- 108 109symbolic procedure adddepli(ex,revdylist)$ 110begin scalar a,b,c,d; 111 for each a in revdylist do << 112 c:=nil; 113 for each b in a do 114 if not my_freeof(ex,b) then c:=b . c; 115 if c or d then d:=c . d; 116 >>; 117 return list(ex,d) 118end$ 119 120%--------------------- 121 122symbolic procedure add_xi_eta_depli(xilist,etalist,revdylist)$ 123begin 124 scalar e1,g,h$ 125 for e1:=1:(length xilist) do << 126 g:=nth(xilist,e1); 127 h:=pnth(g,4); 128 rplaca(h,cadr adddepli(car g,revdylist)) 129 >>; 130 for e1:=1:(length etalist) do << 131 g:=nth(etalist,e1); 132 h:=pnth(g,3); 133 rplaca(h,cadr adddepli(car g,revdylist)) 134 >> 135end$ 136 137%--------------------- 138 139symbolic procedure subtest(uik,sb,xlist,ordok,subordinc)$ 140% uik is a derivative (jet variable) 141% sb is a substitution list without prolongations 142begin 143 scalar el5,el6,el7,el8,el9,el10,sbc$ 144 el5:=combidif(uik); 145 el6:=car el5; el5:=cdr el5; % el6..function name, el5..var.list 146 el7:=nil; el8:=100; el9:=nil; 147 sbc:=sb; 148 while sbc and 149 ((caaar sbc neq el6) or 150 (0 neq <<el7:=comparedif1(cdaar sbc,el5); 151 if el7 and (not zerop el7) and 152 (length(el7)<el8) then 153 <<el8 :=length el7; 154 el9 :=el7; 155 el10:=car sbc 156 >> else el7 157 >>) 158 ) do sbc:=cdr sbc; 159 return 160 if sbc then ((simp!* cadar sbc) . caddar sbc) % simple substitution 161 else 162 if el9 then << % uik is total deriv of car el10 wrt el9 163 uik:=(simp!* cadr el10) . caddr el10; 164 % car uik becomes the expr. to replace the former uik 165 while el9 do << 166 uik:=totdf3(car uik, cdr uik, nth(xlist, car el9), 167 car el9, sb, xlist, ordok, subordinc); 168 el9:=cdr el9 169 >>; 170 uik 171 >> else % no substitution 172 nil 173end$ 174 175%--------------------- 176 177symbolic procedure totdf3(s,depli,x,n,sb,xlist,ordok,subordinc)$ 178% s is the expression to be differentiated wrt. x which is the nth 179% variable in xlist 180% depli is a list of lists, each being the list of jet-variables 181% of order 0,1,2,..., such that s=s(xlist,depli), but 182% as little as possible jet-variables in depli 183% xlist, depli are lisp lists, i.e. without 'list 184% - totdf3 calculates total derivative of s(xlist,depli) w.r.t. x which 185% is the n'th variable, it returns (df(s,x), newdepli) 186% - totdf3 drops jet-variables on which s does not depend 187% - totdf3 automatically does substitutions using the list sb which 188% is updated if prolongations of substitutions are calculated, 189% i.e. sb is destructively changed! 190% - structure of sb: lisp list of lisp lists: 191% ((to_be_replaced_jet_var_name, to_be_replaced_jet_var_deriv_1,..), 192% subst_expr_in_jet_space_coord, list_of_jet_vars_in_subst_expr) 193% - subordinc is a number by how much the order may increase due to 194% substitutions sb. 195% - ordok is the lowest order which must be accurate. If ordok>0 and 196% s is of lower order than ordok then from depli only derivatives 197% of order ordok-1-subordinc to ordok-1 are used. 198begin 199 scalar tdf,el1,el2,el3,el4,el5,newdepli, 200 newdy,dy,ddy$ 201 newdepli:=nil; % the new dependence list 202 newdy:=nil; % the new dep.list due to chain rule 203 ddy:=nil; % ddy .. derivatives of jet-variables 204 % resulting from diff. of lower order 205 %--- Should only terms in the result be acurate that include 206 %--- derivatives of order>=ordok? 207 if ordok>0 then << 208 tdf:=simp!* 0; 209 depli:=copy depli; 210 el2:=length depli; 211 if el2<(ordok-subordinc) then depli:=nil 212 else 213 for el1:=1:(ordok-1-subordinc) do << 214 dy:=pnth(depli,el1); 215 rplaca(dy,nil); 216 >> 217 >> else tdf:=diffsq(s,x); 218 %--- The differentiations wrt. u-derivatives 219 for each el1 in depli do % for each order do 220 <<dy:=union(ddy,el1); ddy:=nil;% dy .. occuring jet-var. of this order 221 while el1 do 222 <<el2:=car el1; el1:=cdr el1;% el2 is one jet-variable of this order 223 el3:=diffsq(s,el2); 224 if null car el3 then % It is tempting to do dy:=delete(el2,dy) but 225 % that can lead to errors because, dy is a list of the derivatives that 226 % appear in the total derivative of s, not necessarily in s itself. That 227 % means, if s depends on a 2nd and 4th order derivative but not on a 3rd 228 % order derivative, then the third order derivative may just have been 229 % added through dy:=union(ddy,el1); so the fact that el3 is nil, i.e. s 230 % does not depend on the 3rd order derivative should not result in 231 % deleting the 3rd order derivative from dy. 232 else << 233 el4:=dif(el2,n); % el4=df(el2,x) 234 %----- Is el4 to be substituted by sb? 235 if el5:=subtest(el4,sb,xlist,ordok,subordinc) then << 236 el4:=car el5; 237 newdepli:=mergedepli(newdepli,cdr el5) 238 >> else <<ddy:=el4 . ddy;el4:=simp!* el4>>; 239 tdf:=addsq(tdf, multsq(el4, el3)) 240 >> 241 >>; 242 newdy:=dy . newdy 243 >>; 244 if ddy then newdy:=ddy . newdy; 245 246 newdepli:=mergedepli(reversip newdy,newdepli); 247 return (tdf . newdepli) 248end$ % of totdf3 249 250%--------------------- 251 252symbolic procedure joinsublists(a)$ 253% It is assumed, a is either nil or a list of lists or nils which 254% have to be joined 255if null a then nil 256 else append(car a,joinsublists(cdr a))$ 257 258%--------------------- 259 260algebraic procedure transeq(eqn,xlist,ylist,sb)$ 261<<for each el1 in sb do eqn:=sub(el1,eqn); 262 for each el1 in ylist do 263 for each el2 in xlist do nodepend el1,el2; 264 eqn>>$ 265 266%--------------------- 267 268symbolic operator drop$ 269symbolic procedure drop(a,vl)$ 270% liefert summe aller terme aus a, die von elementen von vl abhaengen 271begin scalar b$ 272 if not((pairp a) and (car a='plus)) then b:=a 273 else 274 <<vl:=cdr vl; % because vl is an algebraic list 275 for each c in cdr a do 276 if not freeoflist(c,vl) then b:=cons(c,b)$ 277 if b then b:=cons('plus,reverse b)>>$ 278 return b$ 279end$ 280 281%--------------------- 282 283symbolic procedure etamn(u,indxlist,xilist,etalist, 284 ordok,truesub,subordinc,xlist)$ 285% determines etamn recursively 286% At the end, ulist= list of df(u,i,cdr indxlist) for all i 287begin 288 scalar etam,x,h1,h2,h3,h4,h5,ulist,el,r,cplist,depli; 289 290 % Has it already been computed? 291%if nil then 292 if ordok=0 then << 293 h5:=u$ h2:=indxlist$ 294 while h2 do <<h5:=mkid(h5,car h2); h2:=cdr h2>>; 295 h3:=assoc(h5,etamn_al) 296 >>$ 297 if h3 then return cdr h3$ 298 299 if (null indxlist) or ((length indxlist)=1) then 300 <<cplist:=etalist; 301 while u neq cadar cplist do cplist:=cdr cplist; 302 etam:=(caar cplist . caddar cplist) . nil; 303 >> else 304 etam:=etamn(u,cdr indxlist,xilist,etalist,ordok,truesub, 305 subordinc,xlist)$ 306 307 return 308 309 if null indxlist then etam 310 else << 311 ulist:=nil; 312 x:=cdr nth(xilist,car indxlist); % e.g. x := (v3,3,dylist) 313 r:=if null car caar etam then <<depli:=nil;caar etam>> 314 else << 315 h2:=totdf3(caar etam,cdar etam,car x,cadr x,truesub,xlist, 316 ordok,subordinc)$ 317 depli:=cdr h2; 318 car h2 319 >>; 320 etam:=cdr etam; % = reverse ulist 321 cplist:=xilist; 322 323 h3:=nil; 324 while cplist do 325 <<el:=car cplist; % e.g. el=(xi_bz1 bz1 2 ((u))) 326 cplist:=cdr cplist; 327 if (length indxlist)=1 then h1:=dif(u,caddr el) 328 else << 329 h1:=dif(car etam,cadr indxlist); % e.g. h1:=u!`i!`n 330 etam:=cdr etam; 331 >>; 332 333 ulist:=h1 . ulist; 334 if not sqzerop car el then << 335 336 %--- substitution of h1? 337 if h4:=subtest(h1,truesub,xlist,ordok,subordinc) 338 then <<h1:=car h4;depli:=mergedepli(depli,cdr h4)>> 339 else h1:=simp!* h1; 340 r:=subtrsq(r,multsq(h1,<<h2:=totdf3(car el,cadddr el,car x, 341 cadr x,truesub,xlist,0,0)$ 342 if car car h2 then % ie. car h2 neq 0 343 <<if h4 then depli:=mergedepli(depli,cdr h4) 344 else h3:=prepsq h1 . h3; 345 depli:=mergedepli(depli,cdr h2); 346 >>$ 347 car h2 348 >> 349 ) 350 ); 351 >> 352 >>; 353 if h3 then << 354 h3:=list h3; 355 for h2:=1:(length indxlist) do h3:=nil . h3; 356 depli:=mergedepli(depli,h3); 357 >>; 358 % (if not full then drop(r,'list . car revdylist) else r) . 359 % (reverse ulist) 360 h1:=(r . depli) . (reverse ulist)$ 361 if ordok=0 then etamn_al:=cons((h5 . h1),etamn_al); 362 h1 363 >> 364end$ % of etamn 365 366%--------------------- 367 368symbolic procedure callcrack(!*time,cpu,gc,lietrace_,symcon, 369 flist,vl,xilist,etalist,inequ,last_call)$ 370begin 371 scalar g,h,oldbatch_mode,print_old; 372 if !*time then <<terpri()$ 373 write "time to formulate conditions: ", time() - cpu, 374 " ms GC time : ", gctime() - gc," ms"$ 375 >>; 376 if lietrace_ then algebraic << 377 write"Symmetry conditions before CRACK: "; 378 write lisp ('list . symcon); 379 >>; 380 oldbatch_mode:=!*batch_mode$ 381 !*batch_mode:=batch_mode_sub$ 382 print_old:=print_$ 383 if null batch_mode_sub and null print_ then print_:=8$ 384 385 % lex_df=nil gives shorter computation but lex_df=t gives system 386 % that is easier to integrate, eg ODEs it they are in the diff ideal 387 % so computation starts as lex_df=nil and is changed to lex_df=t at 388 % end if necessary 389 if freeof(proc_list_,'try_other_ordering) then 390 proc_list_:=append(proc_list_,list 'try_other_ordering)$ 391 h:=sq!*crack({'list . symcon,'list . inequ,'list . flist,'list . vl}); 392 393 !*batch_mode:=oldbatch_mode$ 394 print_:=print_old$ 395 396 if last_call then return h; 397 398 if h neq list('list) then 399 <<h:=cadr h; 400 symcon:=cdadr h$ 401 for each g in cdaddr h do << 402 xilist :=for each k in xilist collect cons(subsq(car k,{(cadr g . caddr g)}),cdr k)$ %!! 403 etalist:=for each k in etalist collect cons(subsq(car k,{(cadr g . caddr g)}),cdr k)$ %!! 404 inequ :=subst(caddr g, cadr g, inequ); 405%--> Erkennung von 'e, 'x siehe: 406 407% h:=intern car explode cadr e1; 408%write"h=",h;terpri()$ 409% if (h='x) or (h='X) then 410% xilist :=subst(caddr e1, cadr e1, xilist) else 411% if (h='e) or (h='E) or (h="e") or (h="E") then 412% etalist:=subst(caddr e1, cadr e1, etalist) else 413% rederr("One ansatz does not specify XI_ nor ETA_.") 414 >>; 415 416 if lietrace_ then << 417 write"symcon nachher: ",symcon; 418 write"xilist=", xilist; 419 write"etalist=", etalist; 420 >>; 421 flist:=cdr reval cadddr h; 422 if print_ then 423 <<terpri()$ 424 write"Remaining free functions after the last CRACK-run:"; 425 terpri()$ 426 fctprint flist;terpri()$terpri()>>; 427 >>; 428 429 return list(symcon,xilist,etalist,flist,inequ) 430end$ % of callcrack 431 432%--------------------- 433 434put('liepde,'psopfn,'liepde_psopfn)$ 435symbolic procedure liepde_psopfn(inp)$ 436 437 comment 438 439 problem: {{eq1,eq2, ...}, % equations 440 { y1, y2, ...}, % functions 441 { x1, x2, ...} } % variables 442 443 % Equations `eqi' can be given as single differential 444 % expressions which have to vanish or they can be given 445 % in the form df(..,..[,..]) = .. . 446 447 % If the equations are given as single differential 448 % expressions then the program will try to bring it into 449 % the `solved form' ..=.. automatically. 450 451 % The solved forms (either from input or generated within 452 % LIEPDE) will be used for substitutions, to find 453 % all symmetries satisfied by solutions of the equations. 454 % Sufficient conditions for this procedure to be correct, 455 % i.e. to get *all* symmetries of the specified type on the 456 % solution space are: 457 458 % - There are equally many equations and functions. 459 % - Each function is used once for a substitution and 460 % each equation is used once for a substitution. 461 % - All functions differentiated on the left hand sides 462 % (lhs) depend on all variables. 463 % - In no equation df(..,..[,..]) = .. does the right hand 464 % side contain the derivative on the lhs nor any 465 % derivative of it. 466 % - No equation does contain a lhs or the derivative 467 % of a lhs of another equation. 468 469 % These conditions are checked in LIEPDE and execution 470 % is stoped if they are not satisfied, i.e. if the input 471 % was not correct, or if the program was not able to 472 % write the input expressions properly the solved 473 % form ..=.. . One then should find for each function 474 % one derivative which occurs linearly in one equation. 475 % The chosen derivatives should be as high as possible, 476 % at least there must no derivative of them occur in any 477 % equation. An easy way to get the equations in the 478 % desired form is to use 479 % FIRST SOLVE({eq1,eq2,...},{list of derivatives}) 480 481 % NOTE that to improve efficiency it is advisable not to 482 % express lower order derivatives on the left hand side 483 % through higher order derivatives on the right hand side. 484 % SEE also the implications on completeness for the 485 % determination of generalized symmetries with 486 % characteristic functions of a given order, as described 487 % below and the two examples with the Burgers equation. 488 489 symtype: {"point"} % for point symmetries 490 {"contact"} % for contact symmetries, is only 491 % applicable if only one function, 492 % only one equation of order>1 493 {"general",order} % for generalized symmetries of 494 % order `order' which is an integer>0 495 % NOTE: Characteristic functions of 496 % generalized symmetries (i.e. the 497 % eta_.. if xi_..=0) are equivalent 498 % if they are equal on the solution 499 % manifold. Therefore all dependencies 500 % of characteristic functions on 501 % the substituted derivatives and their 502 % derivatives are dropped. This has the 503 % consequence that if, e.g. for the heat 504 % equation df(u,t)=df(u,x,2), df(u,t) is 505 % substituted by df(u,x,2) then 506 % {"general",2) would not include 507 % characteristic functions depending 508 % on df(u,t,x), or df(u,x,3). THEREFORE: 509 % If you want to find all symmetries up 510 % to a given order then 511 % - either avoid substituting lower 512 % order derivatives by expressions 513 % involving higher derivatives, or, 514 % - go up in the order specified in 515 % symtype. 516 % 517 % Example: 518 % 519 % depend u,t,x 520 % liepde({{df(u,t)=df(u,x,2)+df(u,x)**2}, 521 % {u},{t,x}}, 522 % {"general",3},{}) 523 % 524 % will give 10 symmetries + one infinite 525 % family of symmetries whereas 526 % 527 % liepde({{df(u,x,2)=df(u,t)-df(u,x)**2}, 528 % {u},{t,x}}, 529 % {"general",3},{}) 530 % 531 % will give 28 symmetries + one infinite 532 % family of symmetries. 533 534 {xi!_x1 =..., 535 ... 536 eta!_y3=... } % - An ansatz must specify all xi!_.. for 537 % all indep. variables and all eta!_.. 538 % for all dep. variables in terms of 539 % differential expressions which may 540 % involve unknown functions/constants. 541 % The dependencies of the unknown 542 % functions have to declared earlier 543 % using the DEPEND command. 544 % - If the ansatz should consist of the 545 % characteristic functions then set 546 % all xi!_..=0 and assign the charac- 547 % teristic functions to the eta!_.. . 548 549 flist: {- all parameters and functions in the equations which are to 550 be determined, such that there exist symmetries, 551 - if an ansatz has been made in symtype then flist contains 552 all unknown functions and constants in xi!_.. and eta!_..} 553 554 inequ: {all non-vanishing expressions which represent 555 inequalities for the functions in flist} 556 557 Further comments: 558 559 The syntax of input is the usual REDUCE syntax. For example, the 560 derivative of y3 wrt x1 once and x2 twice would be df(y3,x1,x2,2). 561 --> One exception: If in the equations or in the ansatz the dependence 562 of a free function F on a derivative, like df(y3,x1,x2,2) shall be 563 declared then the declaration has to have the form: 564 DEPEND F, Y3!`1!`2!`2 565 - the ! has to preceede each special character, like `, 566 - `i stands for the derivative with respect to the i'th variable in 567 the list of variables (the third list in the problem above) 568 569 If the flag individual_ is t then conditions are investigated for 570 each equation of a system of equations at first individually before 571 conditions resulting from other equations are added. 572 573 If the flag prelim_ is t then preliminary conditions for equations 574 of higher than 1st order are formulated and investigated before the 575 full condition is formulated and investigated by CRACK. 576 577 If the REDUCE switch TIME is set on with ON TIME then times for the 578 individual steps in the calculation are shown. 579 580 Further switches and parameters which can be changed to affect the 581 output and performance of CRACK in solving the symmetry conditions 582 are listed in the file CRINIT.RED. 583; 584 585begin 586 scalar cpu, gc, lietrace_, oldadj, eqlist, ylist, xlist, pointp, 587 contactp, generalp, ansatzp, symord, e1, e2, ordr, sb, 588 dylist, revdylist, xi, eta, eqordr, eqordrcop, no, eqcopy1, 589 truesub, deplist, xilist, etalist, dycopy, freelist, eqlen, 590 dylen, truesubno, minordr, n1, n2, n3, n4, n5, n, h, jetord, 591 allsub, subdy, lhslist, symcon, subordinc, eqn, depli, vl, occli, 592 revdycopy, subordinclist, xicop, etacop, flcop, etapqlist, 593 etapqcop, etapq, oldbatch_mode, allsym, symcon_s, xilist_s, 594 etalist_s, inequ_s, flist_s, truesub_s, oldcollect_sol, 595 oldprint_, flist_slin, flist_snli, return_list, last_call, 596 paralist, proc_list_bak, max_gc_fac_bak, flistorg,problem, 597 symtype,flist,inequ$ 598 599 % lietrace_:=t; 600 backup_reduce_flags()$ 601 cpu:=time()$ gc:=gctime()$ 602 oldadj:=adjust_fnc; adjust_fnc:=nil; 603 oldcollect_sol:=collect_sol; collect_sol:=t; 604 605 %--------- extracting input data 606 problem:=aeval car inp$ 607 symtype:=reval cadr inp$ 608 flist :=reval caddr inp$ 609 inequ :=reval cadddr inp$ 610 eqlist:= 611 if atom cadr problem then list cadr problem else 612 if car cadr problem='list then cdr cadr problem else 613 list cadr problem$ 614 ylist := reval maklist caddr problem; 615 xlist := reval maklist cadddr problem; 616 617 if inequ then inequ:=cdr inequ; 618 if flist then flist:=cdr flist; 619 620 % Is this a non-linear problem? 621 e1:=flist; 622 while e1 and freeof(eqlist,car e1) do e1:=cdr e1; 623 lin_problem:=if e1 then nil else t$ 624 625 eqlen:=length eqlist; 626 627 % if 1+eqlen neq length(ylist) then rederr( <--- taken out 628 % "Number of equations does not match number of unknown functions."); <--- taken out 629 630 for each e1 in cdr ylist do 631 for each e2 in cdr xlist do 632 if my_freeof(e1,e2) then rederr( 633 "Not all functions do depend on all variables."); 634 635 if atom cadr symtype then % default case 636 if cadr symtype = "point" then <<pointp :=t;symord:=0>> else 637 if cadr symtype = "contact" then <<contactp:=t;symord:=1; 638 639 % if eqlen>1 then rederr( <--- taken out 640 % "Contact symmetries only in case of one equation for one function.") <--- taken out 641 >> else 642 if cadr symtype = "general" then <<generalp:=t;symord:=caddr symtype; 643 if (not fixp symord) or (symord<1) then rederr( 644 "The order of the generalized symmetry must be an integer > 0.") 645 >> else rederr("Inconclusive symmetry type.") 646 else << 647 ansatzp:=t; % an ansatz has been made 648 % if length(ylist)+length(xlist) neq length(symtype)+1 then <--- taken out 649 % rederr("Number of assignments in the ansatz is wrong."); <--- taken out 650 651 % find the order of the ansatz 652 symord:=0; 653 % at first the max order of derivatives of elements of ylist 654 for each e1 in cdr symtype do 655 for each e2 in cdr ylist do 656 <<n:=totdeg(caddr e1,e2); 657 if n>symord then symord:=n>>; 658 % then the max order of jet variables occuring in functions to be computed 659 for each e1 in flist do << 660 e2:=fctargs e1$ 661 for each h in e2 do <<n2:=-1+length combidif h;if n2>symord then symord:=n2>>$ 662 >>$ 663 if symord = 0 then pointp :=t else 664 if (symord = 1) and (length(ylist)=2) then contactp:=t else 665 generalp:=t; 666 sb:=nil; 667 for each e1 in flist do if freeof(eqlist,e1) then sb:=cons({'equal,e1,0},sb)$ 668 sb:=cons('list,sb); 669 h:=nil; 670 for each e1 in cdr symtype do << 671 n1:=algebraic(sub(lisp sb,lisp caddr e1))$ % substitution on the rhs of the sym.ansatz 672 if not zerop n1 and 673 (not pairp n1 or 674 (pairp n1 and ((car n1 neq 'equal) or (not zerop caddr n1)))) then << 675 h:=0; 676 write"Your ansatz for the symmetry needs to be homogeneous, i.e. ", 677 "substituting all unknown functions and constants to be computed ", 678 "(which do not occur in the equation) to zero needs to make the ", 679 "symmetry to zero. In your ansatz this is not ", 680 "the case because the list of substitutions:"$ 681 algebraic(write lisp sb)$ 682 write"leaves this right hand side non-zero:"$ 683 algebraic(write lisp {'equal,cadr e1,n1})$ 684 write"To fix your ansatz you could, for example, simply multiply all ", 685 "non-vanishing parts on all right hand sides in your ansatz with one ", 686 "and the same unknown constant, say cc, and add cc to the list of unknowns ", 687 "to be computed and to the list of non-vanishing expressions." 688 >> 689 >> 690 >>$ 691 if h then return nil$ 692 693 problem:=0; 694 695 %---- Are substitutions already given in the input? 696 eqcopy1:=eqlist; 697 while eqcopy1 and (pairp car eqcopy1) and (caar eqcopy1='equal) and 698 (pairp cadar eqcopy1) and (caadar eqcopy1='df) do 699 eqcopy1:=cdr eqcopy1; 700 if null eqcopy1 then truesub:=eqlist; 701 eqcopy1:=nil; 702 703 %--------- initial printout 704 if print_ and logoprint_ then <<terpri()$ 705 write "-----------------------------------------------", 706 "---------------------------"$ terpri()$terpri()$ 707 write"This is LIEPDE - a program for calculating infinitesimal", 708 " symmetries"; terpri()$ 709 write "of single differential equations or systems of de's"; 710 >>; 711 terpri();terpri(); 712 if length xlist=2 then write"The ODE" 713 else write"The PDE"; 714 if length ylist>2 then write"-system"; 715 write " under investigation is :";terpri(); 716 for each e1 in eqlist do algebraic write lisp e1; 717 terpri()$write "for the function(s) : ";terpri()$ 718 terpri()$fctprint cdr reval ylist; 719 terpri()$terpri(); 720 eqlist:=for each e1 in eqlist collect algebraic equ_to_expr(lisp e1); 721 if eqlen > 1 then eqlist:=desort eqlist; 722 723 if !*time then <<terpri()$ 724 terpri()$terpri()$ 725 write"=============== Initializations" ; 726 >>; 727 728 ordr:=0; 729 for each e1 in eqlist do << 730 h:=0; 731 for each e2 in cdr ylist do 732 <<n:=totdeg(e1,e2); 733 if n>h then h:=n>>; 734 eqordr:=h . eqordr; 735 if h>ordr then ordr:=h 736 >>; 737 eqordr:=reversip eqordr; 738 739 if ordr>symord then jetord:=ordr 740 else jetord:=symord$ 741 sb:=subdif1(xlist,ylist,jetord)$ 742 eqlist:=cons('list,eqlist); 743 if ansatzp then eqlist:=list('list,symtype,eqlist); 744 if truesub then eqlist:=list('list,cons('list,truesub),eqlist); 745 if inequ then eqlist:=list('list,cons('list,inequ),eqlist); 746 747 on evallhseqp; 748 eqlist:=transeq(eqlist,xlist,ylist,sb); %<<<<<<----- slow 749 off evallhseqp; 750 751 if inequ then <<inequ :=cdadr eqlist;eqlist:=caddr eqlist>>; 752 if truesub then <<truesub:=cdadr eqlist;eqlist:=caddr eqlist>>; 753 if ansatzp then <<symtype:=cdadr eqlist;eqlist:=cdaddr eqlist>> 754 else eqlist:=cdr eqlist; 755 756 ylist:=cdr ylist; 757 xlist:=cdr xlist; 758 759 if lietrace_ and ansatzp then write"ansatz=",symtype; 760 761 dylist:=ylist . reverse for each e1 in cdr sb collect 762 for each e2 in cdr e1 collect caddr e2; 763 revdylist:=reverse dylist; % dylist with decreasing order 764 765 vl:=xlist; 766 for each e1 in dylist do vl:=append(e1,vl); 767 vl:='list . vl; 768 769 if not ansatzp then 770 deplist:=for n:=0:symord collect nth(dylist,n+1); 771 % list of variables the xi_, eta_ depend on 772 773 xi :=reval algebraic xi!_; 774 eta:=reval algebraic eta!_; 775 n:=0; 776 xilist :=for each e1 in xlist collect 777 <<n:=n+1; 778 if pointp or ansatzp then << 779 h:=mkid('xi_,e1); %!! 780 if not ansatzp then << 781 nodependlist {h}; 782 depnd(h,xlist . deplist); 783 flist:=h . flist; 784 flin_:=h . flin_; 785 depli:=deplist; 786 >> else depli:=nil 787 >> else <<h:=0;depli:=nil>>; 788 {simp h,e1,n,depli} 789 >>; 790 depli:=if (not ansatzp) and (not generalp) then deplist 791 else nil; 792 n:=0; 793 etalist:=for each e1 in ylist collect 794 <<n:=n+1; 795 h:=mkid('eta_,e1); %!! 796 if not ansatzp then << 797 if not generalp then <<nodependlist {h};depnd(h,xlist . deplist)>>; 798 % the generalp-case is done below when substitutions are known 799 flist:=h . flist; 800 flin_:=h . flin_; 801 >>; 802 {simp h,e1,depli} 803 >>; 804 flistorg:=flist; 805 806 if ansatzp then << 807 for each e1 in symtype do << 808 xilist :=for each k in xilist collect cons(subsq(car k,{(cadr e1 . caddr e1)}),cdr k)$ %!! 809 etalist:=for each k in etalist collect cons(subsq(car k,{(cadr e1 . caddr e1)}),cdr k)$ %!! 810 %--> Erkennung von 'e, 'x siehe: 811 % h:=intern car explode cadr e1; 812 %write"h=",h;terpri()$ 813 % if (h='x) or (h='X) then 814 % xilist :=subst(caddr e1, cadr e1, xilist) else 815 % if (h='e) or (h='E) or (h="e") or (h="E") then 816 % etalist:=subst(caddr e1, cadr e1, etalist) else 817 % rederr("One ansatz does not specify XI_ nor ETA_.") 818 >>; 819 add_xi_eta_depli(xilist,etalist,revdylist)$ 820 >>; 821 822 if lietrace_ then write"xilist=",xilist," etalist=",etalist; 823 %---- Determining a substitution list for highest derivatives 824 %---- from eqlist. Substitutions may not be optimal if starting 825 %---- system is not in standard form 826 827 comment: Counting in how many equations each highest 828 derivative occurs. Those which do not occur allow Stephani-Trick, 829 those which do occur once and there linearly are substituted by that 830 equation. 831 832 Because one derivative shall be assigned it must be one of 833 the highest derivatives from each equation used, or one such that 834 no other derivative in the equation is a derivative of it. 835 836 Each equation must be used only once. 837 838 Each derivative must be substituted by only one equation. 839 840 At first determining the number of occurences of each highest 841 derivative. 842 843 The possiblity of substitutions is checked in each total derivative. 844 845 $ 846 847 if truesub then << %--- determination of freelist %and occurlist 848 dycopy:=car revdylist; %--- the highest derivatives 849 while dycopy do 850 <<e1:=car dycopy; dycopy:=cdr dycopy; 851 eqcopy1:=eqlist; 852 while eqcopy1 and my_freeof(car eqcopy1,e1) do 853 eqcopy1:=cdr eqcopy1; 854 855 if null eqcopy1 then freelist :=e1 . freelist 856 %else occurlist:=e1 . occurlist; 857 >> 858 >> else << 859 860 no:=0; % counter of the following repeat-loop 861 % freelist (and occurlist) are determined 862 % only in the first run 863 eqordrcop:=copy eqordr; 864 % for bookkeeping which equation have been used 865 866 repeat << 867 no:=no+1; %--- incrementing the loop counter 868 869 %--- truesubno is the number of substitutions so far found. 870 %--- It is necessary at the end to check whether new substitutions 871 %--- have been found. 872 if null truesub then truesubno:=0 873 else truesubno:=length truesub; 874 %--- substitutions of equations of minimal order are searched first 875 minordr:=1000; %--- minimal order of the so far unused equations 876 for each e1 in eqordrcop do 877 if (e1 neq 0) and (e1<minordr) then minordr:=e1; 878 879 dycopy:=copy nth(dylist,minordr+1); %-- all deriv. of order minordr 880 dylen:=length dycopy; 881 882 allsub:=nil; 883 884 for n1:=1:dylen do %--- checking all deriv. of order minordr 885 <<e1:=nth(dycopy,n1); %--- e1 is the current candidate 886 %--- here test, whether e1 is not a derivative of a lhs of one 887 %--- of the substitutions car e2 found so far 888 h:=combidif(e1); n:=car h; h:=cdr h; 889 e2:=truesub; 890 while e2 and (null comparedif3(cadar e2,n,h)) do e2:=cdr e2; 891 if null e2 then << 892 893 n2:=0; %-- number of equations in which the derivative e1 occurs 894 subdy:=nil; 895 for n3:=1:eqlen do 896 if not my_freeof(nth(eqlist,n3),e1) then 897 % here should also be tested whether derivatives of e1 occur 898 % and not just my_freeof 899 %--> 900 <<n2:=n2+1; 901 if nth(eqordrcop,n3)=minordr then 902 %--- equation is not used yet and of the right order 903 <<e2:=cdr algebraic coeff(lisp nth(eqlist,n3),lisp e1); 904 if hipow!*=1 then 905 subdy:=list(n1,n3,list('equal,e1,list('minus, 906 list('quotient, car e2, cadr e2)))) . subdy 907 >> 908 >>; 909 if n2=0 then if no=1 then freelist:=e1 . freelist else 910 else 911 <<%if no=1 then occurlist:=e1 . occurlist; 912 if subdy then if n2=1 then 913 <<h:=car subdy; 914 truesub:=(caddr h) . truesub; 915 n:=pnth(dycopy ,car h);rplaca(n,0); 916 n:=pnth(eqordrcop,cadr h);rplaca(n,0); 917 >> else 918 allsub:=nconc(allsub,subdy); 919 >> 920 >> 921 >>; 922 923 %---- Taking the remaining known substitutions of highest deriv. 924 h:=subdy:=0; 925 for each h in allsub do 926 if (nth(dycopy , car h) neq 0) and 927 (nth(eqordrcop,cadr h) neq 0) then 928 <<truesub:=(caddr h) . truesub; 929 n:=pnth(dycopy ,car h);rplaca(n,0); 930 n:=pnth(eqordrcop,cadr h);rplaca(n,0); 931 >>; 932 >> until (truesub and (length(truesub)=eqlen)) % complete 933 or (truesubno=length(truesub))$ % no progress 934 allsub:=eqordrcop:=dycopy:=nil; 935 936 if (null truesub) or 937 (eqlen neq length(truesub)) then rederr( 938 "Unable to find all substitutions. Input equations as df(..,..)=..!"); 939 >>; 940 941 lhslist:=for each e1 in truesub collect cadr e1; 942 943 %------ check that functions to be computed do not depend on 944 % lhs of the system 945 chkflist(cons('list,flist),cons('list,lhslist))$ 946 947 %-- Bringing truesub into a specific form: lisp list of lisp lists: 948 % ((to_be_replaced_jet_var_name, to_be_replaced_jet_var_deriv_1,..), 949 % subst_expr_in_jet_space_coord, list_of_jet_vars_in_subst_expr) 950 truesub:=for each e1 in truesub collect 951 cons(combidif cadr e1, adddepli(caddr e1,revdylist))$ 952 953 %--- Checking that no rhs of a substitution contains any lhs or 954 %--- derivative of a lhs 955 h:=t; %--- h=nil if lhs's are derivatives of each other 956 no:=t; %--- no=nil if one lhs can be substituted in a rhs 957 for each e1 in truesub do 958 if h and no then << 959 n1:=caar e1; n2:=cdar e1; dylen:=length n2; 960 for each e2 in truesub do << 961 %--- comparison of two lhs's 962 if not(e1 eq e2) and (n1=caar e2) and 963 comparedif1(n2,cdar e2) then h:=nil; %--- truesub is not ok 964 %--- can the lhs of e1 be substituted on the rhs? 965 dycopy:=caddr e2; 966 for n:=1:dylen do if dycopy then dycopy:=cdr dycopy; 967 for each e3 in dycopy do 968 for each e4 in e3 do 969 if comparedif2(n1,n2,e4) then no:=nil; 970 >> 971 >>; 972 if null h then rederr( 973 "One substitution can be made in the lhs of another substitution!"); 974 if null no then rederr( 975 "One substitution can be made in the rhs of another substitution!"); 976 977 %???????????????????????????????????????????? 978 % %--- Checking that a derivative of each dependent variable is 979 % %--- substituted once. This is a sufficient condition for having 980 % %--- a de-system that is a differential Groebner basis 981 % h:=nil; 982 % for each e1 in lhslist do h:=adjoin(car combidif e1,h); 983 % if length(h) neq length(lhslist) then rederr( 984 % "For at least one function there is more that one substituion!")$ 985 986 %--- Determine of how much the order may increase by a substitution 987 subordinc:=0; 988 subordinclist:=for each h in truesub collect << 989 n:=(length caddr h) - (length car h); 990 if n>subordinc then subordinc:=n; 991 n 992 >>; 993 if lietrace_ then <<terpri()$write"truesub=",truesub; 994 terpri()$write"freelist=",freelist; 995 %terpri()$write"occurlist=",occurlist 996 >>; 997 %--- To avoid non-uniqueness in the case of the investigation of 998 %--- generalized symmetries let the characteristics eta_.. (xi_..=0) 999 %--- not depend on substituted derivatives 1000 if generalp and (null ansatzp) then << 1001 deplist:=ylist . 1002 for each dycopy in cdr deplist collect << 1003 for each h in lhslist do 1004 %---- delete h and derivatives of h 1005 dycopy:=listdifdif1(h,dycopy); 1006 dycopy 1007 >>; 1008 for e1:=1:(length etalist) do << 1009 h:=nth(etalist,e1); 1010 e2:=prepsq car h$ 1011 nodependlist {e2}; 1012 depnd(e2,xlist . deplist); 1013 h:=pnth(h,3); 1014 rplaca(h,deplist) 1015 >> 1016 >>; 1017 % reduced set of solution techniques for preliminary conditions 1018 proc_list_:=delete('multintfac,proc_list_)$ 1019 if !*time then <<terpri()$ 1020 write "time for initializations: ", time() - cpu, 1021 " ms GC time : ", gctime() - gc," ms"$ 1022 cpu:=time()$ gc:=gctime()$ 1023 >>; 1024 %------ Determining first short determining equations and solving them 1025 if prelim_ or individual_ then << 1026 proc_list_bak:=proc_list_$ 1027 proc_list_:='(to_do 1028 separation 1029 subst_level_0 1030 subst_level_03 1031 quick_integration 1032 subst_level_33 1033 gen_separation 1034 % full_integration 1035 )$ 1036 max_gc_fac_bak:=max_gc_fac$ 1037 max_gc_fac:=0; 1038 inverse_trafo_list_incomplete:=nil 1039 >>$ 1040 symcon:=nil; 1041 n1:=0; 1042 if prelim_ and lin_problem then 1043 for each eqn in eqlist do 1044 <<n1:=n1+1; 1045 if !*time then <<terpri()$ 1046 terpri()$terpri()$ 1047 write"=============== Preconditions for the ",n1,". equation" ; 1048 >>; 1049 revdycopy:=revdylist; 1050 for e1:=(nth(eqordr,n1) + 1):ordr do revdycopy:=cdr revdycopy; 1051 n2:=cadr adddepli(eqn,revdycopy); % jet-variables in eqn 1052 vl:=n2; 1053 occli:=lastcar n2; 1054 freelist:=setdiff(car revdycopy,occli); 1055 if pointp and (subordinc=0) then 1056 eqn:=drop(eqn,occli) % dropp all terms without a highest deriv. 1057 else occli:=joinsublists n2$ 1058 % freelist must not contain substituted variables 1059 freelist:=setdiff(freelist,lhslist); 1060 % It must be possible to separate wrt freelist variables 1061 for each n4 in freelist do 1062 if not freeof(depl!*,n4) then freelist:=delete(n4,freelist); 1063 if freelist then << 1064 n:=nth(eqordr,n1); % order of this equation 1065 h:=simp 0; 1066 for each e1 in xilist do 1067 if (cadddr e1) and ((length cadddr e1) > n) then 1068 % xi (=car e1) is of order n 1069 h:=addsq(h, 1070 if sqzerop car e1 then simp 0 1071 else <<n3:=mergedepli(n3,cadddr e1); 1072 multsq(car e1,simpdf {eqn,cadr e1}) 1073 >> 1074 ); 1075 for each e2 in occli do 1076 h:=addsq(h, 1077 multsq(<<n5:=combidif(e2); 1078 n4:=car etamn(car n5,cdr n5,xilist,etalist, 1079 nth(eqordr,n1),truesub,subordinc,xlist); 1080 vl:=mergedepli(vl,cdr n4); 1081 car n4 1082 >>, 1083 simpdf {eqn,e2} 1084 ) 1085 ); 1086 vl:=joinsublists(xlist . vl)$ 1087 1088 if n1>1 then 1089 for each e1 in reverse flistorg do 1090 if not freeof(symcon,e1) then flist:=union({e1},flist); 1091 1092 for each e2 in freelist do 1093 <<e1:=((car diffsq(h,e2)) . 1); 1094 for n2:=1:eqlen do 1095 <<n4:=nth(lhslist,n2); 1096 if not my_freeof(eqn,n4) then 1097 e1:=((car subsq(e1,{(n4 . cadr nth(truesub,n2))})) . 1); 1098 vl:=delete(n4,vl) 1099 >>; 1100 1101 % splitting 1102 if car e1 and (car e1 neq 0) then << % ie. e1<>0 1103 1104 e1:=cdr split_simplify({{'list, mk!*sq e1}, 1105 'list . nil, 1106 'list . flist, 1107 'list . vl, nil 1108 })$ 1109 symcon:=nconc(e1,symcon) 1110 >> 1111 >>; 1112 if symcon and (individual_ or (n1=eqlen)) then << 1113 h:=callcrack(!*time,cpu,gc,lietrace_,symcon, 1114 flist,vl,xilist,etalist,inequ,nil); 1115 symcon :=car h; xilist:=cadr h; 1116 etalist:=caddr h; flist :=cadddr h; 1117 inequ :=cadddr cdr h; h:=nil$ 1118 cpu:=time()$ gc:=gctime()$ 1119 >> 1120 >> 1121 >>; 1122 1123 %------------ Determining the full symmetry conditions 1124 n1:=0; 1125 vl:=nil; 1126 for each eqn in eqlist do 1127 <<n1:=n1+1; 1128 if !*time then <<terpri()$ 1129 terpri()$terpri()$ 1130 write"=============== Full conditions for the ",n1,". equation" ; 1131 >>; 1132 n2:=cadr adddepli(eqn,revdylist); 1133 n3:=n2; % n3 are the variables in the new condition 1134 1135 h:=simp 0; 1136 for each e1 in xilist do 1137 h:=addsq(h,if sqzerop car e1 then simp 0 1138 else <<n3:=mergedepli(n3,cadddr e1); 1139 multsq(car e1,simpdf {eqn,cadr e1}) 1140 >>); 1141 for each e1 in n2 do 1142 for each e2 in e1 do 1143 h:=addsq(h,multsq(<< %<<<<<<----- slow? 1144 n5:=combidif(e2); 1145 n4:=car etamn(car n5,cdr n5,xilist,etalist,0, 1146 truesub,0,xlist); 1147 n3:=mergedepli(n3,cdr n4); 1148 car n4 1149 >>, 1150 simpdf {eqn,e2} %<<<<<<----- slow 1151 ) 1152 ); 1153 h:=(car h . 1); % ie. take numerator 1154 1155 n3:=joinsublists(xlist . n3)$ 1156 for n2:=1:eqlen do 1157 <<n4:=nth(lhslist,n2); 1158 if not my_freeof(eqn,n4) then 1159 h:=car subsq(h,{(n4 . cadr nth(truesub,n2))}) . 1; % take num 1160 n3:=delete(n4,n3) 1161 >>; 1162 vl:=union(vl,n3); 1163 1164 if prelim_ or (n1>1) then 1165 for each e1 in reverse flistorg do 1166 if not freeof(symcon,e1) then flist:=union({e1},flist); 1167 1168 % splitting 1169 if car h and (car h neq 0) then << 1170 h:=cdr split_simplify({{'list,mk!*sq h}, 1171 'list . nil, 1172 'list . flist, 1173 'list . vl, nil 1174 })$ 1175 symcon:=nconc(h,symcon) 1176 >>$ 1177 1178 last_call:=if n1=eqlen then t else nil$ 1179 if (individual_ and lin_problem) or last_call then << 1180 if last_call then << 1181 etamn_al:=nil$ 1182 if prelim_ or individual_ then << 1183 proc_list_:=proc_list_bak$ 1184 max_gc_fac:=max_gc_fac_bak 1185 >> 1186 >>$ 1187 allsym:=callcrack(!*time,cpu,gc,lietrace_,symcon, 1188 flist,vl,xilist,etalist,inequ,last_call); 1189 cpu:=time()$ gc:=gctime()$ 1190 if last_call then flist:=flistorg 1191 else << 1192 symcon :=car allsym; xilist:=cadr allsym; 1193 etalist:=caddr allsym; flist :=cadddr allsym; 1194 inequ :=cadddr cdr allsym; allsym:=nil$ 1195 >> 1196 >> 1197 >>; 1198 1199 eqn:=sb:=e1:=e2:=n:=h:=dylist:=deplist:=symord:=nil;%occurlist:=nil; 1200 lisp 1201 if done_trafo and cdr done_trafo then << 1202 terpri()$ 1203 if cddr done_trafo then 1204 write"The following transformations reverse the transformations" else 1205 write"The following transformation reverses the transformation"$ 1206 terpri()$ 1207 write"performed in the computation:"$ 1208 algebraic write lisp done_trafo$ 1209 if inverse_trafo_list_incomplete then 1210 write"***** The list 'done_trafo' of inverse transformations"$terpri()$ 1211 write" is not complete as at least one transformation"$terpri()$ 1212 write" could not be inverted"$terpri() 1213 >>$ 1214 1215 n1:=0; 1216 % allsym is a list of solutions of crack 1217 if allsym and car allsym='list then allsym:=cdr allsym; 1218 1219 % Sort the symmetries from most general to most special by 1220 % picking the most special first 1221 h:=for each e1 in allsym collect 1222 cons((length cdadr e1)+(length cdaddr e1),e1); 1223 h:=idx_sort h; 1224 allsym:=for each e1 in h collect cdr e1; 1225 while allsym do << 1226 1227 nodependlist ylist; 1228 symcon_s:=cdadar allsym; 1229 1230 xilist_s :=xilist; 1231 etalist_s:=etalist; 1232 inequ_s :=inequ; 1233 truesub_s:=truesub; 1234 paralist :=nil; 1235 for each g in cdaddr car allsym do 1236 if not freeof(eqlist,cadr g) then << 1237 truesub_s :=subst(caddr g, cadr g, truesub_s); 1238 paralist:=g . paralist 1239 >> else << 1240 xilist_s :=for each k in xilist_s collect cons(subsq(car k,{(cadr g . caddr g)}),cdr k)$ %!! 1241 etalist_s:=for each k in etalist_s collect cons(subsq(car k,{(cadr g . caddr g)}),cdr k)$ %!! 1242 inequ_s :=subst(caddr g, cadr g, inequ_s); 1243 >>; 1244 1245 if lietrace_ then << 1246 write"final symcon : ", symcon_s; terpri()$ 1247 write"final xilist = ", xilist_s; terpri()$ 1248 write"final etalist= ", etalist_s; terpri()$ 1249 >>; 1250 flist_s:=cdr reval cadddr car allsym; 1251 allsym:=cdr allsym$ 1252 oldprint_:=print_; print_:=nil; 1253 if print_ then <<terpri()$ 1254 write"Remaining free functions after the last CRACK-run:"; 1255 terpri()$ 1256 fctprint flist_s;terpri()$terpri() 1257 >>; 1258 1259 %------- Calculation finished, simplification of the result 1260 h:=append(for each el in xilist_s collect mk!*sq car el, %!! 1261 for each el in etalist_s collect mk!*sq car el ); %!! 1262 if symcon_s then for each el in symcon_s do h:=cons(el,h); 1263 h:=cons('list,h); 1264 1265 %------- droping redundant constants or functions 1266 flist_slin:=nil; % flist_slin are the new lin. constants and functions 1267 flist_snli:=nil; % flist_snli are the non-lin. parameters in equations 1268 1269 for each e1 in flist_s do 1270 if freeof(paralist,e1) then flist_slin:=e1 . flist_slin 1271 else flist_snli:=e1 . flist_snli; 1272 1273 oldbatch_mode:=!*batch_mode$ !*batch_mode:=nil$ 1274 if print_ then << 1275 write"***** START OF A COMPUTATION TO DROP REDUNDANT *****"$terpri()$ 1276 write"***** CONSTANTS AND FUNCTIONS OF INTEGRATION *****"$terpri()$ 1277 >>$ 1278 sb:=reval dropredundant(h,'list . flist_slin,'list . vl,list('list)); 1279 if print_ then << 1280 write"***** THE COMPUTATION TO DROP REDUNDANT CONSTANTS *****"$terpri()$ 1281 write"***** AND FUNCTIONS OF INTEGRATION FINISHED *****"$terpri()$ 1282 >>$ 1283 !*batch_mode:=oldbatch_mode$ 1284 1285 if sb then << 1286 flist_slin:=cdr cadddr sb; 1287 h:=caddr sb; 1288 sb:=cadr sb; 1289 e1:=nil 1290 >>; 1291 1292 %------- absorbing numerical constants into free constants 1293 if (not freeoflist( xilist_s,flist)) or 1294 (not freeoflist(etalist_s,flist)) then h:=nil else 1295 h:=reval absorbconst(h,'list . flist_slin); 1296 if h then if sb then sb:=append(sb,cdr h) 1297 else sb:='list . cdr h; 1298 1299 %------- doing the substitutions to drop and absorb 1300 if sb then << 1301 if print_ then <<terpri()$ 1302 write"Free constants and/or functions have been rescaled. ">>$ 1303 for each e1 in cdr sb do << 1304 xilist_s :=for each k in xilist_s collect cons(subsq(car k,{(reval cadr e1 . caddr e1)}),cdr k)$ %!! 1305 etalist_s:=for each k in etalist_s collect cons(subsq(car k,{(reval cadr e1 . caddr e1)}),cdr k)$ %!! 1306 symcon_s :=cdr reval cons('list,subst(caddr e1,reval cadr e1,symcon_s)); 1307 >>; 1308 >>; 1309 1310 symcon_s := for each e1 in symcon_s collect 1311 car simplifypdeSQ(simp e1,append(flist_slin,flist_snli),t,nil,nil)$ 1312 print_:=oldprint_; 1313 1314 %------- Computing the prolongation of the symmetry vector 1315 etapqlist:=nil$ 1316 if fixp(prolong_order) and (prolong_order>0) then << 1317 % We can not do "for each e1 in ylist do depnd(e1,list(xlist));" 1318 % because otherwise the formulas generated by etamn are wrong 1319 % on evallhseqp; % left hand sides not needed 1320 sb:=subdif1(cons('list,xlist),cons('list,ylist),prolong_order)$ 1321 % sb is a list of substitutions, like df(y,x) = y`1 1322 % but with lhs=0 because y is x-independent 1323 1324 for each e1 in cdr sb do 1325 for each e2 in cdr e1 do << 1326 h:=combidif(caddr e2); 1327 n4:=mkid('eta_,car h); 1328 for each n2 in cdr h do 1329 n4:=mkid(n4,nth(xlist,n2)); 1330 h:=car etamn(car h,cdr h,xilist_s,etalist_s,0,truesub_s,0,xlist)$ 1331 n3:=(length cdr h) - 1; 1332 if n3>jetord then jetord:=n3$ 1333 etapqlist:=cons(list('equal,n4,mk!*sq car h),etapqlist); 1334 >> 1335 >>$ 1336 revdylist:=nil; 1337 1338 %------- output 1339 if length flist_slin>1 then n:=t 1340 else n:=nil; 1341 xilist_s:=for each el in xilist_s collect 1342 <<e1:=mkid( 'xi_,cadr el); 1343 e1:=list('equal,e1,prepsq car el); %!! 1344 e1>>; 1345 etalist_s:=for each el in etalist_s collect 1346 <<e1:=mkid('eta_,cadr el); 1347 e1:=list('equal,e1,prepsq car el); %!! 1348 e1>>; 1349 %--- Backsubstitution of e.g. u`1`1 --> df(u,x,2) 1350 for each e1 in ylist do depnd(e1,list(xlist)); 1351 on evallhseqp; 1352 sb:=subdif1(cons('list,xlist),cons('list,ylist),jetord)$ 1353 algebraic( 1354 sb:=for each e1 in sb join 1355 for each e2 in e1 collect(rhs e2 = lhs e2)); 1356 off evallhseqp; 1357 xilist_s :=cdr algebraic(sub(sb,cons('list, xilist_s))); 1358 etalist_s:=cdr algebraic(sub(sb,cons('list,etalist_s))); 1359 etapqlist:=cdr algebraic(sub(sb,cons('list,etapqlist))); 1360 xicop := xilist_s$ 1361 etacop :=etalist_s$ 1362 etapqcop:=etapqlist$ 1363 1364 sb:=nil$ 1365 flcop:=flist_slin; 1366 % Drop constants/functions of integration which may have come up 1367 % when integrating parametric functions in the equation: 1368 for each e1 in flist_slin do 1369 if not freeof(eqlist,e1) then flcop:=delete(e1,flcop); 1370 for each e1 in flcop do 1371 <<% if null n2 then n2:=freeof(xicop,e1) and freeof(etacop,e1)$ 1372 if freeof(symcon_s,e1) then << 1373 n1:=n1+1; 1374 xi:=xicop;eta:=etacop;etapq:=etapqcop; 1375 for each e2 in flcop do 1376 if e2 neq e1 then 1377 <<xi:=subst(0,e2,xi);eta:=subst(0,e2,eta);etapq:=subst(0,e2,etapq)>> 1378 else 1379 if null cdr fargs e1 then 1380 <<xi:=subst(1,e2,xi);eta:=subst(1,e2,eta);etapq:=subst(1,e2,etapq)>>; 1381 terpri()$write"-------- ",n1,". Symmetry:";terpri()$ 1382 for each e2 in paralist do algebraic write lisp {'equal,cadr e2,reval caddr e2}; 1383 for each e2 in xi do algebraic write lisp {'equal,cadr e2,reval caddr e2}; 1384 for each e2 in eta do algebraic write lisp {'equal,cadr e2,reval caddr e2}; 1385 for each e2 in etapq do algebraic write lisp {'equal,cadr e2,reval caddr e2}; 1386 if cdr fargs e1 then <<terpri()$write"with ";fctprint list(e1); 1387 terpri()>>$ 1388 xicop :=subst(0,e1,xicop ); 1389 etacop :=subst(0,e1,etacop ); 1390 etapqcop:=subst(0,e1,etapqcop); 1391 flcop:=delete(e1,flcop); 1392 %depl!*:=delete(assoc(e1,depl!*),depl!*)$ 1393 >>; 1394 >>; 1395 1396 if flcop then << 1397 if ((length flcop)+(length flist_snli))>1 then n:=t 1398 else n:=nil; 1399 terpri()$ 1400 if flcop=flist_slin then write"-------- S" 1401 else write"-------- Further s"; 1402 write"ymmetr",if n then "ies:" else "y:"$ 1403 terpri()$ 1404 for each e1 in paralist do algebraic write lisp {'equal,cadr e1,reval caddr e1}; 1405 for each e1 in xicop do algebraic write lisp {'equal,cadr e1,reval caddr e1}; 1406 for each e1 in etacop do algebraic write lisp {'equal,cadr e1,reval caddr e1}; 1407 for each e1 in etapqcop do algebraic write lisp {'equal,cadr e1,reval caddr e1} 1408 >>; 1409 1410 terpri()$ 1411 if flcop then 1412 <<write"with ";fctprint cdr reval ('list . append(flcop,flist_snli))>>; 1413 1414 if null symcon_s then if flcop then 1415 <<write" which ",if n then "are" else "is"," free. "; 1416 terpri()>> else 1417 else 1418 <<h:=print_;print_:=50$ 1419 symcon_s:=for each a in symcon_s collect {'!*sq,a,t}$ 1420 if print_ then 1421 <<terpri()$ 1422 write"which still ",if n then "have" else "has"," to satisfy: "; 1423 terpri()$ 1424 deprint symcon_s; 1425 >> else 1426 <<terpri()$ 1427 write"which ",if n then "have" else "has", 1428 " to satisfy conditions. To see them set "; 1429 terpri()$ 1430 write 1431 "lisp(print_= max. number of terms of an equation to print);"; 1432 terpri()$ 1433 >>; 1434 print_:=h 1435 >>; 1436 return_list:=list('list, 1437 'list . symcon_s, 1438 'list . append(paralist,append(xilist_s,etalist_s)), 1439 'list . append(flist_slin,flist_snli)) . return_list; 1440 nodependlist ylist 1441 >>; 1442 1443 terpri()$ 1444 if (n1=0) and null symcon_s then << 1445 if length eqlist=1 then write"The equation has no symmetry." 1446 else write"The equations have no symmetry."$ 1447 terpri() 1448 >>$ 1449 write"-------- ";terpri(); 1450 1451 for each e1 in ylist do depnd(e1,list(xlist)); 1452 recover_reduce_flags()$ 1453 collect_sol:=oldcollect_sol; 1454 adjust_fnc:=oldadj; 1455 1456 return 'list . return_list 1457end$ % of liepde 1458 1459endmodule$ 1460 1461end$ 1462