MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post
  • Latest Posts Feed
  • I'm particularly interested in data analysis and more specifically in statistical analysis of computer code outputs.

    One of the main activity of this very broad field is named Uncertainty Propagation. In a few words it consists in perturbing the inputs of a computational code in order to understand (and quantify) how these perturbations propagates through the outputs of this code.

    At the core of uncertainty propagation is the ability to generate large numbers of "random" variations of the inputs. Knowing that these entries can be counted in tens, one sees that the first problem consists in generating "random" points in a space of potentially very large dimension.

    Even among my mathematician colleagues an impressive number of them is completely ignorant of the way "random" numbers are generated. I guess that a lot of Mapleprimes' users are too. My purpose is not to give a course on this topic and the affording litterature is vast enough for everyone interested might find informations of any level of complexity.
    Among those who have some knowledge about Pseudo Random Numbers Generators (PRNG), only a few of them know that a PRNG has to pass severe tests ("tests of randomness") before the streams of number it generates might  be qualified as "reasonably random" and therefore this PRNG might be released.

    One of most famous example of a bad PRNG is given by "randu" (IBM 1966, and probably used in Fortran libraries during more than 30 years), this same PRNG that Knuth qualified himself as the "infamous generator".

    These tests of randomness are generally gathered in dedicated libraries and Diehard is probably tone of the most known of them.
    Diehard has originally been developed by George Marsaglia more than twenty years ago and it's still widely ued today.

    I recently decided, not because I have doubts about the quality of the work done by Maplesoft, to test the Maple's PRNG named "Mersenne Twister". First, because it can do no harm to publish quantitative information that allows everyone to know that it is using a proven PRNG; second, because the (very simple) approach used here can fill the gaps I have mentioned above.

    Mersenne Twister (often dubbed mt19937) is considered as a very good PRNG; it is used in a lot of applications (including finance where it is not so rare to sample input spaces of dimensions larger than 1000... ok I know, mt19937 is often considered as a poor candidate for cryptography applications, but it's not my concern here).

    I have thus decided to spend some time to run the Diehard suite of tests on a sequence of integers numbers generated by RandomTools[MersenneTwister].


     

    restart:


    DIEHARD tests suite for Pseudo Random Numbers Generators (PRNG)

    Reference: http://webhome.phy.duke.edu/~rgb/General/dieharder.php

    The installation procedure (Mac OSX) can be found here
        https://gist.github.com/blixt/9abfafdd0ada0f4f6f26
    or here
        http://macappstore.org/dieharder/

    For other operating systems, please search on the web pages.


    dieharder [-h]   # for inline help
    dieharder -l      # to get the lists all the avaliable tests




    A description of the many tests can be found here:
        https://en.wikipedia.org/wiki/Diehard_tests
        https://sites.google.com/site/astudyofentropy/background-information/the-tests/dieharder-test-descriptions
        https://www.stata.com/support/cert/diehard/randnumb_mt64.out

    General theory about PRNG testing can be found here (a reference among many):
        http://liu.diva-portal.org/smash/get/diva2:740158/FULLTEXT01.pdf

    or here (more oriented to the NIST test suite)
        https://www.random.org/analysis/Analysis2005.pdf
        https://nvlpubs.nist.gov/nistpubs/legacy/sp/nistspecialpublication800-22r1a.pdf



    In a terminal window execute the following commands for an exhaustive testing ("-a" option).
    The "-g 202" option means that the generator is replaced by a text format input file
    (use dieharder -h for more details).

    cd //..../Desktop/DIEHARD

    dieharder -g 202 -f SomeAsciiFile -a > //..../Desktop/DIEHARD/TheResultFile.txt

    Be carefull, the complete testing takes several hours (about 5 on my computer)



    __________________________________________________________________________________
     


    Maple's Mersenne Twister Generator

    Maple help page : RandomTools[MersenneTwister][GenerateInteger]
    (see rincluded references to the Mersenne Twister PRNG).

    Note: in the sequel this generator will be dubbed mt19937


    The Mersenne Twister is implemented in many softwares.
    It is higly likely that this PRNG (and the others these softwares propose) have been intensively
    tested with one of the existing PRNG testing libraries.
    Unfortunately only a few editors have made public the results of these tests (probably because
    the implementation in itself is rarely questioned... but a code typo is always a possibility).

    One exception is ths software STATA.
    A summary of the results can be found here
       https://www.stata.com/support/cert/diehard/.
    A complete description of the results of the tests passed is given here
       https://www.stata.com/support/cert/diehard/randnumb_mt64.out

    The classical pattern of the performances of mt19937 can be found here

       http://www2.ic.uff.br/~celso/artigos/pjo6.ps.

    and the table below comes from it (P means "Passed", F means "Failed"):


    ____________________________________________________________________________


    In the Maple code below, a sequence of N UnsignedInt32 numbers is generated from the
    Maple's Mersenne Twister and the result is exported in an ASCII file.
    The Seed is set to 1 (SetState(state=1)) to compare, with a small value of N (let's say N=10)
    the sequence produced by Maple's mt19937 with the the sequence of the same length generated
    by Diehard's mt19937.
    To generate this later sequence and save it in file Diehard_mt19937, just run in a terminan window
    the command (-S 1 means "seed = 1", -t 10 means "a sequence of length 10"):
       dieharder -S 1 -B -o -t 10 > Diehard_mt19937

    About the value of N:

    In http://webhome.phy.duke.edu/~rgb/General/dieharder.php it's recommend that N be at least
    equal to 2.5 million; STATA used N=3 million.
    Other web sources say this value is too small.
    For N=10 million the Maple's mt19937 doesn't pass the tests successfully.
    I used here N=50 million (the resulting ASCII file has size 537 Mo).



    Name of the input file.

    The file generated by Maple is named Maple_mt19937_N=5e7.txt



    One important thing is the preamble of a licit input file.

    This preamble must have 6 lines (the value 10 right to count must be set to the value of N).
    A licit preamble is of the form.

    #==================================================================

    # some text indicating the generator used

    #==================================================================

    type: d

    count: 10

    numbit: 32

    As Maple_mt19937_N=5e7.txt is generated from an ExportMatrix command, this preamble is added
    by hand.
     


    Running multiple Diehard tests

    To run the same tests used to qualify STATA's Mersenne Twister, open a terminal window,
    go to the directory that contains input file Maple_mt19937_N=5e7.txt and run this script:

     for i in {0,1,2,3,4,8,9,10,11,12,13,14,15,16}; do

        dieharder -g 202 -f Maple_mt19937_N=5e7.txt -d $i >> Diehard___Maple_mt19937_N=5e7

     done ;

    The results are then forked in the ASCII file Diehard___Maple_mt19937_N=5e7

     

    with(RandomTools[MersenneTwister]):

    dir := cat("/", currentdir(), "Desktop/DIEHARD/"):
    InputFile := cat(dir, "Maple_mt19937_N=5e7.txt"):

    SetState(state=1);

    N := 5*10^7:

    st := time():
    S := convert([seq(GenerateUnsignedInt32(), i=1..N)], Matrix)^+;
    time()-st;

    S := Vector(4, {(1) = ` 50000000 x 1 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

     

    84.526

    (1)

    st := time():
    ExportMatrix(InputFile, S, format=rectangular, mode=ascii);
    time()-st;

    537066525

     

    61.926

    (2)


    Diehard's results


    Full test suite (about 5 hours of computational time)

    Command :
    dieharder -g 202 -f Maple_mt19937_N=5e7.txt -a > Diehard___ALL___Maple_mt19937_N=5e7


    The results are compared to those obtained for Diehard's mt19937.
    Two ways are used :

      - 1 - In a first stage one generates a stream of PRN and store it in an ASCII file (just as we did with Maple).
             The whole suite of tests is then run on this file.
             Commands (-g 013 codes for mt19937):

             dieharder -S 1 -g 013 -o -t 50000000 > Diehard_mt19937_N=5e7.txt
             dieharder -g 202 -f Diehard_mt19937_N=5e7.txt -a > Diehard___ALL___Diehard_mt19937_N=5e7



      - 2 - The whole suite is run by invoking directectly mt19937 "online"
             Commands :
             dieharder -S 1 -g 013 -t 50000000 -a > Diehard___ALL___Online


    A UNIX diff command has been used to verify that the two files Maple_mt19937_N=5e7.txt and
     Diehard_mt19937_N=5e7.txt were identical (thet were).

    Note that the Diehard doens't responds identically depending on the stream of random numbers comes from a file
    or is generated online (this last [- 2 -] situation seems to give better results).-

    Résumé (114 tests):
       - * - Maple's  and Diehard's  mt19937 respond exactly the same way when the stream of random
              numbers is read from an ASCII file (8 tests failed (******) and 6 weak (**)).
       - * - Diehard's  mt19937 fails 0 test and is weak on 4 tests when the stream is generated online
     

     

    restart:

    dir := currentdir():
    FromMapleFile     := cat(dir, "Diehard___ALL___Maple_mt19937_N=5e7"):
    FromDiehardFile   := cat(dir, "Diehard___ALL___diehard_mt19937_N=5e7"):
    FromDiehardNoFile := cat(dir, "Diehard___ALL___Online"):


    printf("                           ======================|======================|======================|\n"):
    printf("                          |   From Maple's file  | From Diehard's File  | Diehard online test  |\n"):
    printf("==========================|======================|======================|======================|\n"):
    printf("          test       ntup | p.value   Assessment | p.value   Assessment | p.value   Assessment |\n"):
    printf("==========================|======================|======================|======================|\n"):


    for k from 1 to 9 do
      LMF  := readline(FromMapleFile):
      LDF  := readline(FromDiehardFile):
      LDNF := readline(FromDiehardNoFile):
    end do:


    while LMF <> 0 do
      if StringTools:-Search("|", LMF) > 0 then
        res := StringTools:-StringSplit(LMF, "|")[[1, 2, 5, 6]];
        printf("%-20s  %3d | %1.7f ", res[1], parse(res[2]), parse(res[3]));
          if StringTools:-Search("WEAK"  , res[4]) > 0 then printf("    **     |")
        elif StringTools:-Search("FAILED", res[4]) > 0 then printf("  ******   |")
        else printf("  PASSED   |")
        end if:
      end if:
      LMF  := readline(FromMapleFile):

      if StringTools:-Search("|", LDF) > 0 then
        res := StringTools:-StringSplit(LDF, "|")[[5, 6]];
        printf(" %1.7f ", parse(res[1]));
          if StringTools:-Search("  WEAK"  , res[2]) > 0 then printf("     **    |")
        elif StringTools:-Search("  FAILED", res[2]) > 0 then printf("   ******  |")
        else printf("   PASSED  |")
        end if:
      end if:
      LDF  := readline(FromDiehardFile):
                         
      if StringTools:-Search("|", LDNF) > 0 then
        res := StringTools:-StringSplit(LDNF, "|")[[5, 6]];
        printf(" %1.7f ", parse(res[1]));
          if StringTools:-Search("WEAK"  , res[2]) > 0 then printf("     **    |")
        elif StringTools:-Search("FAILED", res[2]) > 0 then printf("   ******    |")
        else printf("   PASSED  |")
        end if:
        printf("\n"):
      end if:
      LDNF := readline(FromDiehardNoFile):


    end do:

                               ======================|======================|======================|
                              |   From Maple's file  | From Diehard's File  | Diehard online test  |
    ==========================|======================|======================|======================|
              test       ntup | p.value   Assessment | p.value   Assessment | p.value   Assessment |
    ==========================|======================|======================|======================|
       diehard_birthdays    0 | 0.9912651   PASSED   | 0.9912651    PASSED  | 0.8284550    PASSED  |
          diehard_operm5    0 | 0.1802226   PASSED   | 0.1802226    PASSED  | 0.5550587    PASSED  |
      diehard_rank_32x32    0 | 0.3099035   PASSED   | 0.3099035    PASSED  | 0.9575440    PASSED  |
        diehard_rank_6x8    0 | 0.2577249   PASSED   | 0.2577249    PASSED  | 0.3915666    PASSED  |
       diehard_bitstream    0 | 0.5519218   PASSED   | 0.5519218    PASSED  | 0.9999462      **    |
            diehard_opso    0 | 0.1456442   PASSED   | 0.1456442    PASSED  | 0.7906533    PASSED  |
            diehard_oqso    0 | 0.4882425   PASSED   | 0.4882425    PASSED  | 0.9574014    PASSED  |
             diehard_dna    0 | 0.0102880   PASSED   | 0.0102880    PASSED  | 0.5149193    PASSED  |
    diehard_count_1s_str    0 | 0.1471956   PASSED   | 0.1471956    PASSED  | 0.9517290    PASSED  |
    diehard_count_1s_byt    0 | 0.1158707   PASSED   | 0.1158707    PASSED  | 0.1568255    PASSED  |
     diehard_parking_lot    0 | 0.1148982   PASSED   | 0.1148982    PASSED  | 0.1611173    PASSED  |
        diehard_2dsphere    2 | 0.9122204   PASSED   | 0.9122204    PASSED  | 0.2056657    PASSED  |
        diehard_3dsphere    3 | 0.9385972   PASSED   | 0.9385972    PASSED  | 0.3620517    PASSED  |
         diehard_squeeze    0 | 0.2686977   PASSED   | 0.2686977    PASSED  | 0.8611266    PASSED  |
            diehard_sums    0 | 0.1602355   PASSED   | 0.1602355    PASSED  | 0.5103248    PASSED  |
            diehard_runs    0 | 0.1235328   PASSED   | 0.1235328    PASSED  | 0.9402086    PASSED  |
            diehard_runs    0 | 0.6341956   PASSED   | 0.6341956    PASSED  | 0.3274267    PASSED  |
           diehard_craps    0 | 0.0243605   PASSED   | 0.0243605    PASSED  | 0.1844482    PASSED  |
           diehard_craps    0 | 0.2952043   PASSED   | 0.2952043    PASSED  | 0.1407422    PASSED  |
     marsaglia_tsang_gcd    0 | 0.0000000   ******   | 0.0000000    ******  | 0.5840531    PASSED  |
     marsaglia_tsang_gcd    0 | 0.0000000   ******   | 0.0000000    ******  | 0.8055035    PASSED  |
             sts_monobit    1 | 0.9397218   PASSED   | 0.9397218    PASSED  | 0.9018886    PASSED  |
                sts_runs    2 | 0.8092469   PASSED   | 0.8092469    PASSED  | 0.2247600    PASSED  |
              sts_serial    1 | 0.2902851   PASSED   | 0.2902851    PASSED  | 0.9223063    PASSED  |
              sts_serial    2 | 0.9541680   PASSED   | 0.9541680    PASSED  | 0.6140772    PASSED  |
              sts_serial    3 | 0.4090798   PASSED   | 0.4090798    PASSED  | 0.2334754    PASSED  |
              sts_serial    3 | 0.5474851   PASSED   | 0.5474851    PASSED  | 0.7370361    PASSED  |
              sts_serial    4 | 0.7282286   PASSED   | 0.7282286    PASSED  | 0.2518826    PASSED  |
              sts_serial    4 | 0.9905724   PASSED   | 0.9905724    PASSED  | 0.6876253    PASSED  |
              sts_serial    5 | 0.8297711   PASSED   | 0.8297711    PASSED  | 0.2123014    PASSED  |
              sts_serial    5 | 0.9092172   PASSED   | 0.9092172    PASSED  | 0.3532615    PASSED  |
              sts_serial    6 | 0.4976615   PASSED   | 0.4976615    PASSED  | 0.9967160      **    |
              sts_serial    6 | 0.9853355   PASSED   | 0.9853355    PASSED  | 0.5537414    PASSED  |
              sts_serial    7 | 0.9675717   PASSED   | 0.9675717    PASSED  | 0.3804243    PASSED  |
              sts_serial    7 | 0.4446567   PASSED   | 0.4446567    PASSED  | 0.0923678    PASSED  |
              sts_serial    8 | 0.7254384   PASSED   | 0.7254384    PASSED  | 0.4544030    PASSED  |
              sts_serial    8 | 0.8984816   PASSED   | 0.8984816    PASSED  | 0.7501155    PASSED  |
              sts_serial    9 | 0.8255134   PASSED   | 0.8255134    PASSED  | 0.4260288    PASSED  |
              sts_serial    9 | 0.6609663   PASSED   | 0.6609663    PASSED  | 0.5622308    PASSED  |
              sts_serial   10 | 0.9984397     **     | 0.9984397      **    | 0.5789212    PASSED  |
              sts_serial   10 | 0.7987434   PASSED   | 0.7987434    PASSED  | 0.8599317    PASSED  |
              sts_serial   11 | 0.5552886   PASSED   | 0.5552886    PASSED  | 0.3546752    PASSED  |
              sts_serial   11 | 0.4417852   PASSED   | 0.4417852    PASSED  | 0.5042245    PASSED  |
              sts_serial   12 | 0.3843880   PASSED   | 0.3843880    PASSED  | 0.6723639    PASSED  |
              sts_serial   12 | 0.1514682   PASSED   | 0.1514682    PASSED  | 0.9428701    PASSED  |
              sts_serial   13 | 0.5396454   PASSED   | 0.5396454    PASSED  | 0.5793677    PASSED  |
              sts_serial   13 | 0.9497671   PASSED   | 0.9497671    PASSED  | 0.3370774    PASSED  |
              sts_serial   14 | 0.3616613   PASSED   | 0.3616613    PASSED  | 0.4372343    PASSED  |
              sts_serial   14 | 0.3996251   PASSED   | 0.3996251    PASSED  | 0.5185021    PASSED  |
              sts_serial   15 | 0.3847188   PASSED   | 0.3847188    PASSED  | 0.3188851    PASSED  |
              sts_serial   15 | 0.1012968   PASSED   | 0.1012968    PASSED  | 0.1631942    PASSED  |
              sts_serial   16 | 0.9974802     **     | 0.9974802      **    | 0.6645914    PASSED  |
              sts_serial   16 | 0.1157822   PASSED   | 0.1157822    PASSED  | 0.3465564    PASSED  |
             rgb_bitdist    1 | 0.4705599   PASSED   | 0.4705599    PASSED  | 0.8627740    PASSED  |
             rgb_bitdist    2 | 0.7578920   PASSED   | 0.7578920    PASSED  | 0.3296790    PASSED  |
             rgb_bitdist    3 | 0.9934502   PASSED   | 0.9934502    PASSED  | 0.5558012    PASSED  |
             rgb_bitdist    4 | 0.3674201   PASSED   | 0.3674201    PASSED  | 0.1607977    PASSED  |
             rgb_bitdist    5 | 0.7930273   PASSED   | 0.7930273    PASSED  | 0.9999802      **    |
             rgb_bitdist    6 | 0.8491477   PASSED   | 0.8491477    PASSED  | 0.3774760    PASSED  |
             rgb_bitdist    7 | 0.1537432   PASSED   | 0.1537432    PASSED  | 0.4715169    PASSED  |
             rgb_bitdist    8 | 0.9454030   PASSED   | 0.9454030    PASSED  | 0.9890644    PASSED  |
             rgb_bitdist    9 | 0.2017856   PASSED   | 0.2017856    PASSED  | 0.0571014    PASSED  |
             rgb_bitdist   10 | 0.9989305     **     | 0.9989305      **    | 0.4575834    PASSED  |
             rgb_bitdist   11 | 0.4441883   PASSED   | 0.4441883    PASSED  | 0.4960057    PASSED  |
             rgb_bitdist   12 | 0.7074388   PASSED   | 0.7074388    PASSED  | 0.6808850    PASSED  |
    rgb_minimum_distance    2 | 0.9604056   PASSED   | 0.9604056    PASSED  | 0.8859729    PASSED  |
    rgb_minimum_distance    3 | 0.5143592   PASSED   | 0.5143592    PASSED  | 0.3266204    PASSED  |
    rgb_minimum_distance    4 | 0.3779106   PASSED   | 0.3779106    PASSED  | 0.3537417    PASSED  |
    rgb_minimum_distance    5 | 0.4861264   PASSED   | 0.4861264    PASSED  | 0.9032057    PASSED  |
        rgb_permutations    2 | 0.9206310   PASSED   | 0.9206310    PASSED  | 0.8052940    PASSED  |
        rgb_permutations    3 | 0.9299743   PASSED   | 0.9299743    PASSED  | 0.2209750    PASSED  |
        rgb_permutations    4 | 0.8330345   PASSED   | 0.8330345    PASSED  | 0.5819945    PASSED  |
        rgb_permutations    5 | 0.2708879   PASSED   | 0.2708879    PASSED  | 0.9276941    PASSED  |
          rgb_lagged_sum    0 | 0.0794660   PASSED   | 0.0794660    PASSED  | 0.9918681    PASSED  |
          rgb_lagged_sum    1 | 0.5279555   PASSED   | 0.5279555    PASSED  | 0.1304600    PASSED  |
          rgb_lagged_sum    2 | 0.0433872   PASSED   | 0.0433872    PASSED  | 0.1149961    PASSED  |
          rgb_lagged_sum    3 | 0.0028004     **     | 0.0028004      **    | 0.2731577    PASSED  |
          rgb_lagged_sum    4 | 0.0000074     **     | 0.0000074      **    | 0.8978870    PASSED  |
          rgb_lagged_sum    5 | 0.1332411   PASSED   | 0.1332411    PASSED  | 0.2065880    PASSED  |
          rgb_lagged_sum    6 | 0.0412128   PASSED   | 0.0412128    PASSED  | 0.7611867    PASSED  |
          rgb_lagged_sum    7 | 0.0225446   PASSED   | 0.0225446    PASSED  | 0.4810145    PASSED  |
          rgb_lagged_sum    8 | 0.0087433   PASSED   | 0.0087433    PASSED  | 0.3120378    PASSED  |
          rgb_lagged_sum    9 | 0.0000000   ******   | 0.0000000    ******  | 0.1334315    PASSED  |
          rgb_lagged_sum   10 | 0.4147842   PASSED   | 0.4147842    PASSED  | 0.2334790    PASSED  |
          rgb_lagged_sum   11 | 0.0206564   PASSED   | 0.0206564    PASSED  | 0.6491578    PASSED  |
          rgb_lagged_sum   12 | 0.0755835   PASSED   | 0.0755835    PASSED  | 0.5332069    PASSED  |
          rgb_lagged_sum   13 | 0.3112028   PASSED   | 0.3112028    PASSED  | 0.4194447    PASSED  |
          rgb_lagged_sum   14 | 0.0000000   ******   | 0.0000000    ******  | 0.2584573    PASSED  |
          rgb_lagged_sum   15 | 0.0890059   PASSED   | 0.0890059    PASSED  | 0.0007064      **    |
          rgb_lagged_sum   16 | 0.2962076   PASSED   | 0.2962076    PASSED  | 0.1344984    PASSED  |
          rgb_lagged_sum   17 | 0.2696070   PASSED   | 0.2696070    PASSED  | 0.2242021    PASSED  |
          rgb_lagged_sum   18 | 0.0826388   PASSED   | 0.0826388    PASSED  | 0.0450341    PASSED  |
          rgb_lagged_sum   19 | 0.0000000   ******   | 0.0000000    ******  | 0.5508302    PASSED  |
          rgb_lagged_sum   20 | 0.0101437   PASSED   | 0.0101437    PASSED  | 0.4290150    PASSED  |
          rgb_lagged_sum   21 | 0.1417859   PASSED   | 0.1417859    PASSED  | 0.1624411    PASSED  |
          rgb_lagged_sum   22 | 0.0160264   PASSED   | 0.0160264    PASSED  | 0.5204838    PASSED  |
          rgb_lagged_sum   23 | 0.0535167   PASSED   | 0.0535167    PASSED  | 0.6571892    PASSED  |
          rgb_lagged_sum   24 | 0.0000000   ******   | 0.0000000    ******  | 0.8578906    PASSED  |
          rgb_lagged_sum   25 | 0.8453426   PASSED   | 0.8453426    PASSED  | 0.3568988    PASSED  |
          rgb_lagged_sum   26 | 0.2113484   PASSED   | 0.2113484    PASSED  | 0.9755715    PASSED  |
          rgb_lagged_sum   27 | 0.1903762   PASSED   | 0.1903762    PASSED  | 0.4356739    PASSED  |
          rgb_lagged_sum   28 | 0.0733066   PASSED   | 0.0733066    PASSED  | 0.8354990    PASSED  |
          rgb_lagged_sum   29 | 0.0000000   ******   | 0.0000000    ******  | 0.1716599    PASSED  |
          rgb_lagged_sum   30 | 0.0932124   PASSED   | 0.0932124    PASSED  | 0.0732090    PASSED  |
          rgb_lagged_sum   31 | 0.0000000   ******   | 0.0000000    ******  | 0.3497910    PASSED  |
          rgb_lagged_sum   32 | 0.0843455   PASSED   | 0.0843455    PASSED  | 0.5441949    PASSED  |
         rgb_kstest_test    0 | 0.4399862   PASSED   | 0.4399862    PASSED  | 0.9766581    PASSED  |
         dab_bytedistrib    0 | 0.0748312   PASSED   | 0.0748312    PASSED  | 0.7035800    PASSED  |
                 dab_dct  256 | 0.0919474   PASSED   | 0.0919474    PASSED  | 0.3985889    PASSED  |
            dab_filltree   32 | 0.1227533   PASSED   | 0.1227533    PASSED  | 0.7390925    PASSED  |
            dab_filltree   32 | 0.6819630   PASSED   | 0.6819630    PASSED  | 0.1773611    PASSED  |
           dab_filltree2    0 | 0.1774773   PASSED   | 0.1774773    PASSED  | 0.2088828    PASSED  |
           dab_filltree2    1 | 0.1718216   PASSED   | 0.1718216    PASSED  | 0.2257006    PASSED  |
            dab_monobit2   12 | 0.9999881     **     | 0.9999881      **    | 0.8084149    PASSED  |

     

     


     

    Download DIEHARD_test_of_MAPLE_MersenneTwister.mw

    A lot of supplementary details are given in the attached file.
    I let the readers discover by themselves if Maple's implementation of the Mersenne Twister PRNG is correct or not.
    Beyond this exercise, I hope this work will be useful to people who could be tempted to test their own generator.

     

     

    I am very pleased to announce that we have released a new version of the free Maple Companion app. For those you may have missed it, the first release of this app gave you a way to take a picture of math using your phone’s camera and upload it into Maple. Instructors have told me they’ve found this very useful in their classes, as they no longer have to deal with transcription errors as students enter problems into Maple.

    So that’s good. But version 2 is a lot better. The Maple Companion now solves math problems directly on your phone. It can handle problems from algebra, precalculus, calculus, linear algebra, differential equations, and more. No need to upload to Maple – students can solve the problem by hand, and then use the app to check their answer, try new operations on the same expression, and even create plots. And if they want to do even more, they can still upload the expression into Maple for more advanced operations and explorations.

    There’s also a built-in math editor, so you can enter problems without the camera, too. And if you use the camera, and it misinterprets part of your expression, you can fix it using the editor instead of having to retake the picture.  Good as the math recognition is, even in the face of some pretty poor handwriting, the ability to tweak the results has proven to be extremely useful.

    There’s lots more we’d like to do with the Maple Companion app over time, and we’d like hear your thoughts, as well. How else can it help students learn?  How else can it act as a complement to Maple? Let us know!

    Visit Maple Companion to learn more, find links to the app stores so you can download the app, and access the feedback form. And if you already have version 1, you can get the new release simply by updating the app on your phone.

     

    This update fixes the problems inadvertently introduced in Maple 2019.2, namely:

    • Maple failed to run the code in the maple.ini/.mapleinit initialization files when loading existing worksheets containing a restart() command
    • Installing some packages from the MapleCloud was unsuccessful

    For anyone who installed the 2019.2 update, installing 2019.2.1 will fix these problems.

    If you are at Maple 2019.1 or earlier, installing this update will bring you straight to Maple 2019.2.1.

    This update is available through Tools>Check for Updates in Maple, and is also available from our website on the Maple 2019.2.1 download page.

    If you are a MapleSim user, please note that these problems do not affect your use of MapleSim. If you use Maple on its own, and if you use Maple command initialization files and/or you need to install a package from the MapleCloud that does not work, please contact Maplesoft Technical Support for assistance.

    We sincerely apologize for the inconvenience and thank you for your patience as we worked through this issue.

    There was a discussion regarding Maple and two light sources a couple of days ago (unfortunately the content of that conversation has been removed by the original poster so if any of you thought you were getting early signs of alzheimers, rest assured the questions/posts were indeed there and you can all put your memories at ease for a little).

    I did indeed go back into Maple to see how far it was that we lost the two light source capability.  I can tell you that as far back as Maple 10.06 dual (or more) light sources was not working.  So I worked my way up - of course the OP of the deleted questions mentioned Maple Vr4 (I think), so I started there.  Maple 6 and Maple 7 is ok.  I was unable to check Maple 8 and 9 as my computers with that software are currently missing (probably buried in my shed somewhere - as long as the wife hasn't got rid of them.. anyways)

    So I'll share with you what once was.  From Maple 7, and an animated gif (which hopefully works) to illustrate two light sources (one brown light, and one mauve or purple light) and what it looked like from Maple 7.

    with(plots):
    setoptions3d(style=patchnogrid,projection=.9):
    a:=plot3d(-x^2-y^2-x*y,x=-1..1,y=-1..1,color=[.5,.9,.9],grid=[15,15]):
    b:=n->display(a,orientation=[n,60]):
    display([seq(b(10*n),n=0..35)],light=[90,-80,1,.2,.1],light=[90,40,1,.5,.8],insequence=true);

    The animation is not working so the uploaded gif file is here for you to check out. .. Couldn't upload the gif file so I've zipped it - that worked

    twolightsource.zip

    I'm only just hearing (haven't experienced) about some serious issues with the 2019.2 updates.  I would recommend waiting for Maplesoft to release an emergency 2019.3 fix update - Maplesoft can NOT leave the last update of 2019 in this state.

    After updating to Maple2019.2, the user initialization files maple.ini (Windows) and .mapleinit (Linux), which are placed into the user's home dierectories, are no longer read in upon startup. The behavior of Maple2019.1 was correct, the files were read.

    In order to estimate parameters of permanent magnet synchronous motor (PMSM) on-line and real-time, an adaptive on-line identification method for motor parameters is proposed. Resistance, inductance and PM flux of PMSM are achieved at the same time in the presented model. By means of Popov’s hyper-stability theory, the model of parameter identification is built in the rotor reference frame. And, PMSM d, q-axis voltage, current and their errors are used to obtain the adaptive laws of parameters. Popov’s hyper-stability theory guarantees stability of the system and convergence of the estimated parameters under certain conditions. The simulation and experimental results illustrate the validity and efficiency of the proposed method.

    We have just released updates to Maple and MapleSim.

    Maple 2019.2 includes corrections and improvements to a variety of areas in the product, including a new “Go to page ____” option in print preview (that am personally quite pleased about), sections are expanded by default when printing or exporting, a fix to a problem using non-executable math with text in document mode that sometimes made it impossible to advance to a new line using Enter, improvements to VectorCalculus, select, abs and other math functions, support for macOS Catalina, and more.  We recommend that all Maple 2019 users install these updates.

    This update is available through Tools>Check for Updates in Maple, and is also available from our website on the Maple 2019.2 download page, where you can also find more details.

    For MapleSim users, the MapleSim 2019.2 family of products includes enhancements in the areas of model development and toolchain connectivity, including substantial enhancements to the MapleSim CAD toolbox.   For more details and download instructions, visit the MapleSim 2019.2 download page.

    Program of Transfer matrix method for solving the vibrational behavior of a cracked beam. For more details and comprehension check the paper entitled ''A transfer matrix method for free vibration analysis and crack identification of stepped beams with multiple edge cracks and different boundary conditions ''
     

    restart; with(LinearAlgebra)NULLNULL

    W11 := A[1]*cos(nu*x)+A[2]*sin(nu*x)+A[3]*cosh(nu*x)+A[4]*sinh(nu*x);

    theta11 := -A[1]*nu*sin(nu*x)+A[2]*nu*cos(nu*x)+A[3]*nu*sinh(nu*x)+A[4]*nu*cosh(nu*x)  NULL

    M11 := EI*(-A[1]*nu^2*cos(nu*x)-A[2]*nu^2*sin(nu*x)+A[3]*nu^2*cosh(nu*x)+A[4]*nu^2*sinh(nu*x))
    S11 := EI*(A[1]*nu^3*sin(nu*x)-A[2]*nu^3*cos(nu*x)+A[3]*nu^3*sinh(nu*x)+A[4]*nu^3*cosh(nu*x))
     

    MD11 := subs(A[1] = 1, A[2] = 0, A[3] = 0, A[4] = 0, W11); MD12 := subs(A[1] = 0, A[2] = 1, A[3] = 0, A[4] = 0, W11); MD13 := subs(A[1] = 0, A[2] = 0, A[3] = 1, A[4] = 0, W11); MD14 := subs(A[1] = 0, A[2] = 0, A[3] = 0, A[4] = 1, W11)
    NULL

    MD21 := subs(A[1] = 1, A[2] = 0, A[3] = 0, A[4] = 0, theta11); MD22 := subs(A[1] = 0, A[2] = 1, A[3] = 0, A[4] = 0, theta11); MD23 := subs(A[1] = 0, A[2] = 0, A[3] = 1, A[4] = 0, theta11); MD24 := subs(A[1] = 0, A[2] = 0, A[3] = 0, A[4] = 1, theta11)
     

    MD31 := subs(A[1] = 1, A[2] = 0, A[3] = 0, A[4] = 0, M11); MD32 := subs(A[1] = 0, A[2] = 1, A[3] = 0, A[4] = 0, M11); MD33 := subs(A[1] = 0, A[2] = 0, A[3] = 1, A[4] = 0, M11); MD34 := subs(A[1] = 0, A[2] = 0, A[3] = 0, A[4] = 1, M11)
     

    MD41 := subs(A[1] = 1, A[2] = 0, A[3] = 0, A[4] = 0, S11); MD42 := subs(A[1] = 0, A[2] = 1, A[3] = 0, A[4] = 0, S11); MD43 := subs(A[1] = 0, A[2] = 0, A[3] = 1, A[4] = 0, S11); MD44 := subs(A[1] = 0, A[2] = 0, A[3] = 0, A[4] = 1, S11)
     

     

    TM := Matrix(4, 4, [[MD11, MD12, MD13, MD14], [MD21, MD22, MD23, MD24], [MD31, MD32, MD33, MD34], [MD41, MD42, MD43, MD44]])

    C := Matrix(4, 4, [[1, 0, 0, 0], [0, 1, c44, 0], [0, 0, 1, 0], [0, 0, 0, 1]])NULL

    NULL    TM2 := subs(x = 0, TM); TM3 := subs(x = L, TM)

    with(MTM)

    TM4 := inv(TM)NULLNULL 

    TM5 := inv(TM2)NULL

    Y11 := MatrixMatrixMultiply(TM3, TM4)

        Y22 := MatrixMatrixMultiply(C, TM)   

    Y33 := MatrixMatrixMultiply(Y11, Y22)

    Y44 := MatrixMatrixMultiply(Y33, TM5)NULL

    BB11 := Y44[3, 3]

    BB12 := Y44[3, 4]

    BB21 := Y44[4, 3]

    BB22 := Y44[4, 4] 

    BB := Matrix(2, 2, [[BB11, BB12], [BB21, BB22]]) 

    NULL

    R11 := det(BB) 

    NULL

    L := .18; L1 := 1; h := 0.5e-2; b := 0.2e-1; rho := 957.5; area = b.h; m := rho*h*b; EI := 0.2682e10*b*h^3mu := ((m.(omega^2))*L^4/EI)^(1/4); x := .5; c44 := .1; c11 := 0NULLNULL

    plot(R11, omega = 1 .. 100)

     

     

     

    ``

    NULL

    NULL


     

    Download transfer.mw

    In the applications I am working on, the information are often represented by hierarchical tables (that is tables where some entries can also be tables, and so on).
    To help people to understand how this information is organized, I have thought to representent this hierarchical table as a tree graph.
    Once this graph built, it becomes very simple to find where a "terminal leaf", that is en entry which is no longer a table, is located in the original table (by location I mean the sequence of indices for which the entry is this "terminal leaf".

    The code provided here is pretension free and I do not doubt a single second  that people here will be able to improve it.
    I published it for i thought other people could face the same kind of problems that I do.


     

    restart

    with(GraphTheory):
    interface(version);

    `Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

    (1)

    gh := proc(T)
      global s, counter, types:
      local  i:
      if type(T, table) then
        for i in [indices(T, nolist)] do
          if type(T[i], table) then
             s := s, op(map(u -> [i, u], [indices(T[i], nolist)] ));
          else
             counter := counter+1:
             types   := types, _Z_||counter = whattype(T[i]);
             s       := s, [i, _Z_||counter];
          end if:
          thisproc(T[i]):
        end do:
      else
        return s
      end if:
    end proc:

    t := table([a1=[alpha=1, beta=2], a2=table([a21=2, a22=table([a221=x, a222=table([a2221={1, 2, 3}, a2222=Matrix(2, 2), a2223=u3, a2224=u4])])]), a3=table([a31=u, a32=v])]);

    global s, counter, types:
    s       := NULL:
    counter := 0:
    types   := NULL:

    ghres := gh(t):
    types := [types]:

    t := table([a1 = [alpha = 1, beta = 2], a3 = table([a32 = v, a31 = u]), a2 = table([a22 = table([a222 = table([a2222 = (Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0})), a2223 = u3, a2221 = {1, 2, 3}, a2224 = u4]), a221 = x]), a21 = 2])])

    (2)


    These 3 lines determine the set of edges of the form ['t', v], that are not been captured by procedure h.
    They correspond to "first level" indices of table t (v in {a1, a2, a3} in the example above)

    L := convert(op~(1, [ghres]), set):     
    R := convert(op~(2, [ghres]), set):
    FirstLevelEdges := map(u -> ['t', u], L union R minus R):


    Complete the set of the edges, build the graph representation TG of table t and draw TG.

    edges := convert~({ghres, FirstLevelEdges[]}, set):
    TG := Graph(edges):

    HighlightVertex(TG, Vertices(TG), white):
    p := DrawGraph(TG, style=tree, root='t'):
     


    The first line is used to change the the "terminal leaves" of names  _Z_n by their type.

    eval(t);

    p       := subs(types, p):
    enlarge := plottools:-transform((x,y) -> [3*x, y]):

    plots:-display(enlarge(p), size=[1000, 400]);

    table([a1 = [alpha = 1, beta = 2], a3 = table([a32 = v, a31 = u]), a2 = table([a22 = table([a222 = table([a2222 = (Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0})), a2223 = u3, a2221 = {1, 2, 3}, a2224 = u4]), a221 = x]), a21 = 2])])

     

     


    This procedure is used to find the "indices path" to a terminal leaf.
    FindLeaf is then applied to all the terminal leaves.

    FindLeaf := proc(TG, leaf)
       local here:
       here := GraphTheory:-ShortestPath(TG, 't', leaf)[1..-2]:
       here := cat(convert(here[1], string), convert(here[2..-1], string)):
       here := StringTools:-SubstituteAll(here, ",", "]["):
       here := parse(here);
    end proc:

    # where is a2221

    printf("%a\n", FindLeaf(TG, a2221));

    t[a2][a22][a222]

     

     


     

    Download Table_Unfolding_2.mw

     


     

    with(LinearAlgebra); restart

    v1 := x^3*a[3]+x^2*a[2]+x*a[1]+a[0]:

    r := diff(v1, x):

    v[i] := subs(x = 0, v1):

    theta[i] := subs(x = 0, r):

    v[j] := subs(x = L, v1):

    theta[j] := subs(x = L, r):

    MD11 := subs(a[0] = 1, a[1] = 0, a[2] = 0, a[3] = 0, v[i]):

    MS11 := subs(a[0] = 0, a[1] = 0, a[2] = 0, a[3] = 0, theta[i]):

    MT11 := subs(a[0] = 1, a[1] = 0, a[2] = 0, a[3] = 0, v[j]):

    MR11 := subs(a[0] = 1, a[1] = 0, a[2] = 0, a[3] = 0, theta[j]):

    TM := Matrix(4, 4, [[MD11, MD12, MD13, MD14], [MS11, MS12, MS13, MS14], [MT11, MT12, MT13, MT14], [MR11, MR12, MR13, MR14]]):

    with(MTM):

    TM1 := inv(TM):

    b := Typesetting:-delayDotProduct(TM1, c):

    c := `<,>`(v[1], theta[1], v[2], theta[2]):

    a[0] := b(1):

    a[1] := b(2):

    a[2] := b(3):

    a[3] := b(4):

    vshape := x^3*a[3]+x^2*a[2]+x*a[1]+a[0]:

    XX1 := collect(vshape, v[1]):

    XX2 := collect(XX1, theta[1]):

    XX3 := collect(XX2, v[2]):

    XX4 := collect(XX3, theta[2]):

    v[1] := 1:

    NN1 := XX4:

    v[1] := 0:

    NN2 := XX4:

    v[1] := 0:

    NN3 := XX4:

    v[1] := 0:

    NN4 := XX4:

    Mat := `<,>`([NN1, NN2, NN3, NN4]):

    Mat1 := diff(Mat, x):

    Mat2 := diff(Mat1, x):

    segma := Mat2:

    KK := EI*expand(Typesetting:-delayDotProduct(segma, transpose(segma))):

    KG := int(KK, x = 0 .. L)

    KG := Matrix(4, 4, {(1, 1) = 12*EI/L^3, (1, 2) = 6*EI/L^2, (1, 3) = -12*EI/L^3, (1, 4) = 6*EI/L^2, (2, 1) = 6*EI/L^2, (2, 2) = 4*EI/L, (2, 3) = -6*EI/L^2, (2, 4) = 2*EI/L, (3, 1) = -12*EI/L^3, (3, 2) = -6*EI/L^2, (3, 3) = 12*EI/L^3, (3, 4) = -6*EI/L^2, (4, 1) = 6*EI/L^2, (4, 2) = 2*EI/L, (4, 3) = -6*EI/L^2, (4, 4) = 4*EI/L})

    (1)

    ``

    NULL


     

    Download stiffnesmatrice.mw

     Stiffness matrix implementation of a finite element beam based on Euler-Bernoulli theory stiffnesmatrice.mw


     

    restart: with(LinearAlgebra):

    # Motion equation (  Vibration of a cracked composite beam using general solution)  Edited by Adjal Yassine #

    ####################################################################

    Motion equation of flexural  vibration in normalized form 

    EI*W^(iv)-m*omega^2*W=0;

    EI*W^iv-m*omega^2*W = 0

    (1)

     

    The general solution form of bending vibration equation

    W1:=A[1]*cosh(mu*x)+A[2]*sinh(mu*x)+A[3]*cos(mu*x)+A[4]*sin(mu*x);

    A[1]*cosh(mu*x)+A[2]*sinh(mu*x)+A[3]*cos(mu*x)+A[4]*sin(mu*x)

    (2)

    where

    E:=2682e6;L:=0.18;h:=0.004;b:=0.02;rho:=2600;area=b*h;m:=rho*h*b;II:=(h*b^3)/12:

    0.2682e10

     

    .18

     

    0.4e-2

     

    0.2e-1

     

    2600

     

    area = 0.8e-4

     

    .20800

    (3)

    mu:=((m*omega^2*L^4/EI)^(1/4)):

     

     Expression of cross-sectional rotation , the bending moment shear  force and torsional moment  are given as follows respectively

    theta1 := (1/L)*(A[1]*mu*sinh(mu*x)+A[2]*mu*cosh(mu*x)-A[3]*mu*sin(mu*x)+A[4]*mu*cos(mu*x));

    (A[1]*mu*sinh(mu*x)+A[2]*mu*cosh(mu*x)-A[3]*mu*sin(mu*x)+A[4]*mu*cos(mu*x))/L

    (4)

    M1:= (EI/L^2)*(A[1]*mu^2*cosh(mu*x)+A[2]*mu^2*sinh(mu*x)-A[3]*mu^2*cos(mu*x)-A[4]*mu^2*sin(mu*x));

    EI*(A[1]*mu^2*cosh(mu*x)+A[2]*mu^2*sinh(mu*x)-A[3]*mu^2*cos(mu*x)-A[4]*mu^2*sin(mu*x))/L^2

    (5)

    S1:= (-EI/L^3)*(A[1]*mu^3*sinh(mu*x)+A[2]*mu^3*cosh(mu*x)+A[3]*mu^3*sin(mu*x)-A[4]*mu^3*cos(mu*x));

    -EI*(A[1]*mu^3*sinh(mu*x)+A[2]*mu^3*cosh(mu*x)+A[3]*mu^3*sin(mu*x)-A[4]*mu^3*cos(mu*x))/L^3

    (6)

     

    W2:=A[5]*cosh(mu*x)+A[6]*sinh(mu*x)+A[7]*cos(mu*x)+A[8]*sin(mu*x);

    A[5]*cosh(mu*x)+A[6]*sinh(mu*x)+A[7]*cos(mu*x)+A[8]*sin(mu*x)

    (7)

     

    theta2:= (1/L)*(A[5]*mu*sinh(mu*x)+A[6]*mu*cosh(mu*x)-A[7]*mu*sin(mu*x)+A[8]*mu*cos(mu*x));

    (A[5]*mu*sinh(mu*x)+A[6]*mu*cosh(mu*x)-A[7]*mu*sin(mu*x)+A[8]*mu*cos(mu*x))/L

    (8)

    M2:= (EI/L^2)*(A[5]*mu^2*cosh(mu*x)+A[6]*mu^2*sinh(mu*x)-A[7]*mu^2*cos(mu*x)-A[8]*mu^2*sin(mu*x));

    EI*(A[5]*mu^2*cosh(mu*x)+A[6]*mu^2*sinh(mu*x)-A[7]*mu^2*cos(mu*x)-A[8]*mu^2*sin(mu*x))/L^2

    (9)

    S2:= -(EI/L^3)*(A[5]*mu^3*sinh(mu*x)+A[6]*mu^3*cosh(mu*x)+A[7]*mu^3*sin(mu*x)-A[8]*mu^3*cos(mu*x));

    -EI*(A[5]*mu^3*sinh(mu*x)+A[6]*mu^3*cosh(mu*x)+A[7]*mu^3*sin(mu*x)-A[8]*mu^3*cos(mu*x))/L^3

    (10)

     

    The boundary conditions at fixed end W1(0)=Theta(0)=0

    X1:=eval(subs(x=0,W1));

    A[1]+A[3]

    (11)

    X2:=eval(subs(x=0,theta1));

    (mu*A[2]+mu*A[4])/L

    (12)

    X2:=collect(X2,mu)*(L/mu);

    A[2]+A[4]

    (13)

     

    The boundary condtions at free end M2(1)=S2(1)=0

    X3:=eval(subs(x=1,M2));

    EI*(A[5]*mu^2*cosh(mu)+A[6]*mu^2*sinh(mu)-A[7]*mu^2*cos(mu)-A[8]*mu^2*sin(mu))/L^2

    (14)

    X3:=collect(X3,mu)*(L^2/mu^2/EI);

    cosh(mu)*A[5]+sinh(mu)*A[6]-cos(mu)*A[7]-sin(mu)*A[8]

    (15)

    X4:=eval(subs(x=1,S2));

    -EI*(A[5]*mu^3*sinh(mu)+A[6]*mu^3*cosh(mu)+A[7]*mu^3*sin(mu)-A[8]*mu^3*cos(mu))/L^3

    (16)

    X4:=collect(X4,mu);

    -EI*(cosh(mu)*A[6]+sinh(mu)*A[5]-cos(mu)*A[8]+sin(mu)*A[7])*mu^3/L^3

    (17)

    X4:=collect(X4,EI)*(L^3/mu^3/EI);

    -cosh(mu)*A[6]-sinh(mu)*A[5]+cos(mu)*A[8]-sin(mu)*A[7]

    (18)

     

    The additional boundary conditions at crack location

    X5:=combine(M1-M2);

    (EI*cosh(mu*x)*mu^2*A[1]-EI*cosh(mu*x)*mu^2*A[5]+EI*sinh(mu*x)*mu^2*A[2]-EI*sinh(mu*x)*mu^2*A[6]-EI*cos(mu*x)*mu^2*A[3]+EI*cos(mu*x)*mu^2*A[7]-EI*sin(mu*x)*mu^2*A[4]+EI*sin(mu*x)*mu^2*A[8])/L^2

    (19)

    X5:=collect(X5,mu);

    (EI*cosh(mu*x)*A[1]-EI*cosh(mu*x)*A[5]+EI*sinh(mu*x)*A[2]-EI*sinh(mu*x)*A[6]-cos(mu*x)*EI*A[3]+A[7]*cos(mu*x)*EI-A[4]*sin(mu*x)*EI+A[8]*sin(mu*x)*EI)*mu^2/L^2

    (20)

    X5:=collect(X5,EI)*(L^2/mu^2/EI);

    A[1]*cosh(mu*x)-A[5]*cosh(mu*x)+A[2]*sinh(mu*x)-A[6]*sinh(mu*x)-A[3]*cos(mu*x)+A[7]*cos(mu*x)-A[4]*sin(mu*x)+A[8]*sin(mu*x)

    (21)

    X6:=combine(S1-S2);

    (-EI*cosh(mu*x)*mu^3*A[2]+EI*cosh(mu*x)*mu^3*A[6]-EI*sinh(mu*x)*mu^3*A[1]+EI*sinh(mu*x)*mu^3*A[5]+EI*cos(mu*x)*mu^3*A[4]-EI*cos(mu*x)*mu^3*A[8]-EI*sin(mu*x)*mu^3*A[3]+EI*sin(mu*x)*mu^3*A[7])/L^3

    (22)

    X6:=collect(X6,mu);

    (-EI*cosh(mu*x)*A[2]+EI*cosh(mu*x)*A[6]-EI*sinh(mu*x)*A[1]+EI*A[5]*sinh(mu*x)+cos(mu*x)*A[4]*EI-cos(mu*x)*A[8]*EI-sin(mu*x)*EI*A[3]+sin(mu*x)*A[7]*EI)*mu^3/L^3

    (23)

    X6:=collect(X6,EI)*(L^3/mu^3/EI);

    -cosh(mu*x)*A[2]+cosh(mu*x)*A[6]-sinh(mu*x)*A[1]+sinh(mu*x)*A[5]+cos(mu*x)*A[4]-cos(mu*x)*A[8]-sin(mu*x)*A[3]+sin(mu*x)*A[7]

    (24)

     

    X7:=combine(W2-(W1+c8*S1));

    (EI*cosh(mu*x)*c8*mu^3*A[2]+EI*sinh(mu*x)*c8*mu^3*A[1]-EI*cos(mu*x)*c8*mu^3*A[4]+EI*sin(mu*x)*c8*mu^3*A[3]-cosh(mu*x)*A[1]*L^3+cosh(mu*x)*A[5]*L^3-sinh(mu*x)*A[2]*L^3+sinh(mu*x)*A[6]*L^3-cos(mu*x)*A[3]*L^3+cos(mu*x)*A[7]*L^3-sin(mu*x)*A[4]*L^3+sin(mu*x)*A[8]*L^3)/L^3

    (25)

    X8:=combine (theta2-(theta1+c44*M1));

    (-EI*cosh(mu*x)*c44*mu^2*A[1]-EI*sinh(mu*x)*c44*mu^2*A[2]+EI*cos(mu*x)*c44*mu^2*A[3]+EI*sin(mu*x)*c44*mu^2*A[4]-L*cosh(mu*x)*mu*A[2]+L*cosh(mu*x)*mu*A[6]-L*sinh(mu*x)*mu*A[1]+L*sinh(mu*x)*mu*A[5]-L*cos(mu*x)*mu*A[4]+L*cos(mu*x)*mu*A[8]+L*sin(mu*x)*mu*A[3]-L*sin(mu*x)*mu*A[7])/L^2

    (26)

     

    The characteristic matrix function of frequency

    FD8:=subs(A[1]=1,A[3]=0,X1);FD12:=subs(A[1]=0,A[3]=0,X1);FD13:=subs(A[1]=0,A[3]=1,X1);FD14:=subs(A[1]=0,A[3]=0,X1);FD15:=subs(A[1]=0,A[3]=0,X1);FD16:=subs(A[1]=0,A[3]=0,X1);FD17:=subs(A[1]=0,A[3]=0,X1);FD18:=subs(A[1]=0,A[3]=0,X1);

    1

     

    0

     

    1

     

    0

     

    0

     

    0

     

    0

     

    0

    (27)

    FD21:=subs(A[2]=0,A[4]=0,X2);FD22:=subs(A[2]=1,A[4]=0,X2);FD23:=subs(A[2]=0,A[4]=0,X2);FD24:=subs(A[2]=0,A[4]=1,X2);FD25:=subs(A[2]=0,A[4]=0,X2);FD26:=subs(A[2]=0,A[4]=0,X2);FD27:=subs(A[2]=0,A[4]=0,X2);FD28:=subs(A[2]=0,A[4]=0,X2);

    0

     

    1

     

    0

     

    1

     

    0

     

    0

     

    0

     

    0

    (28)

     

    FD31:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD32:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD33:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD34:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD35:=subs(A[5]=1,A[6]=0,A[7]=0,A[8]=0,X3);;FD36:=subs(A[5]=0,A[6]=1,A[7]=0,A[8]=0,X3);FD37:=subs(A[5]=0,A[6]=0,A[7]=1,A[8]=0,X3);FD38:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=1,X3);

    0

     

    0

     

    0

     

    0

     

    cosh(mu)

     

    sinh(mu)

     

    -cos(mu)

     

    -sin(mu)

    (29)

    FD41:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD42:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD43:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD44:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD45:=subs(A[5]=1,A[6]=0,A[7]=0,A[8]=0,X4);FD46:=subs(A[5]=0,A[6]=1,A[7]=0,A[8]=0,X4);FD47:=subs(A[5]=0,A[6]=0,A[7]=1,A[8]=0,X4);FD48:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=1,X4);

    0

     

    0

     

    0

     

    0

     

    -sinh(mu)

     

    -cosh(mu)

     

    -sin(mu)

     

    cos(mu)

    (30)

     

    FD51:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD52:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD53:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD54:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD55:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X5);FD56:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X5);FD57:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X5);FD58:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X5);

    cosh(mu*x)

     

    sinh(mu*x)

     

    -cos(mu*x)

     

    -sin(mu*x)

     

    -cosh(mu*x)

     

    -sinh(mu*x)

     

    cos(mu*x)

     

    sin(mu*x)

    (31)

    FD61:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD62:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD63:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD64:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD65:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X6);FD66:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X6);FD67:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X6);FD68:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X6);

    -sinh(mu*x)

     

    -cosh(mu*x)

     

    -sin(mu*x)

     

    cos(mu*x)

     

    sinh(mu*x)

     

    cosh(mu*x)

     

    sin(mu*x)

     

    -cos(mu*x)

    (32)

     

    FD71:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD72:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD73:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD74:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD75:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X7);FD76:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X7);FD77:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X7);FD78:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X7);

    (EI*sinh(mu*x)*c8*mu^3-cosh(mu*x)*L^3)/L^3

     

    (EI*cosh(mu*x)*c8*mu^3-sinh(mu*x)*L^3)/L^3

     

    (EI*sin(mu*x)*c8*mu^3-L^3*cos(mu*x))/L^3

     

    (-EI*cos(mu*x)*c8*mu^3-sin(mu*x)*L^3)/L^3

     

    cosh(mu*x)

     

    sinh(mu*x)

     

    cos(mu*x)

     

    sin(mu*x)

    (33)

    FD81:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD82:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD83:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD84:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD85:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X8);FD86:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X8);FD87:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X8);FD88:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X8);

    (-EI*cosh(mu*x)*c44*mu^2-L*sinh(mu*x)*mu)/L^2

     

    (-EI*sinh(mu*x)*c44*mu^2-L*cosh(mu*x)*mu)/L^2

     

    (EI*cos(mu*x)*c44*mu^2+L*sin(mu*x)*mu)/L^2

     

    (EI*sin(mu*x)*c44*mu^2-L*cos(mu*x)*mu)/L^2

     

    sinh(mu*x)*mu/L

     

    cosh(mu*x)*mu/L

     

    -sin(mu*x)*mu/L

     

    cos(mu*x)*mu/L

    (34)

     

    MM:=matrix(8,8,[[FD11,FD12,FD13,FD14,FD15,FD16,FD17,FD18],[FD21,FD22,FD23,FD24,FD25,FD26,FD27,FD28],[FD31,FD32,FD33,FD34,FD35,FD36,FD37,FD38],[FD41,FD42,FD43,FD44,FD45,FD46,FD47,FD48],[FD51,FD52,FD53,FD54,FD55,FD56,FD57,FD58],[FD61,FD62,FD63,FD64,FD65,FD66,FD67,FD68],[FD71,FD72,FD73,FD74,FD75,FD76,FD77,FD78],[FD81,FD82,FD83,FD84,FD85,FD86,FD87,FD88]]);

    MM := Matrix(8, 8, {(1, 1) = FD11, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 1, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = cosh(mu), (3, 6) = sinh(mu), (3, 7) = -cos(mu), (3, 8) = -sin(mu), (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = -sinh(mu), (4, 6) = -cosh(mu), (4, 7) = -sin(mu), (4, 8) = cos(mu), (5, 1) = cosh(mu*x), (5, 2) = sinh(mu*x), (5, 3) = -cos(mu*x), (5, 4) = -sin(mu*x), (5, 5) = -cosh(mu*x), (5, 6) = -sinh(mu*x), (5, 7) = cos(mu*x), (5, 8) = sin(mu*x), (6, 1) = -sinh(mu*x), (6, 2) = -cosh(mu*x), (6, 3) = -sin(mu*x), (6, 4) = cos(mu*x), (6, 5) = sinh(mu*x), (6, 6) = cosh(mu*x), (6, 7) = sin(mu*x), (6, 8) = -cos(mu*x), (7, 1) = (EI*sinh(mu*x)*c8*mu^3-cosh(mu*x)*L^3)/L^3, (7, 2) = (EI*cosh(mu*x)*c8*mu^3-sinh(mu*x)*L^3)/L^3, (7, 3) = (EI*sin(mu*x)*c8*mu^3-L^3*cos(mu*x))/L^3, (7, 4) = (-EI*cos(mu*x)*c8*mu^3-sin(mu*x)*L^3)/L^3, (7, 5) = cosh(mu*x), (7, 6) = sinh(mu*x), (7, 7) = cos(mu*x), (7, 8) = sin(mu*x), (8, 1) = (-EI*cosh(mu*x)*c44*mu^2-L*sinh(mu*x)*mu)/L^2, (8, 2) = (-EI*sinh(mu*x)*c44*mu^2-L*cosh(mu*x)*mu)/L^2, (8, 3) = (EI*cos(mu*x)*c44*mu^2+L*sin(mu*x)*mu)/L^2, (8, 4) = (EI*sin(mu*x)*c44*mu^2-L*cos(mu*x)*mu)/L^2, (8, 5) = sinh(mu*x)*mu/L, (8, 6) = cosh(mu*x)*mu/L, (8, 7) = -sin(mu*x)*mu/L, (8, 8) = cos(mu*x)*mu/L})

    (35)

    Program end

     

    NULL

     

    ``


     

    Download Vibration_of_a_cracked_composite_beam.mw
     

    restart: with(LinearAlgebra):

    # Motion equation (  Vibration of a cracked composite beam using general solution)  Edited by Adjal Yassine #

    ####################################################################

    Motion equation of flexural  vibration in normalized form 

    EI*W^(iv)-m*omega^2*W=0;

    EI*W^iv-m*omega^2*W = 0

    (1)

     

    The general solution form of bending vibration equation

    W1:=A[1]*cosh(mu*x)+A[2]*sinh(mu*x)+A[3]*cos(mu*x)+A[4]*sin(mu*x);

    A[1]*cosh(mu*x)+A[2]*sinh(mu*x)+A[3]*cos(mu*x)+A[4]*sin(mu*x)

    (2)

    where

    E:=2682e6;L:=0.18;h:=0.004;b:=0.02;rho:=2600;area=b*h;m:=rho*h*b;II:=(h*b^3)/12:

    0.2682e10

     

    .18

     

    0.4e-2

     

    0.2e-1

     

    2600

     

    area = 0.8e-4

     

    .20800

    (3)

    mu:=((m*omega^2*L^4/EI)^(1/4)):

     

     Expression of cross-sectional rotation , the bending moment shear  force and torsional moment  are given as follows respectively

    theta1 := (1/L)*(A[1]*mu*sinh(mu*x)+A[2]*mu*cosh(mu*x)-A[3]*mu*sin(mu*x)+A[4]*mu*cos(mu*x));

    (A[1]*mu*sinh(mu*x)+A[2]*mu*cosh(mu*x)-A[3]*mu*sin(mu*x)+A[4]*mu*cos(mu*x))/L

    (4)

    M1:= (EI/L^2)*(A[1]*mu^2*cosh(mu*x)+A[2]*mu^2*sinh(mu*x)-A[3]*mu^2*cos(mu*x)-A[4]*mu^2*sin(mu*x));

    EI*(A[1]*mu^2*cosh(mu*x)+A[2]*mu^2*sinh(mu*x)-A[3]*mu^2*cos(mu*x)-A[4]*mu^2*sin(mu*x))/L^2

    (5)

    S1:= (-EI/L^3)*(A[1]*mu^3*sinh(mu*x)+A[2]*mu^3*cosh(mu*x)+A[3]*mu^3*sin(mu*x)-A[4]*mu^3*cos(mu*x));

    -EI*(A[1]*mu^3*sinh(mu*x)+A[2]*mu^3*cosh(mu*x)+A[3]*mu^3*sin(mu*x)-A[4]*mu^3*cos(mu*x))/L^3

    (6)

     

    W2:=A[5]*cosh(mu*x)+A[6]*sinh(mu*x)+A[7]*cos(mu*x)+A[8]*sin(mu*x);

    A[5]*cosh(mu*x)+A[6]*sinh(mu*x)+A[7]*cos(mu*x)+A[8]*sin(mu*x)

    (7)

     

    theta2:= (1/L)*(A[5]*mu*sinh(mu*x)+A[6]*mu*cosh(mu*x)-A[7]*mu*sin(mu*x)+A[8]*mu*cos(mu*x));

    (A[5]*mu*sinh(mu*x)+A[6]*mu*cosh(mu*x)-A[7]*mu*sin(mu*x)+A[8]*mu*cos(mu*x))/L

    (8)

    M2:= (EI/L^2)*(A[5]*mu^2*cosh(mu*x)+A[6]*mu^2*sinh(mu*x)-A[7]*mu^2*cos(mu*x)-A[8]*mu^2*sin(mu*x));

    EI*(A[5]*mu^2*cosh(mu*x)+A[6]*mu^2*sinh(mu*x)-A[7]*mu^2*cos(mu*x)-A[8]*mu^2*sin(mu*x))/L^2

    (9)

    S2:= -(EI/L^3)*(A[5]*mu^3*sinh(mu*x)+A[6]*mu^3*cosh(mu*x)+A[7]*mu^3*sin(mu*x)-A[8]*mu^3*cos(mu*x));

    -EI*(A[5]*mu^3*sinh(mu*x)+A[6]*mu^3*cosh(mu*x)+A[7]*mu^3*sin(mu*x)-A[8]*mu^3*cos(mu*x))/L^3

    (10)

     

    The boundary conditions at fixed end W1(0)=Theta(0)=0

    X1:=eval(subs(x=0,W1));

    A[1]+A[3]

    (11)

    X2:=eval(subs(x=0,theta1));

    (mu*A[2]+mu*A[4])/L

    (12)

    X2:=collect(X2,mu)*(L/mu);

    A[2]+A[4]

    (13)

     

    The boundary condtions at free end M2(1)=S2(1)=0

    X3:=eval(subs(x=1,M2));

    EI*(A[5]*mu^2*cosh(mu)+A[6]*mu^2*sinh(mu)-A[7]*mu^2*cos(mu)-A[8]*mu^2*sin(mu))/L^2

    (14)

    X3:=collect(X3,mu)*(L^2/mu^2/EI);

    cosh(mu)*A[5]+sinh(mu)*A[6]-cos(mu)*A[7]-sin(mu)*A[8]

    (15)

    X4:=eval(subs(x=1,S2));

    -EI*(A[5]*mu^3*sinh(mu)+A[6]*mu^3*cosh(mu)+A[7]*mu^3*sin(mu)-A[8]*mu^3*cos(mu))/L^3

    (16)

    X4:=collect(X4,mu);

    -EI*(cosh(mu)*A[6]+sinh(mu)*A[5]-cos(mu)*A[8]+sin(mu)*A[7])*mu^3/L^3

    (17)

    X4:=collect(X4,EI)*(L^3/mu^3/EI);

    -cosh(mu)*A[6]-sinh(mu)*A[5]+cos(mu)*A[8]-sin(mu)*A[7]

    (18)

     

    The additional boundary conditions at crack location

    X5:=combine(M1-M2);

    (EI*cosh(mu*x)*mu^2*A[1]-EI*cosh(mu*x)*mu^2*A[5]+EI*sinh(mu*x)*mu^2*A[2]-EI*sinh(mu*x)*mu^2*A[6]-EI*cos(mu*x)*mu^2*A[3]+EI*cos(mu*x)*mu^2*A[7]-EI*sin(mu*x)*mu^2*A[4]+EI*sin(mu*x)*mu^2*A[8])/L^2

    (19)

    X5:=collect(X5,mu);

    (EI*cosh(mu*x)*A[1]-EI*cosh(mu*x)*A[5]+EI*sinh(mu*x)*A[2]-EI*sinh(mu*x)*A[6]-cos(mu*x)*EI*A[3]+A[7]*cos(mu*x)*EI-A[4]*sin(mu*x)*EI+A[8]*sin(mu*x)*EI)*mu^2/L^2

    (20)

    X5:=collect(X5,EI)*(L^2/mu^2/EI);

    A[1]*cosh(mu*x)-A[5]*cosh(mu*x)+A[2]*sinh(mu*x)-A[6]*sinh(mu*x)-A[3]*cos(mu*x)+A[7]*cos(mu*x)-A[4]*sin(mu*x)+A[8]*sin(mu*x)

    (21)

    X6:=combine(S1-S2);

    (-EI*cosh(mu*x)*mu^3*A[2]+EI*cosh(mu*x)*mu^3*A[6]-EI*sinh(mu*x)*mu^3*A[1]+EI*sinh(mu*x)*mu^3*A[5]+EI*cos(mu*x)*mu^3*A[4]-EI*cos(mu*x)*mu^3*A[8]-EI*sin(mu*x)*mu^3*A[3]+EI*sin(mu*x)*mu^3*A[7])/L^3

    (22)

    X6:=collect(X6,mu);

    (-EI*cosh(mu*x)*A[2]+EI*cosh(mu*x)*A[6]-EI*sinh(mu*x)*A[1]+EI*A[5]*sinh(mu*x)+cos(mu*x)*A[4]*EI-cos(mu*x)*A[8]*EI-sin(mu*x)*EI*A[3]+sin(mu*x)*A[7]*EI)*mu^3/L^3

    (23)

    X6:=collect(X6,EI)*(L^3/mu^3/EI);

    -cosh(mu*x)*A[2]+cosh(mu*x)*A[6]-sinh(mu*x)*A[1]+sinh(mu*x)*A[5]+cos(mu*x)*A[4]-cos(mu*x)*A[8]-sin(mu*x)*A[3]+sin(mu*x)*A[7]

    (24)

     

    X7:=combine(W2-(W1+c8*S1));

    (EI*cosh(mu*x)*c8*mu^3*A[2]+EI*sinh(mu*x)*c8*mu^3*A[1]-EI*cos(mu*x)*c8*mu^3*A[4]+EI*sin(mu*x)*c8*mu^3*A[3]-cosh(mu*x)*A[1]*L^3+cosh(mu*x)*A[5]*L^3-sinh(mu*x)*A[2]*L^3+sinh(mu*x)*A[6]*L^3-cos(mu*x)*A[3]*L^3+cos(mu*x)*A[7]*L^3-sin(mu*x)*A[4]*L^3+sin(mu*x)*A[8]*L^3)/L^3

    (25)

    X8:=combine (theta2-(theta1+c44*M1));

    (-EI*cosh(mu*x)*c44*mu^2*A[1]-EI*sinh(mu*x)*c44*mu^2*A[2]+EI*cos(mu*x)*c44*mu^2*A[3]+EI*sin(mu*x)*c44*mu^2*A[4]-L*cosh(mu*x)*mu*A[2]+L*cosh(mu*x)*mu*A[6]-L*sinh(mu*x)*mu*A[1]+L*sinh(mu*x)*mu*A[5]-L*cos(mu*x)*mu*A[4]+L*cos(mu*x)*mu*A[8]+L*sin(mu*x)*mu*A[3]-L*sin(mu*x)*mu*A[7])/L^2

    (26)

     

    The characteristic matrix function of frequency

    FD8:=subs(A[1]=1,A[3]=0,X1);FD12:=subs(A[1]=0,A[3]=0,X1);FD13:=subs(A[1]=0,A[3]=1,X1);FD14:=subs(A[1]=0,A[3]=0,X1);FD15:=subs(A[1]=0,A[3]=0,X1);FD16:=subs(A[1]=0,A[3]=0,X1);FD17:=subs(A[1]=0,A[3]=0,X1);FD18:=subs(A[1]=0,A[3]=0,X1);

    1

     

    0

     

    1

     

    0

     

    0

     

    0

     

    0

     

    0

    (27)

    FD21:=subs(A[2]=0,A[4]=0,X2);FD22:=subs(A[2]=1,A[4]=0,X2);FD23:=subs(A[2]=0,A[4]=0,X2);FD24:=subs(A[2]=0,A[4]=1,X2);FD25:=subs(A[2]=0,A[4]=0,X2);FD26:=subs(A[2]=0,A[4]=0,X2);FD27:=subs(A[2]=0,A[4]=0,X2);FD28:=subs(A[2]=0,A[4]=0,X2);

    0

     

    1

     

    0

     

    1

     

    0

     

    0

     

    0

     

    0

    (28)

     

    FD31:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD32:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD33:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD34:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD35:=subs(A[5]=1,A[6]=0,A[7]=0,A[8]=0,X3);;FD36:=subs(A[5]=0,A[6]=1,A[7]=0,A[8]=0,X3);FD37:=subs(A[5]=0,A[6]=0,A[7]=1,A[8]=0,X3);FD38:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=1,X3);

    0

     

    0

     

    0

     

    0

     

    cosh(mu)

     

    sinh(mu)

     

    -cos(mu)

     

    -sin(mu)

    (29)

    FD41:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD42:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD43:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD44:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD45:=subs(A[5]=1,A[6]=0,A[7]=0,A[8]=0,X4);FD46:=subs(A[5]=0,A[6]=1,A[7]=0,A[8]=0,X4);FD47:=subs(A[5]=0,A[6]=0,A[7]=1,A[8]=0,X4);FD48:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=1,X4);

    0

     

    0

     

    0

     

    0

     

    -sinh(mu)

     

    -cosh(mu)

     

    -sin(mu)

     

    cos(mu)

    (30)

     

    FD51:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD52:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD53:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD54:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD55:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X5);FD56:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X5);FD57:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X5);FD58:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X5);

    cosh(mu*x)

     

    sinh(mu*x)

     

    -cos(mu*x)

     

    -sin(mu*x)

     

    -cosh(mu*x)

     

    -sinh(mu*x)

     

    cos(mu*x)

     

    sin(mu*x)

    (31)

    FD61:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD62:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD63:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD64:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD65:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X6);FD66:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X6);FD67:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X6);FD68:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X6);

    -sinh(mu*x)

     

    -cosh(mu*x)

     

    -sin(mu*x)

     

    cos(mu*x)

     

    sinh(mu*x)

     

    cosh(mu*x)

     

    sin(mu*x)

     

    -cos(mu*x)

    (32)

     

    FD71:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD72:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD73:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD74:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD75:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X7);FD76:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X7);FD77:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X7);FD78:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X7);

    (EI*sinh(mu*x)*c8*mu^3-cosh(mu*x)*L^3)/L^3

     

    (EI*cosh(mu*x)*c8*mu^3-sinh(mu*x)*L^3)/L^3

     

    (EI*sin(mu*x)*c8*mu^3-L^3*cos(mu*x))/L^3

     

    (-EI*cos(mu*x)*c8*mu^3-sin(mu*x)*L^3)/L^3

     

    cosh(mu*x)

     

    sinh(mu*x)

     

    cos(mu*x)

     

    sin(mu*x)

    (33)

    FD81:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD82:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD83:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD84:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD85:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X8);FD86:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X8);FD87:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X8);FD88:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X8);

    (-EI*cosh(mu*x)*c44*mu^2-L*sinh(mu*x)*mu)/L^2

     

    (-EI*sinh(mu*x)*c44*mu^2-L*cosh(mu*x)*mu)/L^2

     

    (EI*cos(mu*x)*c44*mu^2+L*sin(mu*x)*mu)/L^2

     

    (EI*sin(mu*x)*c44*mu^2-L*cos(mu*x)*mu)/L^2

     

    sinh(mu*x)*mu/L

     

    cosh(mu*x)*mu/L

     

    -sin(mu*x)*mu/L

     

    cos(mu*x)*mu/L

    (34)

     

    MM:=matrix(8,8,[[FD11,FD12,FD13,FD14,FD15,FD16,FD17,FD18],[FD21,FD22,FD23,FD24,FD25,FD26,FD27,FD28],[FD31,FD32,FD33,FD34,FD35,FD36,FD37,FD38],[FD41,FD42,FD43,FD44,FD45,FD46,FD47,FD48],[FD51,FD52,FD53,FD54,FD55,FD56,FD57,FD58],[FD61,FD62,FD63,FD64,FD65,FD66,FD67,FD68],[FD71,FD72,FD73,FD74,FD75,FD76,FD77,FD78],[FD81,FD82,FD83,FD84,FD85,FD86,FD87,FD88]]);

    MM := Matrix(8, 8, {(1, 1) = FD11, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 1, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = cosh(mu), (3, 6) = sinh(mu), (3, 7) = -cos(mu), (3, 8) = -sin(mu), (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = -sinh(mu), (4, 6) = -cosh(mu), (4, 7) = -sin(mu), (4, 8) = cos(mu), (5, 1) = cosh(mu*x), (5, 2) = sinh(mu*x), (5, 3) = -cos(mu*x), (5, 4) = -sin(mu*x), (5, 5) = -cosh(mu*x), (5, 6) = -sinh(mu*x), (5, 7) = cos(mu*x), (5, 8) = sin(mu*x), (6, 1) = -sinh(mu*x), (6, 2) = -cosh(mu*x), (6, 3) = -sin(mu*x), (6, 4) = cos(mu*x), (6, 5) = sinh(mu*x), (6, 6) = cosh(mu*x), (6, 7) = sin(mu*x), (6, 8) = -cos(mu*x), (7, 1) = (EI*sinh(mu*x)*c8*mu^3-cosh(mu*x)*L^3)/L^3, (7, 2) = (EI*cosh(mu*x)*c8*mu^3-sinh(mu*x)*L^3)/L^3, (7, 3) = (EI*sin(mu*x)*c8*mu^3-L^3*cos(mu*x))/L^3, (7, 4) = (-EI*cos(mu*x)*c8*mu^3-sin(mu*x)*L^3)/L^3, (7, 5) = cosh(mu*x), (7, 6) = sinh(mu*x), (7, 7) = cos(mu*x), (7, 8) = sin(mu*x), (8, 1) = (-EI*cosh(mu*x)*c44*mu^2-L*sinh(mu*x)*mu)/L^2, (8, 2) = (-EI*sinh(mu*x)*c44*mu^2-L*cosh(mu*x)*mu)/L^2, (8, 3) = (EI*cos(mu*x)*c44*mu^2+L*sin(mu*x)*mu)/L^2, (8, 4) = (EI*sin(mu*x)*c44*mu^2-L*cos(mu*x)*mu)/L^2, (8, 5) = sinh(mu*x)*mu/L, (8, 6) = cosh(mu*x)*mu/L, (8, 7) = -sin(mu*x)*mu/L, (8, 8) = cos(mu*x)*mu/L})

    (35)

    Program end

     

    NULL

     

    ``


     

    Download Vibration_of_a_cracked_composite_beam.mwVibration_of_a_cracked_composite_beam.mwVibration_of_a_cracked_composite_beam.mw

     

    Splitting PDE parameterized symmetries

    and Parameter-continuous symmetry transformations

    The determination of symmetries for partial differential equation systems (PDE) is relevant in several contexts, the most obvious of which is of course the determination of the PDE solutions. For instance, generally speaking, the knowledge of a N-dimensional Lie symmetry group can be used to reduce the number of independent variables of PDE by N. So if PDE depends only on N independent variables, that amounts to completely solving it. If only N-1 symmetries are known or can be successfully used then PDE becomes and ODE; etc., all advantageous situations. In Maple, a complete set of symmetry commands, to perform each step of the symmetry approach or several of them in one go, is part of the PDEtools  package.

     

    Besides the dependent and independent variables, PDE frequently depends on some constant parameters, and besides the PDE symmetries for arbitrary values of those parameters, for some particular values of them, PDE transforms into a completely different problem, admitting different symmetries. The question then is: how can you determine those particular values of the parameters and the corresponding different symmetries? That was the underlying subject of a recent question in Mapleprimes. The answer to those questions is relatively simple and yet not entirely obvious for most of us, motivating this post, organized briefly around one example.

     

    To reproduce the input/output below you need Maple 2019 and to have installed the Physics Updates v.449 or higher.

     

    Consider the family of Korteweg-de Vries equation for u(x, t)involving three constant parameters a, b, q. For convenience (simpler input and more readable output) use the diff_table  and declare  commands

    with(PDEtools)

    U := diff_table(u(x, t))

    pde := b*U[]*U[x]+a*U[x]+q*U[x, x, x]+U[t] = 0

    b*u(x, t)*(diff(u(x, t), x))+a*(diff(u(x, t), x))+q*(diff(diff(diff(u(x, t), x), x), x))+diff(u(x, t), t) = 0

    (1)

    declare(U[])

    ` u`(x, t)*`will now be displayed as`*u

    (2)

    This pde admits a 4-dimensional symmetry group, whose infinitesimals - for arbitrary values of the parameters a, b, q- are given by

    I__1 := Infinitesimals(pde, [u], specialize_Cn = false)

    [_xi[x](x, t, u) = (1/3)*_C1*x+_C3*t+_C4, _xi[t](x, t, u) = _C1*t+_C2, _eta[u](x, t, u) = (1/3)*((-2*b*u-2*a)*_C1+3*_C3)/b]

    (3)

    Looking at pde (1) as a nonlinear problem in u, a, b and q, it splits into four cases for some particular values of the parameter:

    pde__cases := casesplit(b*u(x, t)*(diff(u(x, t), x))+a*(diff(u(x, t), x))+q*(diff(diff(diff(u(x, t), x), x), x))+diff(u(x, t), t) = 0, parameters = {a, b, q}, caseplot)

    `========= Pivots Legend =========`

     

    p1 = q

     

    p2 = b*u(x, t)+a

     

    p3 = b

     

     

    `casesplit/ans`([diff(diff(diff(u(x, t), x), x), x) = -(b*u(x, t)*(diff(u(x, t), x))+a*(diff(u(x, t), x))+diff(u(x, t), t))/q], [q <> 0]), `casesplit/ans`([diff(u(x, t), x) = -(diff(u(x, t), t))/(b*u(x, t)+a), q = 0], [b*u(x, t)+a <> 0]), `casesplit/ans`([u(x, t) = -a/b, q = 0], [b <> 0]), `casesplit/ans`([diff(u(x, t), t) = 0, a = 0, b = 0, q = 0], [])

    (4)

    The legend above indicates the pivots and the tree of cases, depending on whether each pivot is equal or different from 0. At the end there is the algebraic sequence of cases. The first case is the general case, for which the symmetry infinitesimals were computed as I__1 above, but clearly the other three cases admit more general symmetries. Consider for instance the second case, pass the ignoreparameterizingequations to ignore the parameterizing equation q = 0, and you get

    I__2 := Infinitesimals(pde__cases[2], ignore)

    `* Partial match of  'ignore' against keyword 'ignoreparameterizingequations'`

     

    [_xi[x](x, t, u) = _F3(x, t, u), _xi[t](x, t, u) = Intat(((b*u+a)*(D[1](_F3))(_a, ((b*u+a)*t-x+_a)/(b*u+a), u)-_F1(u, ((b*u+a)*t-x)/(b*u+a))*b+(D[2](_F3))(_a, ((b*u+a)*t-x+_a)/(b*u+a), u))/(b*u+a)^2, _a = x)+_F2(u, ((b*u+a)*t-x)/(b*u+a)), _eta[u](x, t, u) = _F1(u, ((b*u+a)*t-x)/(b*u+a))]

    (5)

    These infinitesimals are indeed much more general than I__1, in fact so general that (5) is almost unreadable ... Specialize the three arbitrary functions into something "easy" just to be able follow - e.g. take _F1 to be just the + operator, _F2 the * operator and _F3 = 1

    eval(I__2, [_F1 = `+`, _F2 = `*`, _F3 = 1])

    [_xi[x](x, t, u) = 1, _xi[t](x, t, u) = Intat(-(u+((b*u+a)*t-x)/(b*u+a))*b/(b*u+a)^2, _a = x)+u*((b*u+a)*t-x)/(b*u+a), _eta[u](x, t, u) = u+((b*u+a)*t-x)/(b*u+a)]

    (6)

    simplify(value([_xi[x](x, t, u) = 1, _xi[t](x, t, u) = Intat(-(u+((b*u+a)*t-x)/(b*u+a))*b/(b*u+a)^2, _a = x)+u*((b*u+a)*t-x)/(b*u+a), _eta[u](x, t, u) = u+((b*u+a)*t-x)/(b*u+a)]))

    [_xi[x](x, t, u) = 1, _xi[t](x, t, u) = (b^3*t*u^4+((3*a*t-x)*u^3-u^2*x-t*x*u)*b^2+((3*a^2*t-2*a*x)*u^2-a*u*x-a*t*x+x^2)*b+a^2*u*(a*t-x))/(b*u+a)^3, _eta[u](x, t, u) = (b*u^2+(b*t+a)*u+a*t-x)/(b*u+a)]

    (7)

    This symmetry is of course completely different than [_xi[x](x, t, u) = (1/3)*_C1*x+_C3*t+_C4, _xi[t](x, t, u) = _C1*t+_C2, _eta[u](x, t, u) = ((-2*b*u-2*a)*_C1+3*_C3)/(3*b)]computed for the general case.

     

    The symmetry (7) can be verified against pde__cases[2] or directly against pde after substituting q = 0.

    [_xi[x](x, t, u) = (1/3)*_C1*x+_C3*t+_C4, _xi[t](x, t, u) = _C1*t+_C2, _eta[u](x, t, u) = (1/3)*((-2*b*u-2*a)*_C1+3*_C3)/b]

    (8)

    SymmetryTest([_xi[x](x, t, u) = 1, _xi[t](x, t, u) = (b^3*t*u^4+((3*a*t-x)*u^3-u^2*x-t*x*u)*b^2+((3*a^2*t-2*a*x)*u^2-a*u*x-a*t*x+x^2)*b+a^2*u*(a*t-x))/(b*u+a)^3, _eta[u](x, t, u) = (b*u^2+(b*t+a)*u+a*t-x)/(b*u+a)], pde__cases[2], ignore)

    `* Partial match of  'ignore' against keyword 'ignoreparameterizingequations'`

     

    {0}

    (9)

    SymmetryTest([_xi[x](x, t, u) = 1, _xi[t](x, t, u) = (b^3*t*u^4+((3*a*t-x)*u^3-u^2*x-t*x*u)*b^2+((3*a^2*t-2*a*x)*u^2-a*u*x-a*t*x+x^2)*b+a^2*u*(a*t-x))/(b*u+a)^3, _eta[u](x, t, u) = (b*u^2+(b*t+a)*u+a*t-x)/(b*u+a)], subs(q = 0, pde))

    {0}

    (10)

    Summarizing: "to split PDE symmetries into cases according to the values of the PDE parameters, split the PDE into cases with respect to these parameters (command PDEtools:-casesplit ) then compute the symmetries for each case"

     

    Parameter continuous symmetry transformations

     

    A different, however closely related question, is whether pde admits "symmetries with respect to the parameters a, b and q", so whether exists continuous transformations of the parameters a, b and q that leave pde invariant in form.

     

    Beforehand, note that since the parameters are constants with regards to the dependent and independent variables (here u(x, t)), such continuous symmetry transformations cannot be used directly to compute a solution for pde. They can, however, be used to reduce the number of parameters. And in some contexts, that is exactly what we need, for example to entirely remove the splitting into cases due to their presence, or to proceed applying a solving method that is valid only when there are no parameters (frequently the case when computing exact solutions to "PDE & Boundary Conditions").

     

    To compute such "continuous symmetry transformations of the parameters" that leave pde invariant one can always think of these parameters as "additional independent variables of pde". In terms of formulation, that amounts to replacing the dependency in the dependent variable, i.e. replace u(x, t) by u(x, t, a, b, q)

     

    pde__xtabq := subs((x, t) = (x, t, a, b, q), pde)

    b*u(x, t, a, b, q)*(diff(u(x, t, a, b, q), x))+a*(diff(u(x, t, a, b, q), x))+q*(diff(diff(diff(u(x, t, a, b, q), x), x), x))+diff(u(x, t, a, b, q), t) = 0

    (11)

    Compute now the infinitesimals: note there are now three additional ones, related to continuous transformations of "a,b,"and q - for readability, avoid displaying the redundant functionality x, t, a, b, q, u on the left-hand sides of these infinitesimals

    Infinitesimals(pde__xtabq, displayfunctionality = false)

    [_xi[x] = (1/3)*(_F4(a, b, q)*q+_F3(a, b, q))*x/q+_F6(a, b, q)*t+_F7(a, b, q), _xi[t] = _F4(a, b, q)*t+_F5(a, b, q), _xi[a] = _F1(a, b, q), _xi[b] = _F2(a, b, q), _xi[q] = _F3(a, b, q), _eta[u] = (1/3)*((b*u+a)*_F3(a, b, q)-2*((b*u+a)*_F4(a, b, q)+(3/2)*u*_F2(a, b, q)+(3/2)*_F1(a, b, q)-(3/2)*_F6(a, b, q))*q)/(b*q)]

    (12)

    This result is more general than what is convenient for algebraic manipulations, so specialize the seven arbitrary functions of a, b, q and keep only the first symmetry that result from this specialization: that suffices to illustrate the removal of any of the three parameters a, b, or q

    S := Library:-Specialize_Fn([_xi[x] = (1/3)*(_F4(a, b, q)*q+_F3(a, b, q))*x/q+_F6(a, b, q)*t+_F7(a, b, q), _xi[t] = _F4(a, b, q)*t+_F5(a, b, q), _xi[a] = _F1(a, b, q), _xi[b] = _F2(a, b, q), _xi[q] = _F3(a, b, q), _eta[u] = (1/3)*((b*u+a)*_F3(a, b, q)-2*((b*u+a)*_F4(a, b, q)+(3/2)*u*_F2(a, b, q)+(3/2)*_F1(a, b, q)-(3/2)*_F6(a, b, q))*q)/(b*q)])[1 .. 1]

    [_xi[x] = 0, _xi[t] = 0, _xi[a] = 1, _xi[b] = 0, _xi[q] = 0, _eta[u] = -1/b]

    (13)

    To remove the parameters, as it is standard in the symmetry approach, compute a transformation to canonical coordinates, with respect to the parameter a. That means a transformation that changes the list of infinitesimals, or likewise its infinitesimal generator representation,

    InfinitesimalGenerator(S, [u(x, t, a, b, q)])

    proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc

    (14)

    into [_xi[x] = 0, _xi[t] = 0, _xi[a] = 1, _xi[b] = 0, _xi[q] = 0, _eta[u] = 0] or its equivalent generator representation  proc (f) options operator, arrow; diff(f, a) end proc

    That same transformation, when applied to pde__xtabq, entirely removes the parameter a.

    The transformation is computed using CanonicalCoordinates and the last argument indicates the "independent variable" (in our case a parameter) that the transformation should remove. We choose to remove a

    CanonicalCoordinates(S, [u(x, t, a, b, q)], [upsilon(xi, tau, alpha, beta, chi)], a)

    {alpha = a, beta = b, chi = q, tau = t, xi = x, upsilon(xi, tau, alpha, beta, chi) = (b*u(x, t, a, b, q)+a)/b}

    (15)

    declare({alpha = a, beta = b, chi = q, tau = t, xi = x, upsilon(xi, tau, alpha, beta, chi) = (b*u(x, t, a, b, q)+a)/b})

    ` u`(x, t, a, b, q)*`will now be displayed as`*u

     

    ` upsilon`(xi, tau, alpha, beta, chi)*`will now be displayed as`*upsilon

    (16)

    Invert this transformation in order to apply it

    solve({alpha = a, beta = b, chi = q, tau = t, xi = x, upsilon(xi, tau, alpha, beta, chi) = (b*u(x, t, a, b, q)+a)/b}, {a, b, q, t, x, u(x, t, a, b, q)})

    {a = alpha, b = beta, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*beta-alpha)/beta}

    (17)

    The next step is not necessary, but just to understand how all this works, verify its action over the infinitesimal generator proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc

    ChangeSymmetry({a = alpha, b = beta, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*beta-alpha)/beta}, proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi])

    proc (f) options operator, arrow; diff(f, alpha) end proc

    (18)

    Now that we see the transformation (17) is the one we want, just use it to change variables in pde__xtabq

    PDEtools:-dchange({a = alpha, b = beta, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*beta-alpha)/beta}, pde__xtabq, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi], simplify)

    upsilon(xi, tau, alpha, beta, chi)*(diff(upsilon(xi, tau, alpha, beta, chi), xi))*beta+chi*(diff(diff(diff(upsilon(xi, tau, alpha, beta, chi), xi), xi), xi))+diff(upsilon(xi, tau, alpha, beta, chi), tau) = 0

    (19)

    As expected, this result depends only on two parameters, beta, and chi, and the one equivalent to a (that is alpha, see the transformation used (17)), is not present anymore.

    To remove b or q we use the same steps, (15), (17) and (19), just changing the parameter to be removed, indicated as the last argument  in the call to CanonicalCoordinates . For example, to eliminate b (represented in the new variables by beta), input

    CanonicalCoordinates(S, [u(x, t, a, b, q)], [upsilon(xi, tau, alpha, beta, chi)], b)

    {alpha = b, beta = a, chi = q, tau = t, xi = x, upsilon(xi, tau, alpha, beta, chi) = (b*u(x, t, a, b, q)+a)/b}

    (20)

    solve({alpha = b, beta = a, chi = q, tau = t, xi = x, upsilon(xi, tau, alpha, beta, chi) = (b*u(x, t, a, b, q)+a)/b}, {a, b, q, t, x, u(x, t, a, b, q)})

    {a = beta, b = alpha, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*alpha-beta)/alpha}

    (21)

    PDEtools:-dchange({a = beta, b = alpha, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*alpha-beta)/alpha}, pde__xtabq, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi], simplify)

    upsilon(xi, tau, alpha, beta, chi)*(diff(upsilon(xi, tau, alpha, beta, chi), xi))*alpha+chi*(diff(diff(diff(upsilon(xi, tau, alpha, beta, chi), xi), xi), xi))+diff(upsilon(xi, tau, alpha, beta, chi), tau) = 0

    (22)

    and as expected this result does not contain "beta. "To remove a second parameter, the whole cycle is repeated starting with computing infinitesimals, for instance for (22). Finally, the case of function parameters is treated analogously, by considering the function parameters as additional dependent variables instead of independent ones.

     


     

    Download How_to_split_symmetries_into_cases_(II).mw

    Edgardo S. Cheb-Terrab
    Physics, Differential Equations and Mathematical Functions, Maplesoft

    Application developed using Maple and MapleSim. You can observe the vector analysis using Maple and the simulation using MapleSim. Also included a video of the result. It is a simple structure. A pole fastened by two cables and a force applied to the top. The results are to calculate tensions one and two. Consider the mass of each rope. In spanish.

    POSTE_PARADO.zip

    Lenin Araujo Castillo

    Ambassador of Maple

     

    I have noticed that there exists a Stack Exchange site for mathematica, and not for maple. My discussions with the part of Stack Exchange that handle the creation of a new Stack Exchange community have said that I must accrue a certain level of interest in the subject in order for it to be approved, and so I thought I would begin here to see if there is suffice level of interest.

    This would not diminish the use of the Maple Primes forum, and an additional proposal, in consideration of the years of dedication that have gone into this domain, be to pool the data between the two, make reputation points the same on both, perhaps even user profiles and questions answered already linkable, and all of the questions already addressed here showing up in the search on both domains.

    I am proposing this simply because I want to encourage the use of maple, and have noted that Stack Exchange is very popular. 

    So I am posting this to get overall feedback from other Maple users, as to what their opinion is regarding this proposal, and advice on whether it should and how it ought to be pursued.

    1 2 3 4 5 6 7 Last Page 2 of 269