Applications, Examples and Libraries

Mass, spring, conveyor belt, and dry friction

by: Rouben Rostamian 7711 Product: Maple 2019

December 19 2019

12 13

Here is an animation of a mass-spring system where the mass slides horizontally on a steadily moving conveyor belt.

The contact between the block of mass and the belt is of the dry friction kind (Coulomb's friction). Consequently the block periodically sticks to the belt and moves forward with it until the force of the stretching spring overcomes the force of friction and yanks it back, making it to slip against the belt. In the animation the block is shown in a dark color while slipping, and a light color while sticking.

The fully executed Maple worksheet can be slow to load and requires a good deal of memory. Therefore I have attached two versions which are identical in all respects except that in one of them I have removed the Maple output to make is easy to load if your computer has limitations.

Download worksheet (no Maple output) block-sliding-on-conveyor-belt-stripped.mw

Doiwnload worksheet (with Maple output) (sorry, exceeds MaplePrime's size limit)

The solution to another problem of Putnam Mathematical...

by: Kitonum 20084

December 12 2019

0 0

Here is two solutions with Maple of the problem A2 of Putnam Mathematical Competition 2019 . The first solution is entirely based on the use of the geometry package; the second solution does not use this package. Since the triangle is defined up to similarity, without loss of generality, we can set its vertices A(0,0) , B(1,0) , C(x0,y0) and then calculate the parameters x0, y0 using the conditions of the problem.

The problem

A2: In the triangle ∆ABC, let G be the centroid, and let I be the center of the
inscribed circle. Let α and β be the angles at the vertices A and B, respectively.
Suppose that the segment IG is parallel to AB and that β = 2*arctan(1/3). Find α.

>

# Solution 1 with the geometry package
restart;

# Calculation

with(geometry):
local I:
point(A,0,0): point(B,1,0): point(C,x0,y0):
assume(y0>0,-y0*(-1+x0-((1-x0)^2+y0^2)^(1/2))+y0*((x0^2+y0^2)^(1/2)+x0) <> 0):
triangle(t,[A,B,C]):
incircle(ic,t, 'centername'=I):
Cn:=coordinates(I):
centroid(G,t):
CG:=coordinates(G):
a:=-expand(tan(2*arctan(1/3))):
solve({Cn[2]=CG[2],y0/(x0-1)=a}, explicit);
point(C,eval([x0,y0],%)):
answer=FindAngle(line(AB,[A,B]),line(AC,[A,C]));

# Visualization (G is the point of medians intersection)

triangle(t,[A,B,C]):
incircle(ic,t, 'centername'=I):
centroid(G,t):
segment(s,[I,G]):
median(mB,B,t): median(mC,C,t):
draw([A(symbol=solidcircle,color=black),B(symbol=solidcircle,color=black),C(symbol=solidcircle,color=black),I(symbol=solidcircle,color=green),G(symbol=solidcircle,color=blue),t(color=black),s(color=red,thickness=2),ic(color=green),mB(color=blue,thickness=0),mC(color=blue,thickness=0)], axes=none, size=[800,500], printtext=true,font=[times,20]);

Warning, The imaginary unit, I, has been renamed _I

Warning, solve may be ignoring assumptions on the input variables.

>

# Solution 2 by a direct calculation

# Calculation

restart;
local I;
sinB:=y0/sqrt(x0^2+y0^2):
cosB:=x0/sqrt(x0^2+y0^2):
Sol1:=eval([x,y],solve({y=-(x-1)/3,y=(sinB/(1+cosB))*x}, {x,y})):
tanB:=expand(tan(2*arctan(1/3))):
Sol2:=solve({y0/3=Sol1[2],y0=-tanB*(x0-1)},explicit);
A:=[0,0]: B:=[1,0]: C:=eval([x0,y0],Sol2[2]):
AB:=<(B-A)[]>: AC:=<(C-A)[]>:
answer=arccos(AB.AC/sqrt(AB.AB)/sqrt(AC.AC));

# Visualization

with(plottools): with(plots):
ABC:=curve([A,B,C,A]):
I:=simplify(eval(Sol1,Sol2[2]));
c:=circle(I,eval(Sol1[2],Sol2[2]),color=green):
G:=(A+B+C)/~3;
IG:=line(I,G,color=red,thickness=2):
P:=pointplot([A,B,C,I,G], color=[black$3,green,blue], symbol=solidcircle):
T:=textplot([[A[],"A"],[B[],"B"],[C[],"C"],[I[],"I"],[G[],"G"]], font=[times,20], align=[left,below]):
M:=plot([[(C+t*~((A+B)/2-C))[],t=0..1],[(B+t*~((A+C)/2-B))[],t=0..1]], color=blue, thickness=0):
display(ABC,c,IG,P,T,M, scaling=constrained, axes=none,size=[800,500]);

Warning, The imaginary unit, I, has been renamed _I

>

Download Putnam.mw

Putnam Mathematical Competition 2019

by: vv 12453 Product: Maple 2018

December 11 2019

5 2

Maple can easily solve the B4 problem of the Putnam Mathematical Competition 2019 link

B4. Let F be the set of functions f(x,y) that are twice continuously differentiable for x≥1, y≥1 and that satisfy the following two equations:

For each f2F, let

Determine m(f), and show that it is independent of the choice of f.

>	# Solution

>	pdsolve({ xdiff(f(x,y),x)+ydiff(f(x,y),y) = xyln(xy), x^2diff(f(x,y),x,x)+y^2diff(f(x,y),y,y) = xy });

${f(x, y) = (1/2)*(x*y+2*_C1)*ln(x*y)-(1/2)*x*y-2*_C1*ln(x)+_C2}$

(1)

>	f:=unapply(rhs(%[]), x,y);

(2)

>	h := f(s+1, s+1) - f(s+1, s) - f(s, s+1) + f(s, s);

(1/2)*((s+1)^2+2*_C1)*ln((s+1)^2)-(1/2)*(s+1)^2-(s*(s+1)+2*_C1)*ln(s*(s+1))+s*(s+1)+(1/2)*(s^2+2*_C1)*ln(s^2)-(1/2)*s^2

(3)

>	minimize(h, s=1..infinity);

(4)

>	answer = simplify(%);

(5)

>

Download putnam2019-b4.mw

Triangulation, hatching and texturing of simple...

by: mmcdara 6149 Product: Maple

December 09 2019

7 2

Hi,

As an amusement, I decided several months ago to develop some procedures to fill a simple polygon* by hatches or simple textures.

* A simple polygon is a polygon whose sides either do not intersect or either have a common vertex.

This work started with the very simple observation that if I was capable to hatch or texture a triangle, I will be too to hatch or texture any simple polygon once triangulated.
I also did some work to extend this work to non-simple polygons but there remains some issues to fix (which explains while it is not deliverd here).

The main ideat I used for hatching and texturing is based upon the description of each triangles by a set of 3 inequalities that any interior point must verify.
A hatch of this triangle is thius a segment whose each point is interior.
The closely related idea is used for texturing. Given a simple polygon, periodically replicated to form the texture, the set of points of each replicate that are interior to a given triangle must verify a set of inequalities (the 3 that describe the triangle, plus N if the pattern of the texture is a simple polygon with N sides).

Unfortunately I never finalise this work.
Recently @Christian Wolinski asked a question about texturing that reminded me this "ancient" work of mine.
So I decided to post it as it is, programatically imperfect, lengthy to run, and most of all french-written for a large part.
I guess it is a quite unorthodox way to proceed but some here could be interested by this work to take it back and improve it.

The module named "trianguler" (= triangulate) is a translation into Maple of Frederic Legrand's Python code (full reference given in the worksheet).
I added my own procedure "hachurer" (= hatching) to this module.
The texturing part is not included in this module for it was still in development.

A lot of improvements can be done that I could list, but I prefer not to be too intrusive in your evaluation of this work. So make your own idea about it and do not hesitate to ask me any informations you might need (in particular translation questions).

PS: this work has been done with Maple 2015.2

>

restart:

Reference: http://www.f-legrand.fr/scidoc/docmml/graphie/geometrie/polygone/polygone.html
                    (in french)
                    reference herein : M. de Berg, O. Cheong, M. van Kreveld, M. Overmars,
                                                 Computational geometry, (Springer, 2010)

Direct translation of the Frederic Legrand's Python code

Meaning of the different french terms

voisin_sommet (n, i, di)
        let L the list [1, ..., n] where n is the number of vertices
        This procedure returns the index of the neighbour (voisin) of the vertex (sommet) i when L is rotated by di

equation_droite (P0, P1, M)
        Let P0 and P1 two vertices and M an arbitrary point.
        Let (P0, P1) the vector originated at P0 and ending at P1 (idem for (P0, M)) and the unitary vector in the Z direction.
        This procedure returns (P0, P1) o (P0, M) .

point_dans_triangle (triangle, M) P1, P2]
        This procedure returns "true" if point M is within (strictly) the triangle "triangle" and "false" if not.

sommet_distance_maximale (polygone, P0, P1, P2, indices)
        Given a polygon (polygone) threes vertices P0, P1 and P2 and a list of indices , this procedure returns
        the vertex of the polygon "polygone" which verifies: 1/ this vertex is strictly within
        the triangle [P0, P1, P2] and 2/ it is the farthest from side [P1, P2] amid all the vertices that verifies point 1/.
        If there is no such vertex the procedure returns NULL.

sommet_gauche (polygone)
        This procedure returns the index of the leftmost ("gauche" means "left") vertex in polygon "polygone".
        If more than one vertices have the same minimum abscissa value then only the first one is returned.

nouveau_polygone(polygone,i_debut,i_fin)
        This procedure creates a new polygon from index i_debut (could be translated by i_first) to i_end (i_last)

trianguler_polygone_recursif(polygone)
        This procedure recursively divides a polygon in two parts A and B from its leftmost vertex.
         If A (B) is a triangle the list "liste_triangles" (mening "list of triangles") is augmented by A (B);
         if not the procedure executes recursively on A and B.

trianguler_polygone(polygone)
         This procedure triangulates the polygon "polygon"

hachurer(shapes, hatch_angle, hatch_number, hatch_color)
         This procedure generates stes of hatches of different angles, colors and numbers

Limitations:
   1/ "polygone" is a simply connected polygon
   2/ two different sides S and S', either do not intersect or either have a common vertex

>

trianguler := module()
export voisin_sommet, equation_droite, interieur_forme, point_dans_triangle, sommet_distance_maximale,
       sommet_gauche, nouveau_polygone, trianguler_polygone_recursif, trianguler_polygone, hachurer:

#-------------------------------------------------------------------
voisin_sommet := (n, i, di) -> ListTools:-Rotate([$1..n], di)[i]:

#-------------------------------------------------------------------
equation_droite := proc(P0, P1, M) (P1[1]-P0[1])*(M[2]-P0[2]) - (P1[2]-P0[2])*(M[1]-P0[1]) end proc:

#-------------------------------------------------------------------
interieur_forme := proc(forme, M)
  local N:
  N := numelems(forme);
  { seq( equation_droite(forme[n], forme[piecewise(n=N, 1, n+1)], M) >= 0, n=1..N) }
end proc:

#-------------------------------------------------------------------
point_dans_triangle := proc(triangle, M)
  `and`(
          is( equation_droite(triangle[1], triangle[2], M) > 0 ),
          is( equation_droite(triangle[2], triangle[3], M) > 0 ),
          is( equation_droite(triangle[3], triangle[1], M) > 0 )
       )
end proc:

#-------------------------------------------------------------------
sommet_distance_maximale := proc(polygone, P0, P1, P2, indices)
  local n, distance, j, i, M, d;

  n        := numelems(polygone):
  distance := 0:
  j        := NULL:

  for i from 1 to n do
    if `not`(member(i, indices)) then
      M := polygone[i];
      if point_dans_triangle([P0, P1, P2], M) then
        d := abs(equation_droite(P1, P2, M)):
        if d > distance then
          distance := d:
          j := i
        end if:
      end if:
    end if:
  end do:
  return j:
end proc:

#-------------------------------------------------------------------
sommet_gauche := polygone -> sort(polygone, key=(x->x[1]), output=permutation)[1]:

#-------------------------------------------------------------------
nouveau_polygone := proc(polygone, i_debut, i_fin)
  local n, p, i:

  n := numelems(polygone):
  p := NULL:
  i := i_debut:

  while i <> i_fin do
    p := p, polygone[i]:
    i := voisin_sommet(n, i, 1)
  end do:
  p := [p, polygone[i_fin]]
end proc:

#-------------------------------------------------------------------
trianguler_polygone_recursif := proc(polygone)
  local n, j0, j1, j2, P0, P1, P2, j, polygone_1, polygone_2:
  global liste_triangles:
  n := numelems(polygone):
  j0 := sommet_gauche(polygone):
  j1 := voisin_sommet(n, j0, +1):
  j2 := voisin_sommet(n, j0, -1):
  P0 := polygone[j0]:
  P1 := polygone[j1]:
  P2 := polygone[j2]:
  j := sommet_distance_maximale(polygone, P0, P1, P2, [j0, j1, j2]):

  if `not`(j::posint) then
    liste_triangles := liste_triangles, [P0, P1, P2]:
    polygone_1      := nouveau_polygone(polygone,j1,j2):
    if numelems(polygone_1) = 3 then
      liste_triangles := liste_triangles, polygone_1:
    else
      thisproc(polygone_1)
    end if:

  else
    polygone_1 := nouveau_polygone(polygone, j0, j ):
    polygone_2 := nouveau_polygone(polygone, j , j0):
    if numelems(polygone_1) = 3 then
      liste_triangles := liste_triangles, polygone_1:
    else
      thisproc(polygone_1)
    end if:
    if numelems(polygone_2) = 3 then
      liste_triangles := liste_triangles, polygone_2:
    else
      thisproc(polygone_2)
    end if:
  end if:

  return [liste_triangles]:
end proc:

#-------------------------------------------------------------------
trianguler_polygone := proc(polygone)
  trianguler_polygone_recursif(polygone):
  return liste_triangles:
end proc:

#-------------------------------------------------------------------
hachurer := proc(shapes, hatch_angle::list, hatch_number::list, hatch_color::list)

local A, La, Lp;
local N, P, _sides, L_sides, Xshape, ch, rel, p_rel, n, sol, p_range:
local AllHatches, window, p, _segment:
local NT, ka, N_Hatches, p_range_t, nt, shape, p_hatches, WhatHatches:

#-----------------------------------------------------------------
# internal functions:
#
# La(x, y, alpha, p) is the implicit equation of a straight line of angle alpha relatively
#                    to the horizontal axis and intercept p
#
# Lp(x, y, P) is the implicit equation of a straight line passing through points P[1] and P[2]
#
# interieur_triangle(triangle, M)

La := (x, y, alpha, p) -> cos(alpha)*x - sin(alpha)*y + p;
Lp := proc(x, y, P::list) (x-P[1][1])*(P[2][2]-P[1][2]) - (y-P[1][2])*(P[2][1]-P[1][1] = 0) end proc;

p_range    := [+infinity, -infinity]:
NT         := numelems(shapes):

AllHatches := NULL:

for ka from 1 to numelems(hatch_angle) do
  A         := hatch_angle[ka]:
  N_Hatches := hatch_number[ka]:
  p_range_t := NULL:
  _sides    := []:
  L_sides   := []:
  rel       := []:
  for nt from 1 to NT do

    shape := shapes[nt]:
    # _sides : two points description of the sides of the shape
    # L_sides : implicit equations of the straight lines that support the sides

    N        := [1, 2, 3];
    P        := [2, 3, 1];
    _sides   := [ _sides[] , [ seq([shape[n], shape[P[n]]], n in N) ] ];
    L_sides := [ L_sides[], Lp~(x, y, _sides[-1]) ];

    # Inequalities that define the interior of the shape

    rel := [ rel[], trianguler:-interieur_forme(shape, [x, y]) ];

    # Given the orientation of the hatches we search here the extreme values of
    # the intercept p for wich a straight line of equation La(x, y, alpha, p)
    # cuts the shape.

    p_rel := NULL:

    for n from 1 to numelems(L_sides[-1]) do
      sol   := solve({La(x, y, A, q), lhs(L_sides[-1][n])} union rel[-1], [x, y]);
      p_rel := p_rel, `if`(sol <> [], [rhs(op(1, %)), rhs(op(3, %))], [+infinity, -infinity]);
    end do:
    p_range_t := p_range_t, evalf(min(op~(1, [p_rel]))..max(op~(2, [p_rel])));
    p_range   := evalf(min(op~(1, [p_rel]), op(1, p_range))..max(op~(2, [p_rel]), op(2, p_range)));

  end do: # end of the loop over triangles

  p_range_t := [p_range_t]:
  p_hatches := [seq(p_range, (op(2, p_range)-op(1, p_range))/N_Hatches)]:
  # Building of the hatches
  #
  # This construction is far from being optimal.
  # Here again the main goal was to obtain the hatches with a minimum effort
  # if algorithmic development.

  window      := min(op~(1..shape))..max(op~(1..shape)):
  WhatHatches := map(v -> map(u -> if verify(u, v, 'interval'('closed') ) then u end if, p_hatches), p_range_t):

  for nt from 1 to NT do
    for p in WhatHatches[nt] do
      _segment := []:
      for n from 1 to numelems(L_sides[nt]) do
         _segment := _segment, evalf( solve({La(x, y, A, p), lhs(L_sides[nt][n])} union rel[nt], [x, y]) );
      end do;
      map(u -> u[], [_segment]);
      AllHatches := AllHatches, plot(map(u -> rhs~(u), %), color=hatch_color[ka]):
    end do:
  end do;

end do: # end of loop over hatch angles

plots:-display(
  PLOT(POLYGONS(polygone, COLOR(RGB, 1$3))),
  AllHatches,
  scaling=constrained
)

end proc:

end module:

Legrand's example (see reference above)

>	global liste_triangles: liste_triangles := NULL:

>	polygone := [[0,0],[0.5,-1],[1.5,-0.2],[2,-0.5],[2,0],[1.5,1],[0.3,0],[0.5,1]]: trianguler:-trianguler_polygone(polygone): PLOT(seq(POLYGONS(u, COLOR(RGB, rand()/10^12, rand()/10^12, rand()/10^12)), u in liste_triangles), VIEW(0..2, -2..2))

>	trianguler:-hachurer([liste_triangles], [-Pi/4, Pi/4], [40, 40], [red, blue])

>	F := (P, a, b) -> map(p -> [p[1]+a, p[2]+b], P):

>

MOTIF := [[0, 0], [0.05, 0], [0.05, 0.05], [0, 0.05]];
motifs := [ seq(seq(F(MOTIF, 0+i*0.075, 0+j*0.075), i=0..26), j=-14..13) ]:

plots:-display(
  plot([polygone[], polygone[1]], color=red, filled),
  map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):

texture    := NULL:
rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):

