1% Raffaele Vitolo, 09/10/09 2% This is the computation for (higher) symmetries of Burgers 3 4load_package(cdiff); 5 6 7 8% The following instructions initialize the total derivatives. The first 9% string is the name of the vector field, 10% the second item is the list of even variables 11% (note that u1, u2, ... are u_x, u_xx, ...), 12% the third item is the list of non-commuting variables 13% (`ext' stands for `external' like in external (wedge) product). 14 15super_vectorfield(ddx,{x,t,u,u1,u2,u3,u4,u5,u6,u7, 16u8,u9,u10,u11,u12,u13,u14,u15,u16,u17}, 17{ext 1,ext 2,ext 3,ext 4,ext 5,ext 6,ext 7,ext 8,ext 9,ext 10,ext 1811,ext 12,ext 13,ext 14,ext 15,ext 16,ext 17,ext 18,ext 19,ext 20,ext 1921,ext 22,ext 23,ext 24,ext 25,ext 26,ext 27,ext 28,ext 29,ext 30, 20ext 31,ext 32,ext 33,ext 34,ext 35,ext 36,ext 37,ext 38,ext 39,ext 40, 21ext 41,ext 42,ext 43,ext 44,ext 45,ext 46,ext 47,ext 48,ext 49,ext 50, 22ext 51,ext 52,ext 53,ext 54,ext 55,ext 56,ext 57,ext 58,ext 59,ext 60, 23ext 61,ext 62,ext 63,ext 64,ext 65,ext 66,ext 67,ext 68,ext 69,ext 70, 24ext 71,ext 72,ext 73,ext 74,ext 75,ext 76,ext 77,ext 78,ext 79,ext 80 25}); 26 27 28{20,80} 29 30 31super_vectorfield(ddt,{x,t,u,u1,u2,u3,u4,u5,u6,u7, 32u8,u9,u10,u11,u12,u13,u14,u15,u16,u17}, 33{ext 1,ext 2,ext 3,ext 4,ext 5,ext 6,ext 7,ext 8,ext 9,ext 10,ext 3411,ext 12,ext 13,ext 14,ext 15,ext 16,ext 17,ext 18,ext 19,ext 20,ext 3521,ext 22,ext 23,ext 24,ext 25,ext 26,ext 27,ext 28,ext 29,ext 30, 36ext 31,ext 32,ext 33,ext 34,ext 35,ext 36,ext 37,ext 38,ext 39,ext 40, 37ext 41,ext 42,ext 43,ext 44,ext 45,ext 46,ext 47,ext 48,ext 49,ext 50, 38ext 51,ext 52,ext 53,ext 54,ext 55,ext 56,ext 57,ext 58,ext 59,ext 60, 39ext 61,ext 62,ext 63,ext 64,ext 65,ext 66,ext 67,ext 68,ext 69,ext 70, 40ext 71,ext 72,ext 73,ext 74,ext 75,ext 76,ext 77,ext 78,ext 79,ext 80 41}); 42 43 44{20,80} 45 46 47 48% Specification of the vectorfield ddx. 49% The meaning of the first index is the parity of variables. 50% In particular here we have just even variables. 51% The second index parametrizes the second item (list) 52% in the super_vectorfield declaration. 53 54ddx(0,1):=1$ 55 56 57ddx(0,2):=0$ 58 59 60ddx(0,3):=u1$ 61 62 63ddx(0,4):=u2$ 64 65 66ddx(0,5):=u3$ 67 68 69ddx(0,6):=u4$ 70 71 72ddx(0,7):=u5$ 73 74 75ddx(0,8):=u6$ 76 77 78ddx(0,9):=u7$ 79 80 81ddx(0,10):=u8$ 82 83 84ddx(0,11):=u9$ 85 86 87ddx(0,12):=u10$ 88 89 90ddx(0,13):=u11$ 91 92 93ddx(0,14):=u12$ 94 95 96ddx(0,15):=u13$ 97 98 99ddx(0,16):=u14$ 100 101 102ddx(0,17):=u15$ 103 104 105ddx(0,18):=u16$ 106 107 108ddx(0,19):=u17$ 109 110 111ddx(0,20):=letop$ 112 113 114 115% Specification of the vectorfield ddt 116% In the evolutionary case we never have more than one time derivative 117% other derivatives are u_txxx ... 118 119ddt(0,1):=0$ 120 121 122ddt(0,2):=1$ 123 124 125ddt(0,3):=ut$ 126 127 128ddt(0,4):=ut1$ 129 130 131ddt(0,5):=ut2$ 132 133 134ddt(0,6):=ut3$ 135 136 137ddt(0,7):=ut4$ 138 139 140ddt(0,8):=ut5$ 141 142 143ddt(0,9):=ut6$ 144 145 146ddt(0,10):=ut7$ 147 148 149ddt(0,11):=ut8$ 150 151 152ddt(0,12):=ut9$ 153 154 155ddt(0,13):=ut10$ 156 157 158ddt(0,14):=ut11$ 159 160 161ddt(0,15):=ut12$ 162 163 164ddt(0,16):=ut13$ 165 166 167ddt(0,17):=ut14$ 168 169 170ddt(0,18):=letop$ 171 172 173ddt(0,19):=letop$ 174 175 176ddt(0,20):=letop$ 177 178 179 180% The equation -- this is also used to specify internal variables. 181% For evolutionary equations internal variables are of the type 182% (t,x,u,u_x,u_xx,...) 183 184ut:=u2+2*u*u1; 185 186 187ut := 2*u*u1 + u2 188 189 190ut1:=ddx ut; 191 192 193 2 194ut1 := 2*u*u2 + 2*u1 + u3 195 196ut2:=ddx ut1; 197 198 199ut2 := 2*u*u3 + 6*u1*u2 + u4 200 201ut3:=ddx ut2; 202 203 204 2 205ut3 := 2*u*u4 + 8*u1*u3 + 6*u2 + u5 206 207ut4:=ddx ut3; 208 209 210ut4 := 2*u*u5 + 10*u1*u4 + 20*u2*u3 + u6 211 212ut5:=ddx ut4; 213 214 215 2 216ut5 := 2*u*u6 + 12*u1*u5 + 30*u2*u4 + 20*u3 + u7 217 218ut6:=ddx ut5; 219 220 221ut6 := 2*u*u7 + 14*u1*u6 + 42*u2*u5 + 70*u3*u4 + u8 222 223ut7:=ddx ut6; 224 225 226 2 227ut7 := 2*u*u8 + 16*u1*u7 + 56*u2*u6 + 112*u3*u5 + 70*u4 + u9 228 229ut8:=ddx ut7; 230 231 232ut8 := 2*u*u9 + 18*u1*u8 + u10 + 72*u2*u7 + 168*u3*u6 + 252*u4*u5 233 234ut9:=ddx ut8; 235 236 237 2 238ut9 := 2*u*u10 + 20*u1*u9 + u11 + 90*u2*u8 + 240*u3*u7 + 420*u4*u6 + 252*u5 239 240ut10:=ddx ut9; 241 242 243ut10 := 244 2452*u*u11 + 22*u1*u10 + u12 + 110*u2*u9 + 330*u3*u8 + 660*u4*u7 + 924*u5*u6 246 247ut11:=ddx ut10; 248 249 250ut11 := 2*u*u12 + 24*u1*u11 + 132*u10*u2 + u13 + 440*u3*u9 + 990*u4*u8 251 252 2 253 + 1584*u5*u7 + 924*u6 254 255ut12:=ddx ut11; 256 257 258ut12 := 2*u*u13 + 26*u1*u12 + 572*u10*u3 + 156*u11*u2 + u14 + 1430*u4*u9 259 260 + 2574*u5*u8 + 3432*u6*u7 261 262ut13:=ddx ut12; 263 264 265ut13 := 2*u*u14 + 28*u1*u13 + 2002*u10*u4 + 728*u11*u3 + 182*u12*u2 + u15 266 267 2 268 + 4004*u5*u9 + 6006*u6*u8 + 3432*u7 269 270ut14:=ddx ut13; 271 272 273ut14 := 2*u*u15 + 30*u1*u14 + 6006*u10*u5 + 2730*u11*u4 + 910*u12*u3 274 275 + 210*u13*u2 + u16 + 10010*u6*u9 + 12870*u7*u8 276 277 278% Test for verifying the commutation of total derivatives. 279% Highest order defined terms yield some `letop' 280% which means `careful' in Dutch and is treated as a new variable. 281operator ev; 282 283 284 285for i:=1:17 do write ev(0,i):=ddt(ddx(0,i))-ddx(ddt(0,i)); 286 287 288ev(0,1) := 0 289 290ev(0,2) := 0 291 292ev(0,3) := 0 293 294ev(0,4) := 0 295 296ev(0,5) := 0 297 298ev(0,6) := 0 299 300ev(0,7) := 0 301 302ev(0,8) := 0 303 304ev(0,9) := 0 305 306ev(0,10) := 0 307 308ev(0,11) := 0 309 310ev(0,12) := 0 311 312ev(0,13) := 0 313 314ev(0,14) := 0 315 316ev(0,15) := 0 317 318ev(0,16) := 0 319 320ev(0,17) := letop - 2*u*u16 - 32*u1*u15 - 16016*u10*u6 - 8736*u11*u5 321 322 - 3640*u12*u4 - 1120*u13*u3 - 240*u14*u2 - u17 - 22880*u7*u9 323 324 2 325 - 12870*u8 326 327 328%% This is the list of variables with respect to their grading, 329%% starting from degree ONE. 330 331graadlijst:={{u},{u1},{u2},{u3},{u4},{u5}, 332{u6},{u7},{u8},{u9},{u10},{u11},{u12},{u13},{u14},{u15},{u16},{u17}}; 333 334 335graadlijst := {{u}, 336 337 {u1}, 338 339 {u2}, 340 341 {u3}, 342 343 {u4}, 344 345 {u5}, 346 347 {u6}, 348 349 {u7}, 350 351 {u8}, 352 353 {u9}, 354 355 {u10}, 356 357 {u11}, 358 359 {u12}, 360 361 {u13}, 362 363 {u14}, 364 365 {u15}, 366 367 {u16}, 368 369 {u17}} 370 371 372% This is the list of all monomials of degree 0, 1, 2, ... 373% which can be constructed from the above list of elementary variables 374% with their grading. 375 376grd0:={1}; 377 378 379grd0 := {1} 380 381grd1:= mkvarlist1(1,1)$ 382 383 384grd2:= mkvarlist1(2,2)$ 385 386 387grd3:= mkvarlist1(3,3)$ 388 389 390grd4:= mkvarlist1(4,4)$ 391 392 393grd5:= mkvarlist1(5,5)$ 394 395 396grd6:= mkvarlist1(6,6)$ 397 398 399grd7:= mkvarlist1(7,7)$ 400 401 402grd8:= mkvarlist1(8,8)$ 403 404 405grd9:= mkvarlist1(9,9)$ 406 407 408grd10:= mkvarlist1(10,10)$ 409 410 411grd11:= mkvarlist1(11,11)$ 412 413 414grd12:= mkvarlist1(12,12)$ 415 416 417grd13:= mkvarlist1(13,13)$ 418 419 420grd14:= mkvarlist1(14,14)$ 421 422 423grd15:= mkvarlist1(15,15)$ 424 425 426grd16:= mkvarlist1(16,16)$ 427 428 429 430% Initialize a counter for the vector of arbitrary constants 431operator c,equ; 432 433 434 435ctel:=0; 436 437 438ctel := 0 439 440 441% we assume a generating function of degree <= 5 442 443sym:= 444(for each el in grd0 sum (c(ctel:=ctel+1)*el))+ 445(for each el in grd1 sum (c(ctel:=ctel+1)*el))+ 446(for each el in grd2 sum (c(ctel:=ctel+1)*el))+ 447(for each el in grd3 sum (c(ctel:=ctel+1)*el))+ 448(for each el in grd4 sum (c(ctel:=ctel+1)*el))+ 449(for each el in grd5 sum (c(ctel:=ctel+1)*el))$ 450 451 452 453% This is the equation ell_B(sym)=0, where B=0 is Burgers'equation 454% and sym is the generating function. From now on all equations 455% are arranged in a single vector whose name is `equ'. 456 457equ 1:=ddt(sym)-ddx(ddx(sym))-2*u*ddx(sym)-2*u1*sym ; 458 459 460 2 461equ(1) := 2*(4*c(19)*u1*u4 + 10*c(19)*u2*u3 + 3*c(18)*u*u1*u3 + 3*c(18)*u*u2 462 463 2 2 464 - c(18)*u1*u4 + 3*c(17)*u1 *u2 - c(17)*u2*u3 + 2*c(16)*u *u1*u2 465 466 2 3 2 467 - 2*c(16)*u*u1*u3 - c(16)*u1 *u2 + c(15)*u*u1 - c(15)*u*u2 468 469 2 2 3 5 470 - 2*c(15)*u1 *u2 - 3*c(14)*u *u1*u2 - 3*c(14)*u*u1 - c(13)*u *u1 471 472 3 2 2 473 - 10*c(13)*u *u1 + 3*c(12)*u1*u3 + 3*c(12)*u2 + 2*c(11)*u*u1*u2 474 475 3 2 3 476 - c(11)*u1*u3 + c(10)*u1 - c(10)*u2 - 2*c(9)*u*u1*u2 - c(9)*u1 477 478 4 2 2 479 - c(8)*u *u1 - 6*c(8)*u *u1 + 2*c(7)*u1*u2 - c(6)*u1*u2 480 481 3 2 2 2 482 - c(5)*u *u1 - 3*c(5)*u*u1 - c(3)*u *u1 - c(3)*u1 - c(2)*u*u1 483 484 - c(1)*u1) 485 486 487% This is the list of variables, to be passed to the equation solver. 488 489vars:={x,t,u,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,u15,u16,u17}; 490 491 492vars := {x, 493 494 t, 495 496 u, 497 498 u1, 499 500 u2, 501 502 u3, 503 504 u4, 505 506 u5, 507 508 u6, 509 510 u7, 511 512 u8, 513 514 u9, 515 516 u10, 517 518 u11, 519 520 u12, 521 522 u13, 523 524 u14, 525 526 u15, 527 528 u16, 529 530 u17} 531 532 533% This is the number of initial equation(s) 534 535tel:=1; 536 537 538tel := 1 539 540 541% The following procedure uses multi_coeff (from the package `tools'). 542% It gets all coefficients of monomials appearing in the initial equation(s). 543% The coefficients are put into the vector equ after the initial equations. 544 545procedure splitvars i; 546begin; 547ll:=multi_coeff(equ i,vars); 548equ(tel:=tel+1):=first ll; 549ll:=rest ll; 550for each el in ll do equ(tel:=tel+1):=second el; 551end; 552 553 554splitvars 555 556 557% This command initialize the equation solver. 558% It passes the equation(s) togeher with their number `tel', 559% the constants'vector `c', its length `ctel', 560% an arbitrary constant `f' that may appear in computations. 561 562initialize_equations(equ,tel,{},{c,ctel,0},{f,0,0}); 563 564 565 566% Run the procedure splitvars in order to obtain equations on coefficiens 567% of each monomial. 568 569splitvars 1; 570 571 572 573% Next command tells the solver the total number of equations obtained 574% after running splitvars. 575 576put_equations_used tel; 577 578 57924 580 581 582% It is worth to write down the equations for the coefficients. 583 584for i:=2:tel do write equ i; 585 586 5870 588 589 - 2*c(1) 590 591 - 2*c(2) 592 593 - 2*c(3) 594 595 - 2*c(3) 596 597 - 6*c(5) 598 599 - 2*c(5) 600 6012*(2*c(7) - c(6)) 602 603 - 12*c(8) 604 605 - 2*c(8) 606 6072*(c(10) - c(9)) 608 6094*(c(11) - c(9)) 610 6112*(3*c(12) - c(10)) 612 6132*(3*c(12) - c(11)) 614 615 - 20*c(13) 616 617 - 2*c(13) 618 6192*(c(15) - 3*c(14)) 620 6212*(2*c(16) - 3*c(14)) 622 6232*(3*c(17) - c(16) - 2*c(15)) 624 6252*(3*c(18) - c(15)) 626 6272*(3*c(18) - 2*c(16)) 628 6292*(10*c(19) - c(17)) 630 6312*(4*c(19) - c(18)) 632 633 634% This command solves the equations for the coefficients. 635% Note that we have to skip the initial equations! 636 637for i:=2:tel do integrate_equation i; 638 639 640equ(2) = 0 641 642equ(3): Solved for c(1) 643 644equ(4): Solved for c(2) 645 646equ(5): Solved for c(3) 647 648equ(6) = 0 649 650equ(7): Solved for c(5) 651 652equ(8) = 0 653 654equ(9): Solved for c(7) 655 656equ(10): Solved for c(8) 657 658equ(11) = 0 659 660equ(12): Solved for c(10) 661 662equ(13): Solved for c(11) 663 664equ(14): Solved for c(12) 665 666equ(15) = 0 667 668equ(16): Solved for c(13) 669 670equ(17) = 0 671 672equ(18): Solved for c(15) 673 674equ(19): Solved for c(16) 675 676equ(20): Solved for c(17) 677 678equ(21): Solved for c(18) 679 680equ(22) = 0 681 682equ(23): Solved for c(19) 683 684equ(24) = 0 685 686 687write "sym:=",sym; 688 689 690 3 2 2 691sym:=(12*c(14)*u *u1 + 18*c(14)*u *u2 + 36*c(14)*u*u1 + 12*c(14)*u*u3 692 693 2 694 + 30*c(14)*u1*u2 + 3*c(14)*u4 + 12*c(9)*u *u1 + 12*c(9)*u*u2 695 696 2 697 + 12*c(9)*u1 + 4*c(9)*u3 + 12*c(6)*u*u1 + 6*c(6)*u2 + 12*c(4)*u1)/12 698 699 700; 701 702end; 703 704Tested on x86_64-pc-windows CSL 705Time (counter 1): 16 ms 706 707End of Lisp run after 0.01+0.06 seconds 708real 0.21 709user 0.00 710sys 0.07 711