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 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.

    Integral Transforms (revamped) and PDEs

     

    Integral transforms, implemented in Maple as the inttrans  package, are special integrals that appear frequently in mathematical-physics and that have remarkable properties. One of the main uses of integral transforms is for the computation of exact solutions to ordinary and partial differential equations with initial/boundary conditions. In Maple, that functionality is implemented in dsolve/inttrans  and in pdsolve/boundary conditions .

     

    During the last months, we have been working heavily on several aspects of these integral transform functions and this post is about that. This is work in progress, in collaboration with Katherina von Bulow

     

    The integral transforms are represented by the commands of the inttrans  package:

    with(inttrans)

    [addtable, fourier, fouriercos, fouriersin, hankel, hilbert, invfourier, invhilbert, invlaplace, invmellin, laplace, mellin, savetable, setup]

    (1)

    Three of these commands, addtable, savetable, and setup (this one is new, only present after installing the Physics Updates) are "administrative" commands while the others are computational representations for integrals. For example,

    FunctionAdvisor(integral_form, fourier)

    [fourier(a, b, z) = Int(a/exp(I*b*z), b = -infinity .. infinity), MathematicalFunctions:-`with no restrictions on `(a, b, z)]

    (2)

    FunctionAdvisor(integral_form, mellin)

    [mellin(a, b, z) = Int(a*b^(z-1), b = 0 .. infinity), MathematicalFunctions:-`with no restrictions on `(a, b, z)]

    (3)

    For all the integral transform commands, the first argument is the integrand, the second one is the dummy integration variable of a definite integral and the third one is the evaluation point. (also called transform variable). The integral representation is also visible using the convert network

    laplace(f(t), t, s); % = convert(%, Int)

    laplace(f(t), t, s) = Int(f(t)*exp(-s*t), t = 0 .. infinity)

    (4)

    Having in mind the applications of these integral transforms to compute integrals and exact solutions to PDE with boundary conditions, five different aspects of these transforms received further development:

    • 

    Compute Derivatives: Yes or No

    • 

    Numerical Evaluation

    • 

    Two Hankel Transform Definitions

    • 

    More integral transform results

    • 

    Mellin and Hankel transform solutions for Partial Differential Equations with boundary conditions


    The project includes having all these tranforms available at user level (not ready), say as FourierTransform for inttrans:-fourier, so that we don't need to input with(inttrans) anymore. Related to these changes we also intend to have Heaviside(0) not return undefined anymore, and return itself instead, unevaluated, so that one can set its value according to the problem/preferred convention (typically 0, 1/2 or 1) and have all the Maple library following that choice.

    The material presented in the following sections is reproducible already in Maple 2019 by installing the latest Physics Updates (v.435 or higher),

    Compute derivatives: Yes or No.

     

    For historical reasons, previous implementations of these integral transform commands did not follow a standard paradigm of computer algebra: "Given a function f(x), the input diff(f(x), x) should return the derivative of f(x)". The implementation instead worked in the opposite direction: if you were to input the result of the derivative, you would receive the derivative representation. For example, to the input laplace(-t*f(t), t, s) you would receive d*laplace(f(t), t, s)/ds. This is particularly useful for the purpose of using integral transforms to solve differential equations but it is counter-intuitive and misleading; Maple knows the differentiation rule of these functions, but that rule was not evident anywhere. It was not clear how to compute the derivative (unless you knew the result in advance).

     

    To solve this issue, a new command, setup, has been added to the package, so that you can set "whether or not" to compute derivatives, and the default has been changed to computederivatives = true while the old behavior is obtained only if you input setup(computederivatives = false). For example, after having installed the Physics Updates,

    Physics:-Version()

    `The "Physics Updates" version in the MapleCloud is 435 and is the same as the version installed in this computer, created 2019, October 1, 12:46 hours, found in the directory /Users/ecterrab/maple/toolbox/2019/Physics Updates/lib/`

    (1.1)

    the current settings can be queried via

    setup(computederivatives)

    computederivatives = true

    (1.2)

    and so differentiating returns the derivative computed

    (%diff = diff)(laplace(f(t), t, s), s)

    %diff(laplace(f(t), t, s), s) = -laplace(f(t)*t, t, s)

    (1.3)

    while changing this setting to work as in previous releases you have this computation reversed: you input the output (1.3) and you get the corresponding input

    setup(computederivatives = false)

    computederivatives = false

    (1.4)

    %diff(laplace(f(t), t, s), s) = -laplace(t*f(t), t, s)

    %diff(laplace(f(t), t, s), s) = diff(laplace(f(t), t, s), s)

    (1.5)

    Reset the value of computederivatives

    setup(computederivatives = true)

    computederivatives = true

    (1.6)

    %diff(laplace(f(t), t, s), s) = -laplace(t*f(t), t, s)

    %diff(laplace(f(t), t, s), s) = -laplace(f(t)*t, t, s)

    (1.7)

    In summary: by default, derivatives of integral transforms are now computed; if you need to work with these derivatives as in  previous releases, you can input setup(computederivatives = false). This setting can be changed any time you want within one and the same Maple session, and changing it does not have any impact on the performance of intsolve, dsolve and pdsolve to solve differential equations using integral transforms.

    ``

    Numerical Evaluation

     

    In previous releases, integral transforms had no numerical evaluation implemented. This is in the process of changing. So, for example, to numerically evaluate the inverse laplace transform ( invlaplace  command), three different algorithms have been implemented: Gaver-Stehfest, Talbot and Euler, following the presentation by Abate and Whitt, "Unified Framework for Numerically Inverting Laplace Transforms", INFORMS Journal on Computing 18(4), pp. 408–421, 2006.

     

    For example, consider the exact solution to this partial differential equation subject to initial and boundary conditions

    pde := diff(u(x, t), x) = 4*(diff(u(x, t), t, t))

    iv := u(x, 0) = 0, u(0, t) = 1