for ref in liste_triangles do
  ref;
  #
  # the three lines below are used to define REF counter clockwise
  #
  g           := trianguler:-sommet_gauche(ref):
  bas         := sort(op~(2, ref), output=permutation);
  REF         := ref[[g, op(map(u -> if u<>g then u end if, bas))]];
  rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
  texture_ref := map(u -> plots:-inequal(rel_ref union u, x=0..2, y=-1..1, color=blue, 'nolines'), rel_motifs):
  texture     := texture, texture_ref:
end do:

plots:-display(
  plot([polygone[], polygone[1]], color=red, scaling=constrained),
  texture
)

>

MOTIF := [[0, 0], [0.05, 0], [0.05, 0.05], [0, 0.05]];
motifs := [ seq(seq(F(MOTIF, piecewise(j::odd, 0.05, 0.1)+i*0.1, 0+j*0.05), i=-0.2..20), j=-20..20) ]:
plots:-display(
  plot([polygone[], polygone[1]], color=red, filled),
  map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):

texture    := NULL:
rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):

for ref in liste_triangles do
  ref;
  g := trianguler:-sommet_gauche(ref):
  bas := sort(op~(2, ref), output=permutation);
  REF := ref[[g, op(map(u -> if u<>g then u end if, bas))]];
  rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
  texture_ref := map(u -> plots:-inequal(rel_ref union u, x=0..2, y=-1..1, color=blue, 'nolines'), rel_motifs):
  texture     := texture, texture_ref:
end do:

plots:-display(
  plot([polygone[], polygone[1]], color=red, scaling=constrained),
  texture
)

>

MOTIF := [[0, 0], [0.4, 0], [0.4, 0.14], [0, 0.14]]:
motifs := [ seq(seq(F(MOTIF, piecewise(j::odd, 0.4, 0.2)+i*0.4, 0+j*0.14), i=-1..4), j=-8..7) ]:

plots:-display(
  plot([polygone[], polygone[1]], color=red, filled),
  map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):

palettes := ColorTools:-PaletteNames():
ColorTools:-GetPalette("HTML"):

couleurs := [SandyBrown, PeachPuff, Peru, Linen, Bisque, Burlywood, Tan, AntiqueWhite,      NavajoWhite, BlanchedAlmond, PapayaWhip, Moccasin, Wheat]:

nc   := numelems(couleurs):
roll := rand(1..nc):

motifs_nb      := numelems(motifs):
motifs_couleur := [ seq(cat("HTML ", couleurs[roll()]), n=1..motifs_nb) ]:

texture    := NULL:
rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):

for ref in liste_triangles do
  ref;
  g := trianguler:-sommet_gauche(ref):
  bas := sort(op~(2, ref), output=permutation);
  REF := ref[[g, op(map(u -> if u<>g then u end if, bas))]];
  rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
  texture_ref := map(n -> plots:-inequal(rel_ref union rel_motifs[n], x=0..2, y=-1..1, color=motifs_couleur[n], 'nolines'), [$1..motifs_nb]):
  texture     := texture, texture_ref:
end do:

plots:-display(
  plot([polygone[], polygone[1]], color=red, scaling=constrained),
  texture
)

>

MOTIF := [ seq(0.1*~[cos(Pi/6+Pi/3*i), sin(Pi/6+Pi/3*i)], i=0..5) ]:
motifs := [ seq(seq(F(MOTIF, i*0.2*cos(Pi/6)+piecewise(j::odd, 0, 0.08), j*0.3*sin(Pi/6)), i=0..12), j=-6..6) ]:

plots:-display(
  plot([polygone[], polygone[1]], color=red, filled),
  map(u -> plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):

motifs_nb      := numelems(motifs):
motifs_couleur := [ seq(`if`(n::even, yellow, brown) , n=1..motifs_nb) ]:

texture    := NULL:
rel_motifs := map(u -> trianguler:-interieur_forme(u, [x, y]), motifs):

for ref in liste_triangles do
  ref;
  g := trianguler:-sommet_gauche(ref):
  bas := sort(op~(2, ref), output=permutation);
  REF := ref[[g, op(map(u -> if u<>g then u end if, bas))]];
  rel_ref     := trianguler:-interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
  texture_ref := map(n -> plots:-inequal(rel_ref union rel_motifs[n], x=0..2, y=-1..1, color=motifs_couleur[n], 'nolines'), [$1..motifs_nb]):
  texture     := texture, texture_ref:
end do:

plots:-display(
  plot([polygone[], polygone[1]], color=red, scaling=constrained),
  texture
)

>

Download Triangulation_Hatching_Texturing.mw

Gingerbread House with 3D Geometry

by: TechnicalSupport 670 Product: Maple 2019

December 05 2019

12 0

We recently had a question about using some of the plotting commands in Maple to draw things. We were feeling creative and thought why not take it a step further and draw something in 3D.

Using the geom3d, plottools, and plots packages we decided to make a gingerbread house.

To make the base of the house we decided to use 2 cubes, as these would give us additional lines and segments for the icing on the house.

point(p__1,[2,3,2]):
point(p__2,[3,3,2]):
cube(c1,p__1,2):
cube(c2,p__2,2):
base:=draw([c1,c2],color=tan);

Using the same cubes but changing the style to be wireframe and point we made some icing lines and decorations for the gingerbread house.

base_decor1:=draw([c1,c2],style=wireframe,thickness=3,color=red,transparency=0.2):
base_decor2:=draw([c1,c2],style=wireframe,thickness=10,color=green,linestyle=dot):
base_decor3:=draw([c1,c2],style=point,thickness=2,color="Silver",symbol=sphere):
base_decor:=display(base_decor1,base_decor2,base_decor3);

To create the roof we found the vertices of the cubes and used those to find the top corners of the base.

v1:=vertices(c1):
v2:=vertices(c2):
pc1:=seq(point(pc1||i,v1[i]),i=1..nops(v1)):
pc2:=seq(point(pc2||i,v2[i]),i=1..nops(v2)):
topCorners:=[pc1[5],pc1[6],pc2[1],pc2[2]]:
d1:=draw(topCorners):

Using these top corners we found the midpoints (where the peak of the roof would be) and added the roof height to the coordinates.

midpoint(lc1,topCorners[1],topCorners[2]):
detail(lc1);

point(cc1,[-(2*sqrt(3))/3 + 2, (2*sqrt(3))/3 + 3+1, 2]):
d3:=draw(cc1):

midpoint(lc2,topCorners[3],topCorners[4]):
detail(lc2);

point(cc2,[(2*sqrt(3))/3 + 3, (2*sqrt(3))/3 + 3+1, 2]):
d4:=draw(cc2):

With the midpoints and vertices at the front and rear of the house we made two triangles for the attic of the gingerbread house.

triangle(tf,[topCorners[1],topCorners[2],cc1]):
front:=draw(tf,color=brown):

triangle(tb,[topCorners[3],topCorners[4],cc2]):
back:=draw(tb,color=tan):

Using these same points again we made more triangles to be the roof.

triangle(trl,[cc1,cc2,pc1[5]]):
triangle(trh,[pc2[2],pc1[6],cc1]):
triangle(tll,[cc1,cc2,pc2[2]]):
triangle(tlh,[pc2[1],pc1[5],cc2]):
roof:=draw([trl,trh,tll,tlh],color="Chocolate");

Our gingerbread house now had four walls, a roof, and icing, but no door. Creating the door was as easy as making a parallelepiped, but what is a door without more icing?

door:=display(plottools:-parallelepiped([1,0,0],[0,1.2,0],[0,0,0.8],[0.8,1.9,1.6]),color="DarkRed"):
door_decor1:=display(plottools:-parallelepiped([1,0,0],[0,1.2,0],[0,0,0.8],[0.8,1.9,1.6]),color="Gold",style=line):
door_decor2:=display(plottools:-parallelepiped([1,0,0],[0,1.2,0],[0,0,0.8],[0.8,1.9,1.6]),color="Silver", style=line,linestyle=dot,thickness=5):
door_decor:=display(door_decor1,door_decor2):

Now having a door we could have left it like this, but what better way to decorate a gingerbread house than with candy canes? Naturally, if we were going to have one candy cane we were going to have lots of candy canes. To facilitate this we made a candy cane procedure.

candy_pole:=proc(c:=[0,0,0], {segR:=0.1}, {segH:=0.1}, {segn:=7}, {tilt_theta:=0}, {theta:=0}, {arch:=true}, {flip:=false})
local cane1,cane2,cane_s,cane_c,cane0,cane,i,cl,cd,ch, cane_a,tmp,cane_ac,cane_a1,cane00,cane01,cane02,cane_a1s,tmp2,cane_a2s:
uses plots,geom3d:
cl:=c[1]:
cd:=c[2]:
ch:=c[3]:
cane1:=plottools:-cylinder([cd, ch, cl], segR, segH,style=surface):

cane2:=display(plottools:-rotate(cane1,Pi/2,[[cd,ch,cl],[cd+1,ch,cl]])):
cane_s:=[cane2,seq(display(plottools:-translate(cane2,0,segH*i,0)),i=1..segn-1)]:
cane_c:=seq(ifelse(type(i,odd),red,white),i=1..segn):

cane0:=display(cane_s,color=[cane_c]):

if arch then

cane_a:=plottools:-translate(cane2,0,segH*segn-segH/2,0):
tmp:=i->plottools:-rotate(cane_a,i*Pi/24, [ [cd,ch+segH*segn-segH/2,cl+segR*2] , [cd+1,ch+segH*segn-segH/2,cl+segR*2] ] ):

cane_ac:=seq(ifelse(type(i,odd),red,white),i=1..24):

                cane_a1s:=seq(plottools:-translate(tmp(i),0,segH*i/12,segR*i/4),i=1..12):

tmp2:=i->plottools:-rotate(cane_a1s[12],i*Pi/24,[[cd,ch+segH*segn-0.05,cl+segR*2],[cd+1,ch+segH*segn-0.05,cl+segR*2]]):

cane_a2s:=seq(plottools:-translate(tmp2(i),0,-segH*i/500,segR*i/4),i=1..12):
cane_a1:=display(cane_a1s,cane_a2s,color=[cane_ac]):
cane00:=display(cane0,cane_a1);

                if flip then

cane01:=plottools:-rotate(cane00,tilt_theta,[[cd,ch,cl],[cd+1,ch,cl]]):
cane02:=plottools:-rotate(cane01,theta,[[cd,ch,cl],[cd,ch+1,cl]]):
cane:=plottools:-reflect(cane01,[[cd,ch,cl],[cd,ch+1,cl]]):

                else

cane01:=plottools:-rotate(cane00,tilt_theta,[[cd,ch,cl],[cd+1,ch,cl]]):
cane:=plottools:-rotate(cane01,theta,[[cd,ch,cl],[cd,ch+1,cl]]):

                end if:

                return cane:

else

                cane:=plottools:-rotate(cane0,tilt_theta,[[cd,ch,cl],[cd+1,ch,cl]]):

                return cane:

end if:

end proc:

With this procedure we decided to add candy canes to the front, back, and sides of the gingerbread house. In addition we added two candy poles.

Candy Canes in front of the house:

cane1:=candy_pole([1.2,0,2],segn=9,arch=false):
cane2:=candy_pole([2.8,0,2],segn=9,arch=false):
cane3:=candy_pole([2.7,0.8,3.3],segn=9,segR=0.05,tilt_theta=-Pi/7):
cane4:=candy_pole([1.3,0.8,3.3],segn=9,segR=0.05,tilt_theta=-Pi/7,flip=true):
front_canes:=display(cane1,cane2,cane3,cane4):

Candy Canes at the back of the house:

caneb3:=candy_pole([1.5,4.2,2.5],segn=15,segR=0.05,tilt_theta=-Pi/3,flip=true):
caneb4:=candy_pole([2.5,4.2,2.5],segn=15,segR=0.05,tilt_theta=-Pi/3):}
back_canes:=display(caneb3,caneb4):

Candy Canes at the left of the house:

canel1:=candy_pole([0.8,1.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
canel2:=candy_pole([0.8,2.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
canel3:=candy_pole([0.8,4,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
left_canes:=display(canel1,canel2,canel3):

Candy Canes at the right of the house:

caner1:=candy_pole([3.2,1.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
caner2:=candy_pole([3.2,2.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
caner3:=candy_pole([3.2,4,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
right_canes:=display(caner1,caner2,caner3):

canes:=display(front_canes,back_canes,right_canes,left_canes):

With these canes in place all that was left was to create the ground and display our Gingerbread House.

ground:=display(plottools:-parallelepiped([5,0,0],[0,0.5,0],[0,0,4],[0,1.35,0]),color="DarkGreen"):

display([door,door_decor,d1,base,base_decor,d3,d4,front,back,roof,ground,canes],orientation=[-100,0,95]);

You can download the full worksheet creating the gingerbread house below:

Geometry_Gingerbread.mw

Maple and two light sources

by: Christopher2222 5785 Product: Maple 7

November 24 2019

2 4

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

Transfer matrix method for vibration analysis

by: Adjal yassine 25 Product: Maple

November 11 2019

0 0

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

Download transfer.mw

Representing a hierarchical table as a tree graph

by: mmcdara 6149 Product: Maple

November 08 2019

2 0

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);

(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

Stiffness matrix implementation of a finite element...

by: Adjal yassine 25

November 07 2019

0 5

>

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]):

>

$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)

>

Download stiffnesmatrice.mw

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

Free vibration of a cracked composite beam

by: Adjal yassine 25 Product: Maple

November 06 2019

5 0

>	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

>	EIW^(iv)-momega^2*W=0;

(1)

The general solution form of bending vibration equation

>	W1:=A[1]cosh(mux)+A[2]sinh(mux)+A[3]cos(mux)+A[4]sin(mux);

(2)

where

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

(3)

>	mu:=((momega^2L^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]musinh(mux)+A[2]mucosh(mux)-A[3]musin(mux)+A[4]mucos(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^2cosh(mux)+A[2]mu^2sinh(mux)-A[3]mu^2cos(mux)-A[4]mu^2sin(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^3sinh(mux)+A[2]mu^3cosh(mux)+A[3]mu^3sin(mux)-A[4]mu^3cos(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(mux)+A[6]sinh(mux)+A[7]cos(mux)+A[8]sin(mux);

(7)

>

>	theta2:= (1/L)(A[5]musinh(mux)+A[6]mucosh(mux)-A[7]musin(mux)+A[8]mucos(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^2cosh(mux)+A[6]mu^2sinh(mux)-A[7]mu^2cos(mux)-A[8]mu^2sin(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^3sinh(mux)+A[6]mu^3cosh(mux)+A[7]mu^3sin(mux)-A[8]mu^3cos(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));

(11)

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

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

(12)

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

(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);

(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);

(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);

(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);

(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);

(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);

(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);

(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);

(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

(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

(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

>

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

>	EIW^(iv)-momega^2*W=0;

(1)

The general solution form of bending vibration equation

>	W1:=A[1]cosh(mux)+A[2]sinh(mux)+A[3]cos(mux)+A[4]sin(mux);

(2)

where

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

(3)

>	mu:=((momega^2L^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]musinh(mux)+A[2]mucosh(mux)-A[3]musin(mux)+A[4]mucos(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^2cosh(mux)+A[2]mu^2sinh(mux)-A[3]mu^2cos(mux)-A[4]mu^2sin(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^3sinh(mux)+A[2]mu^3cosh(mux)+A[3]mu^3sin(mux)-A[4]mu^3cos(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(mux)+A[6]sinh(mux)+A[7]cos(mux)+A[8]sin(mux);

(7)

>

>	theta2:= (1/L)(A[5]musinh(mux)+A[6]mucosh(mux)-A[7]musin(mux)+A[8]mucos(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^2cosh(mux)+A[6]mu^2sinh(mux)-A[7]mu^2cos(mux)-A[8]mu^2sin(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^3sinh(mux)+A[6]mu^3cosh(mux)+A[7]mu^3sin(mux)-A[8]mu^3cos(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));

(11)

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

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

(12)

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

(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);

(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);

(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);

(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);

(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);

(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);

(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);

(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);

(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

(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

(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

>

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

by: ecterrab 13431 Product: Maple

October 18 2019

5 6

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 involving three constant parameters . For convenience (simpler input and more readable output) use the diff_table and declare commands

>

(1)

>

(2)

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

>

(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(bu(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)$

`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 above, but clearly the other three cases admit more general symmetries. Consider for instance the second case, pass the to ignore the parameterizing equation , and you get

>

[_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 , 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

>

[_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)

>

[_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 computed for the general case.

The symmetry (7) can be verified against or directly against pde after substituting .

(8)

>

(9)

>

(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 ), 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 by

>

(11)

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

>

[_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 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

>

(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,

>

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

(14)

into or its equivalent generator representation

That same transformation, when applied to , 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

>

(15)

>

(16)

Invert this transformation in order to apply it

>

(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])

(18)

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

>	$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)$

(19)

As expected, this result depends only on two parameters, , and , and the one equivalent to a (that is , 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 ), input

>

(20)

>

(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)$

(22)

and as expected this result does not contain 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

Pole simulation with tensions in MapleSim

by: Lenin Araujo Castillo 1790 Products: Maple 2018 , MapleSim 2018

October 14 2019

2 0

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

Integral Transforms (revamped) and PDE

by: ecterrab 13431 Product: Maple

October 01 2019

8 14

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:

>

(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,

>

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

(2)

>

[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) = 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 anymore. Related to these changes we also intend to have not return 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 , the input should return the derivative of ". 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 you would receive . 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 . For example, after having installed the Physics Updates,

>

`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

>

(1.2)

and so differentiating returns the derivative computed

>

(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

>

(1.4)

>

(1.5)

Reset the value of computederivatives

>

(1.6)

>

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

>

Note that these two conditions are not entirely compatible: the solution returned cannot be valid for and simultaneously. However, a solution discarding that point does exist and is given by

>

u(x, t) = -invlaplace(exp(-(1/2)*s^(1/2)*t)/s, s, x)+1

(2.1)

Verifying the solution, one condition remains to be tested

>

[0, 0, -invlaplace(exp(-(1/2)*s^(1/2)*t)/s, s, 0)]

(2.2)

Since we now have numerical evaluation rules, we can test that what looks different from 0 in the above is actually 0.

>

-invlaplace(exp(-(1/2)*s^(1/2)*t)/s, s, 0)

(2.3)

Add a small number to the initial value of t to skip the point

>

The default method used is the method of Euler sums and the numerical evaluation is performed as usual using the evalf command. For example, consider

>

The Laplace transform of F is given by

>

(1/2)*2^(1/2)*Pi^(1/2)*exp(-(1/2)/s)/s^(3/2)

(2.4)

and the inverse Laplace transform of LT in inert form is

>

%invlaplace((1/2)*2^(1/2)*Pi^(1/2)*exp(-(1/2)/s)/s^(3/2), s, t)

(2.5)

At we have

>

%invlaplace((1/2)*2^(1/2)*Pi^(1/2)*exp(-(1/2)/s)/s^(3/2), s, 1)

(2.6)

>

(2.7)

This result is consistent with the one we get if we first compute the exact form of the inverse Laplace transform at t = 1:

>

%invlaplace((1/2)*2^(1/2)*Pi^(1/2)*exp(-(1/2)/s)/s^(3/2), s, 1) = sin(2^(1/2))

(2.8)

>

(2.9)

In addition to the standard use of evalf to numerically evaluate inverse Laplace transforms, one can invoke each of the three different methods implemented using the MathematicalFunctions:-Evalf command

>

(2.10)

>

(2.11)

>

(2.12)

>

(2.13)

Regarding the method we use by default: from a numerical experiment with varied problems we have concluded that our implementation of the Euler (sums) method is faster and more accurate than the other two.

Two Hankel transform definitions

In previous Maple releases, the definition of the Hankel transform was given by

hankel(f(t), t, s, nu) = Int(f(t)*sqrt(s*t)*BesselJ(nu, s*t), t = 0 .. infinity)

where is the function. This definition, sometimes called alternative definition of the Hankel transform, has the inconvenience of the square root in the integrand, complicating the form of the hankel transform for the Laplacian in cylindrical coordinates. On the other hand, the definition more frequently used in the literature is

hankel(f(t), t, s, nu) = Int(f(t)*t*BesselJ(nu, s*t), t = 0 .. infinity)

With it, the Hankel transform of diff(u(r, t), r, r)+(diff(u(r, t), r))/r+diff(u(r, t), t, t) is given by the simple ODE form d^2*`&Hopf;`(k, t)/dt^2-k^2*`&Hopf;`(k, t) . Not just this transform but several other ones acquire a simpler form with the definition that does not have a square root in the integrand.

So the idea is to align Maple with this simpler definition, while keeping the previous definition as an alternative. Hence, by default, when you load the inttrans package, the new definition in use for the Hankel transform is

>

hankel(f(t), t, s, nu) = Int(f(t)*t*BesselJ(nu, s*t), t = 0 .. infinity)

(3.1)

You can change this default so that Maple works with the alternative definition as in previous releases. For that purpose, use the new inttrans:-setup command (which you can also use to query about the definition in use at any moment):

>

(3.2)

This change in definition is automatically taken into account by other parts of the Maple library using the Hankel transform. For example, the differentiation rule with the new definition is

>

(3.3)

This differentiation rule resembles (is connected to) the differentiation rule for BesselJ, and this is another advantage of the new definition.

>

%diff(BesselJ(nu, z), z) = -BesselJ(nu+1, z)+nu*BesselJ(nu, z)/z

(3.4)

Furthermore, several transforms have acquired a simpler form, as for example:

>

%hankel(exp(I*a*r)/r, r, k, 0) = 1/(-a^2+k^2)^(1/2)

(3.5)

Let's compare: make the definition be as in previous releases.

>

(3.6)

>

hankel(f(t), t, s, nu) = Int(f(t)*(s*t)^(1/2)*BesselJ(nu, s*t), t = 0 .. infinity)

(3.7)

The differentiation rule with the previous (alternative) definition was not as simple:

>

(3.8)

And the transform (3.5) was also not so simple:

>

%hankel(exp(I*a*r)/r, r, k, 0) = (I*a*hypergeom([3/4, 3/4], [3/2], a^2/k^2)*GAMMA(3/4)^4+Pi^2*k*hypergeom([1/4, 1/4], [1/2], a^2/k^2))/(k*Pi*GAMMA(3/4)^2)

(3.9)

Reset to the new default value of the definition.

>

(3.10)

>

hankel(f(t), t, s, nu) = Int(f(t)*t*BesselJ(nu, s*t), t = 0 .. infinity)

(3.11)

More integral transform results

The revision of the integral transforms includes also filling gaps: previous transforms that were not being computed are now computed. Still with the Hankel transform, consider the operators

>

`D/t` := proc (u) options operator, arrow; (diff(u, t))/t end proc

Being able to transform these operators into algebraic expressions or differential equations of lower order is key for solving PDE problems with Boundary Conditions.

Consider, for instance, this ODE

>

(4.1)

>

((diff(diff(diff(u(t), t), t), t))*t^3-12*(diff(diff(u(t), t), t))*t^2+57*(diff(u(t), t))*t-105*u(t))/t^3

(4.2)

Its Hankel transform is a simple algebraic expression

>

(4.3)

An example with formula_plus

>

((diff(diff(diff(diff(u(t), t), t), t), t))*t^4+38*(diff(diff(diff(u(t), t), t), t))*t^3+477*(diff(diff(u(t), t), t))*t^2+2295*(diff(u(t), t))*t+3465*u(t))/t^4

(4.4)

>

(4.5)

In the case of hankel , not just differential operators but also several new transforms are now computable

>

piecewise(nu = 0, Dirac(k)/k, nu/k^2)

(4.6)

>

piecewise(And(nu = 0, m = 0), Dirac(k)/k, 2^(m+1)*k^(-m-2)*GAMMA(1+(1/2)*m+(1/2)*nu)/GAMMA((1/2)*nu-(1/2)*m))

(4.7)

>

Mellin and Hankel transform solutions for Partial Differential Equations with Boundary Conditions

In previous Maple releases, the Fourier and Laplace transforms were used to compute exact solutions to PDE problems with boundary conditions. Now, Mellin and Hankel transforms are also used for that same purpose.

Example:

>

pde := x^2*(diff(u(x, y), x, x))+x*(diff(u(x, y), x))+diff(u(x, y), y, y) = 0

iv := u(x, 0) = 0, u(x, 1) = piecewise(0 <= x and x < 1, 1, 1 < x, 0)

>

(5.1)

As usual, you can let pdsolve choose the solving method, or indicate the method yourself:

>

pde := diff(u(r, t), r, r)+(diff(u(r, t), r))/r+diff(u(r, t), t, t) = -Q__0*q(r)
iv := u(r, 0) = 0

>

u(r, t) = Q__0*(-hankel(exp(-s*t)*hankel(q(r), r, s, 0)/s^2, s, r, 0)+hankel(hankel(q(r), r, s, 0)/s^2, s, r, 0))

(5.2)

It is sometimes preferable to see these solutions in terms of more familiar integrals. For that purpose, use

>

u(r, t) = Q__0*(-(Int(exp(-s*t)*(Int(q(r)*r*BesselJ(0, r*s), r = 0 .. infinity))*BesselJ(0, r*s)/s, s = 0 .. infinity))+Int((Int(q(r)*r*BesselJ(0, r*s), r = 0 .. infinity))*BesselJ(0, r*s)/s, s = 0 .. infinity))

(5.3)

An example where the hankel transform is computable to the end:

>

pde := c^2*(diff(u(r, t), r, r)+(diff(u(r, t), r))/r) = diff(u(r, t), t, t)
iv := u(r, 0) = A*a/(a^2+r^2)^(1/2), (D[2](u))(r, 0) = 0

>

u(r, t) = (1/2)*A*a*((-c^2*t^2+(2*I)*a*c*t+a^2+r^2)^(1/2)+(-c^2*t^2-(2*I)*a*c*t+a^2+r^2)^(1/2))/((-c^2*t^2-(2*I)*a*c*t+a^2+r^2)^(1/2)*(-c^2*t^2+(2*I)*a*c*t+a^2+r^2)^(1/2))

(5.4)

>

Download Integral_Transforms_(revamped).mw

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

Graph Theory and Pokémon Application

by: vdbustos 40 Product: Maple 2019

September 19 2019

8 0

Hi there!

One of my favorite videogames is pokémon as you can probably guess from the title. As a player I always wanted to optimize my chances of obtaining the rarest and best pokémon in the game. I have been working on an application that aims to use graph theory to analyze the game Pokémon Blue. The application explores the following questions:

Which is the rarest pokémon in the game?
Where can I find an specific pokémon and with what probabilities?
What is the place with most different species of wild pokémon?

I also included algorithms for the following: Given a certain desired team

Find the minimum amount of places to visit to catch them and return the list of the places the player will need to visit.
What are the routes with best probabilities to catch each pokémon from my desired team?

Check out my application at: https://www.maplesoft.com/applications/view.aspx?SID=154565.

The following are some of the results obtained in the app:

What is the most common pokémon?

I did not only considered the amount of places a pokémon can appear in but also the probabilities of it appearing in each place.

What are the connections between pokémon and places?

In my graph, I connected a pokémon and a place if such pokémon could be caught in that place. The following is an example for the pokémon Pidgey. The weights of the edges are the probabilities of finding Pidgey in each route.

Viceversa, I did the same for how a route is connected to the pokémons in it:

Map of the Game
I also generated a colour coded version for the map of the game: where blue means that the place is a water route, brown means it's a cave and green means it's a tall-grass route.
It's amazing what Maple's graph theory toolbox can do.

Arabic Cipher

by: Fereydoon_Shekofte 223 Product: Maple 18

September 15 2019

3 11

This is a simple encryption method to hide text messages

Mentioned in Arabic manuscrips with more than hundreds years old ...

PRINCIPLE :

Just the place of letters in the sentence rearranged as described below :

For example "ABCDE" we pick up the First letter "A" from the left and write it as the last letter in the Right "......A"

but this time we pick up the letter "E" as the last letter from Right and place it at the Left Side of the previous one ".....EA"

and this cycle continue until for rest letters ... "CDBEA" .

by this way the text become hard to discover !

It is Amazing that for decoding this message you should repeat the same rearrangment algorithm several times until the readable text appears as the first "ABCDE"

EXample :

"AlbertEinstein"

"iEntsrteebilnA"

"eterbsitlnnEAi"

"tilsnbnrEeAtie"

"rnEbenAstliiet"

"sAtnleibiEentr"

"biieElennttArs"

"nenltEteAirisb"

"etAEitrlinsebn"

"lritnisEeAbtne"

"EseiAnbttinrel"

"tbtniAnireeslE"

"inrAeienstlbEt"

"nesitelAbrEnti"

"AlbertEinstein"

the same text appeared after 14 step cycle

Arabic Cipher

>

ArabicCipher := proc (x) options operator, arrow; StringTools[Permute](x, [seq(1+iquo(StringTools[Length](x), 2)+((1/2)*i+(1/2)*irem(i, 2))*(-1)^(i+irem(StringTools[Length](x), 2)), i = 0 .. StringTools[Length](x)-1)]) end proc

proc (x) options operator, arrow; StringTools[Permute](x, [seq(1+iquo(StringTools[Length](x), 2)+((1/2)*i+(1/2)*irem(i, 2))*(-1)^(i+irem(StringTools[Length](x), 2)), i = 0 .. StringTools[Length](x)-1)]) end proc

(1.1)

>

"iEntsrteebilnA", "eterbsitlnnEAi", "tilsnbnrEeAtie", "rnEbenAstliiet", "sAtnleibiEentr", "biieElennttArs", "nenltEteAirisb", "etAEitrlinsebn", "lritnisEeAbtne", "EseiAnbttinrel", "tbtniAnireeslE", "inrAeienstlbEt", "nesitelAbrEnti", "AlbertEinstein"

(1.2)

>

"nSohoedkyoefrteeF", "yokedferotheoeSFn", "otrheefodeeSkFony", "deoefSekeFhorntyo", "eFkheoSrfnetoyeod", "fnreStooeyhekoFde", "eyohoetkSoeFrdnef", "SoketFerodhnoeyfe", "odrhenFoteeykfoeS", "teoeFynkefhoredSo", "efkhnoyrFeedoSeot", "FereydoonShekofte"

(1.3)

>

Download Arabic_Cipher.mw

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post