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
  • A lot of scientific software propose packages enabling drawing figures in XKCD style/
    Up to now I thought this was restricted to open products (R, Python, ...) but I recently discovered Matlab and even Mathematica were doing same.

    Layton S (2012). “XKCDIFY! Adding flair to boring Matlab Axes one plot at a time.” Last accessed on December 08, 2014, URL

    Woods S (2012). “xkcd-style graphs.” Last accessed on December 08, 2014, URL 11355#11355.


    So why not Maple?

    As a regular user of R, I could have visualize the body of the corresponding procedures to see how these drawings were made and just translate theminto Maple.
    But copying for the sake of copying is not of much interest.
    So I started to develop some primitives for "XKCD-drawing" lines, polygons, circles and even histograms.
    My goal is not to write an XKCD package (I don't have the skills for that) but just to arouse the interest of (maybe) a few people here who could continue this preliminary work

    A main problem is the one of the XKCD fonts: no question to redefine them in Maple and I guess using them in a commercial code is not legal (?). So no XKCD font in this first work, nor even the funny guy who appears recurently on the drawings (but it could be easily constructed in Maple).

    In a recent post (Plot styling - experimenting with Maple's plotting...) Samir Khan proposed a few styles made of several plotting options,  some of which he named "Excel style" or "Oscilloscope style"... maybe a future "XKCD xtyle" in Maple ?

    This work has been done with Maple 2015 and reuses an old version of a 1D-Kriging procedure 





    The principle is always the same:
        1/   Let L a straight line which is either defined by its two ending points (xkvd_hline) or taken as the default [0, 0], [1, 0] line.
              For xkvd_hline the given line L is firstly rotate to be aligned with the horizontal axis.

        2/   Let P1, ..., PN N points on L. Each Pn writes [xn, yn]

        3/   A random perturbation rn is added yo the values y1, ..., yN

        4/   A stationnary random process RP, with gaussian correlation function is used to build a smooth curve passing through the points
              (x1, y1+r1), ..., (xN, yN+rN) (procedure KG where "KG" stands for "Kriging")

        5/   The result is drawn or mapped to some predefined shape :

        6/   A procedure xkcd_func is also provided to draw functions defined by an explicit relation.

    KG := proc(X, Y, psi, sigma)
      local NX, DX, K, mu, k, y:
      NX := numelems(X);
      DX := < seq(Vector[row](NX, X), k=1..NX) >^+:
      K  := (sigma, psi) -> evalf( sigma^2 *~ exp~(-((DX - DX^+) /~ psi)^~2) ):
      mu := add(Y) / NX;
      k  := (x, sigma, psi) -> evalf( convert(sigma^2 *~ exp~(-((x -~ X ) /~ psi)^~2), Vector[row]) ):
      y  := mu + k(x, sigma, psi) . (K(sigma, psi))^(-1) . convert((Y -~ mu), Vector[column]):
      return y
    end proc:

    xkcd_hline := proc(p1::list, p2::list, a::nonnegative, lc::positive, col)
      # p1 : first ending point
      # p2 : second ending point
      # a  : amplitude of the random perturbations
      # lc : correlation length
      # col: color
      local roll, NX, LX, X, Z:
      roll := rand(-1.0 .. 1.0):
      NX   := 10:
      LX   := p2[1]-p1[1]:
      X    := [seq(p1[1]..p2[1], LX/(NX-1))]:
      Z    := [p1[2], seq(p1[2]+a*roll(), k=1..NX-1)]:
      return plot(KG(X, Z, lc*LX, 1), x=min(X)..max(X), color=col, scaling=constrained):
    end proc:

    xkcd_line := proc(L::list, a::nonnegative, lc::positive, col, {lsty::integer:=1})
      # L  : list which contains the two ending point
      # a  : amplitude of the random perturbations
      # lc : correlation length
      # col: color
      local T, roll, NX, DX, DY, LX, A, m, M, X, Z, P:
      T    := (a, x0, y0, l) ->
                   (x,y) -> [ x0 + l * (x*cos(a)-y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ]
      roll := rand(-1.0 .. 1.0):
      NX   := 5:
      DX   := L[2][1]-L[1][1]:
      DY   := L[2][2]-L[1][2]:
      LX := sqrt(DX^2+DY^2):
      if DX <> 0 then
         A := arcsin(DY/LX):
         A:= Pi/2:
      end if:
      X := [seq(0..1, 1/(NX-1))]:
      Z := [ seq(a*roll(), k=1..NX)]:
      P := plot(KG(X, Z, lc, 1), x=0..1, color=col, scaling=constrained, linestyle=lsty):
      return T(A, op(L[1]), LX)(P)
    end proc:

    xkcd_func := proc(f, r::list, NX::posint, a::positive, lc::positive, col)
      # f  : function to draw
      # r  : plot range
      # NX : number of equidistant "nodes" in the range r (boundaries included)
      # a  : amplitude of the random perturbations
      # lc : correlation length
      # col: color
      local roll, F, LX, Pf, Xf, Zf:
      roll := rand(-1.0 .. 1.0):
      F    := unapply(f, indets(f, name)[1]);
      LX   := r[2]-r[1]:
      Pf   := [seq(r[1]..r[2], LX/(NX-1))]:
      Xf   := Pf +~ [seq(a*roll(), k=1..numelems(Pf))]:
      Zf   := F~(Pf) +~ [seq(a*roll(), k=1..numelems(Pf))]:
      return plot(KG(Xf, Zf, lc*LX, 1), x=min(Xf)..max(Xf), color=col):
    end proc:

    xkcd_hist := proc(H, ah, av, ax, ay, lch, lcv, lcx, lcy, colh, colxy)
      # H   : Histogram
      # ah  : amplitude of the random perturbations on the horizontal boundaries of the bins
      # av  : amplitude of the random perturbations on the vertical boundaries of the bins
      # ax  : amplitude of the random perturbations on the horizontal axis
      # ay  : amplitude of the random perturbations on the vertical axis
      # lch : correlation length on the horizontal boundaries of the bins
      # lcv : correlation length on the vertical boundaries of the bins
      # lcx : correlation length on the horizontal axis
      # lcy : correlation length on the vertical axis
      # colh: color of the histogram
      # col : color of the axes
      local data, horiz, verti, horizontal_lines, vertical_lines, po, rpo, p1, p2:
      data  := op(1..-2, op(1, H)):
      verti := sort( [seq(data[n][3..4][], n=1..numelems([data]))] , key=(x->x[1]) );
      verti := verti[1],
                    n -> if verti[n][2] > verti[n+1][2] then
                          end if,
      horiz := seq(data[n][[4, 3]], n=1..numelems([data])):

      horizontal_lines := NULL:
      for po in horiz do
        horizontal_lines := horizontal_lines, xkcd_hline(po[1], po[2], ah, lch, colh):
      end do:

      vertical_lines := NULL:
      for po in [verti] do
        rpo := po[[2, 1]]:
        vertical_lines := vertical_lines, xkcd_hline([0, rpo[2]], rpo, av, lcv, colh):
      end do:

      p1 := [2*verti[1][1]-verti[2][1], 0]:
      p2 := [2*verti[-1][1]-verti[-2][1], 0]:

          xkcd_hline(p1, p2, ax, lcx, colxy),
          T(xkcd_hline([0, 0], [1.2*max(op~(2, [verti])), 0], ay, lcy, colxy)),
    end proc:

    xkcd_polyline := proc(L::list, a::nonnegative, lc::positive, col)
      # xkcd_polyline reduces to xkcd_line whebn L has 2 elements
      # L  : List of points
      # a  : amplitude of the random perturbations
      # lc : correlation length
      # col: color
      local T, roll, NX, n, DX, DY, LX, A, m, M, X, Z, P:
      T    := (a, x0, y0, l) ->
                   (x,y) -> [ x0 + l * (x*cos(a)-y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ]
      roll := rand(-1.0 .. 1.0):
      NX   := 5:
      for n from 1 to numelems(L)-1 do
        DX   := L[n+1][1]-L[n][1]:
        DY   := L[n+1][2]-L[n][2]:
        LX := sqrt(DX^2+DY^2):
        if DX <> 0 then
          A := evalf(arcsin(abs(DY)/LX)):
          if DX >= 0 and DY <= 0 then A := -A end if:
          if DX <= 0 and DY >  0 then A := Pi-A end if:
          if DX <= 0 and DY <= 0 then A := Pi+A end if:
          A:= Pi/2:
          if DY < 0 then A := 3*Pi/2 end if:
        end if:
        if n=1 then
          X := [seq(0..1, 1/(NX-1))]:
          Z := [seq(a*roll(), k=1..NX)]:
          X := [0    , seq(1/(NX-1)..1, 1/(NX-1))]:
          Z := [Z[NX], seq(a*roll(), k=1..NX-1)]:
        end if:
        P    := plot(KG(X, Z, lc, 1), x=0..1, color=col, scaling=constrained):
        P||n := T(A, op(L[n]), LX)(P):
      end do;
      return seq(P||n, n=1..numelems(L)-1)
    end proc:

    xkcd_circle := proc(a::nonnegative, lc::positive, r::positive, cent::list, col)
      # a   : amplitude of the random perturbations
      # lc  : correlation length
      # r   : redius of the circle
      # cent: center of the circle
      # col : color
      local roll, NX, LX, X, Z, xkg, A:
      roll := rand(-1.0 .. 1.0):
      NX   := 10:
      X    := [seq(0..1, 1/(NX-1))]:
      Z    := [0, seq(a*roll(), k=1..NX-1)]:
      xkg  := KG(X, Z, lc, 1):
      A    := Pi*roll():
      return plot([cent[1]+r*(1+xkg)*cos(2*Pi*x+A), cent[2]+r*(1+xkg)*sin(2*Pi*x+A), x=0..1], color=col)
    end proc:

    T := plottools:-transform((x,y) -> [y, x]):

    # Axes plot

    x_axis := xkcd_hline([0, 0], [10, 0], 0.03, 0.5, black):
    y_axis := xkcd_hline([0, 0], [10, 0], 0.03, 0.5, black):


    # A simple function

    f := 1+10*(x/5-1)^2:
    F := xkcd_func(f, [0.5, 9.5], 6, 0.3, 0.4, red):



    # An histogram

    S := Sample(Normal(0,1),100):
    H := Histogram(S, maxbins=6):
    xkcd_hist(H,   0, 0.02, 0.001, 0.01,   1, 0.1, 0.01, 1,   red, black)


    # Axes plus grid with two red straight lines

    r := rand(-0.1 .. 0.1):

    x_axis := xkcd_line([[-2, 0], [10, 0]], 0.01, 0.2, black):
    y_axis := xkcd_line([[0, -2], [0, 10]], 0.01, 0.2, black):
    d1     := xkcd_line([[-1, 1], [9, 9]] , 0.01, 0.2, red):
    d2     := xkcd_line([[-1, 9], [9, -1]], 0.01, 0.2, red):
      x_axis, y_axis,
      seq( xkcd_line([[-2+0.3*r(), u+0.3*r()], [10+0.3*r(), u+0.3*r()]], 0.005, 0.5, gray), u in [seq(-1..9, 2)]),
      seq( xkcd_line([[u+0.3*r(), -2+0.3*r()], [u+0.3*r(), 10+0.3*r()]], 0.005, 0.5, gray), u in [seq(-1..9, 2)]),
      d1, d2,


    # Axes and a couple of polygonal lines

    d1 := xkcd_polyline([[0, 0], [1, 3], [3, 5], [7, 1], [9, 7]], 0.01, 1, red):
    d2 := xkcd_polyline([[0, 9], [2, 8], [5, 2], [8, 3], [10, -1]], 0.01, 1, blue):

      x_axis, y_axis,
      d1, d2,


    # A few polygonal shapes

      xkcd_polyline([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]], 0.01, 1, red),
      xkcd_polyline([[1/3, 1/3], [2/3, 1/3], [2/3, 4/3], [-1, 4/3], [1/3, 1/3]], 0.01, 1, blue),
      xkcd_polyline([[-1/3, -1/3], [4/3, 1/2], [1/2, 1/2], [1/2,-1], [-1/3, -1/3]], 0.01, 1, green),


    # A few circles

    cols  := [red, green, blue, gold, black]:                                # colors
    cents := convert( Statistics:-Sample(Uniform(-1, 3), [5, 2]), listlist): # centers
    radii := Statistics:-Sample(Uniform(1/2, 2), 5):                         # radii
    lcs   := Statistics:-Sample(Uniform(0.2, 0.7), 5):                       # correlations lengths

        xkcd_circle(0.02, lcs[n], radii[n], cents[n], cols[n]),


    # A 3D drawing

    x_axis := xkcd_line([[0, 0], [5, 0]], 0.01, 0.2, black):
    y_axis := xkcd_line([[0, 0], [4, 2]], 0.01, 0.2, black):
    z_axis := xkcd_line([[0, 0], [0, 5]], 0.01, 0.2, black):

    f1 := 4*cos(x/6)-1:
    F1 := xkcd_func(f1, [0.5, 5], 6, 0.001, 0.8, red):
    F2 := xkcd_line([[0.5, eval(f1, x=0.5)], [3, 4]], 0.01, 0.2, red):
    f3 := 4*cos((x-2)/6):
    F3 := xkcd_func(f3, [3, 7], 6, 0.001, 0.8, red):
    F4 := xkcd_line([[5, eval(f1, x=5)], [7, eval(f3, x=7)]], 0.01, 0.2, red):

    dx := xkcd_line([[2, 1], [4, 1]], 0.01, 0.2, gray, lsty=3):
    dy := xkcd_line([[2, 0], [4, 1]], 0.01, 0.2, gray, lsty=3):
    dz := xkcd_line([[4, 1], [4, 3]], 0.01, 0.2, gray, lsty=3):

    po := xkcd_circle(0.02, 0.3, 0.1, [4, 3], blue):

    # Numerical value come from "probe info + copy/paste"

    nvect   := xkcd_polyline([[4, 3], [4.57, 4.26], [4.35, 4.1], [4.57, 4.26], [4.58, 4.02]], 0.01, 1, blue):
    tg1vect := xkcd_polyline([[4, 3], [4.78, 2.59], [4.49, 2.87], [4.78, 2.59], [4.46, 2.57]], 0.01, 1, blue):
    tg2vect := xkcd_polyline([[4, 3], [4.79, 3.35], [4.70, 3.13], [4.79, 3.35], [4.46, 3.35]], 0.01, 1, blue):
    rec1    := xkcd_polyline([[4.118, 3.286], [4.365, 3.396], [4.257, 3.108]], 0.01, 1, blue):
    rec2    := xkcd_polyline([[4.257, 3.108], [4.476, 2.985], [4.259, 2.876]], 0.01, 1, blue):

      x_axis, y_axis, z_axis,
      F1, F2, F3, F4,
      dx, dy, dz,
      nvect, tg1vect, tg2vect, rec1, rec2,


    # Arrow

    d1 := xkcd_polyline([[0, 0], [1, 0], [0.9, 0.05], [1, 0], [0.9, -0.05]], 0.01, 1, red):

    T := (a, x0, y0, l) ->
                   (x,y) -> [ x0 + l * (x*cos(a)-y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ]

      seq( T(2*Pi*n/10, 0.5, 0, 1/2)(
                      [[0, 0], [1, 0], [0.9, 0.05], [1, 0], [0.9, -0.05]],
                      ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12])






    An attempt to find the equation of an ellipse inscribed in a given triangle. 
    The program works on the basis of the ELS procedure.  After the procedure works, the  solutions are filtered.
    ELS procedure solves the system of equations f1, f2, f3, f4, f5 for the coefficients of the second-order curve.
    The equation f1 corresponds to the condition that the side of the triangle intersects t a curve of the second order at one point.
    The equation f2 corresponds to the condition that the point x1,x2  belongs to a curve of the second order.
    Equation f3 corresponds to the condition that the side of the triangle is tangent to the second order curve at the point x1,x2.
    The equation f4 is similar to the equation f2, and the equation f5 is similar to the equation f3.
    For example

    @ianmccr posted here about making help for a package using makehelp. Here I show how to do this with the HelpTools package.

    The attached worksheet shows how to create the help database for the Orbitals package available at the Applications Center or in the Maple Cloud. The help pages were created as worksheets - start using an existing help page as a model - use View/Open Page As Worksheet and then save from there. The topics and other information are entered by adding Attributes under File/Document Properties - for example for the realY help page these are:

    Active=false means a regular help page; Active=true means an example worksheet.

    There may be several aliases, for example the cartesion help page also describes the fullcartesian command and so Alias is: Orbitals[fullcartesian],cartesian,fullcartesian and Keyword is: Orbitals,cartesian,fullcartesian

    Once a worksheet is created for each help page they are assembled into the help database with the attached file. More information is in the attached file

    MacDude posted a worksheet for creating help files using the makehelp command that I found very useful.  However, his worksheet created a table of contents that organized the command help pages under a single folder.  I wanted to create a help file table of contents with the commands organized into sub-folders under the main application folder. ie. Package Name,Folder Name, command. The help page for makehelp is not very informative; in particular the example that shows (purportedly) how to override the existing help pages with your own help file is very misleading.  Also, the description of the parameters to the browser option of the makehelp command is too vague. I needed an error message to tell me that the browser option expects parameters in the form List(Name,String).  In the end, I was able to get a folder/sub-folder structure using a structure [`Folder Name`,`Subfolder Name`,"Command"].  Sub-folders are sorted alphabetically. If anyone has a need, the attached worksheet shows how I created the table of contents structure.

    I like tweaking plots to get the look and feel I want, and luckily Maple has many plotting options that I often play with. Here, I visualize the same data several times, but each time with different styling.

    First, some data.

    data_1 := [[0,0],[1,2],[2,1.3],[3,6]]:
    data_2 := [[0.5,3],[1,1],[2,5],[3,2]]:
    data_3 := [[-0.5,3],[1.3,1],[2.5,5],[4.5,2]]:

    This is the default look.

    plot([data_1, data_2, data_3])

    I think the darker background on this plot makes it easier to look at.

    gray_grid :=
     background      = "LightGrey"
    ,color           = [ ColorTools:-Color("RGB",[150/255, 40 /255, 27 /255])
                        ,ColorTools:-Color("RGB",[0  /255, 0  /255, 0  /255])
                        ,ColorTools:-Color("RGB",[68 /255, 108/255, 179/255]) ]
    ,axes            = frame
    ,axis[2]         = [color = black, gridlines = [10, thickness = 1, color = ColorTools:-Color("RGB", [1, 1, 1])]]
    ,axis[1]         = [color = black, gridlines = [10, thickness = 1, color = ColorTools:-Color("RGB", [1, 1, 1])]]
    ,axesfont        = [Arial]
    ,labelfont       = [Arial]
    ,size            = [400*1.78, 400]
    ,labeldirections = [horizontal, vertical]
    ,filled          = false
    ,transparency    = 0
    ,thickness       = 5
    ,style           = line:
    plot([data_1, data_2, data_3], gray_grid);

    I call the next style Excel, for obvious reasons.

    excel :=
     background      = white
    ,color           = [ ColorTools:-Color("RGB",[79/255,  129/255, 189/255])
                        ,ColorTools:-Color("RGB",[192/255, 80/255,   77/255])
                        ,ColorTools:-Color("RGB",[155/255, 187/255,  89/255])]
    ,axes            = frame
    ,axis[2]         = [gridlines = [10, thickness = 0, color = ColorTools:-Color("RGB",[134/255,134/255,134/255])]]
    ,font            = [Calibri]
    ,labelfont       = [Calibri]
    ,size            = [400*1.78, 400]
    ,labeldirections = [horizontal, vertical]
    ,transparency    = 0
    ,thickness       = 3
    ,style           = point
    ,symbol          = [soliddiamond, solidbox, solidcircle]
    ,symbolsize      = 15:
    plot([data_1, data_2, data_3], excel)

    This style makes the plot look a bit like the oscilloscope I have in my garage.

    dark_gridlines :=
     background      = ColorTools:-Color("RGB",[0,0,0])
    ,color           = white
    ,axes            = frame
    ,linestyle       = [solid, dash, dashdot]
    ,axis            = [gridlines = [10, linestyle = dot, color = ColorTools:-Color("RGB",[0.5, 0.5, 0.5])]]
    ,font            = [Arial]
    ,size            = [400*1.78, 400]:
    plot([data_1, data_2, data_3], dark_gridlines);

    The colors in the next style remind me of an Autumn morning.

    autumnal :=
     background      =  ColorTools:-Color("RGB",[236/255, 240/255, 241/255])
    ,color           = [  ColorTools:-Color("RGB",[144/255, 54/255, 24/255])
                         ,ColorTools:-Color("RGB",[105/255, 108/255, 51/255])
                         ,ColorTools:-Color("RGB",[131/255, 112/255, 82/255]) ]
    ,axes            = frame
    ,font            = [Arial]
    ,size            = [400*1.78, 400]
    ,filled          = true
    ,axis[2]         = [gridlines = [10, thickness = 1, color = white]]
    ,axis[1]         = [gridlines = [10, thickness = 1, color = white]]
    ,symbol          = solidcircle
    ,style           = point
    ,transparency    = [0.6, 0.4, 0.2]:
    plot([data_1, data_2, data_3], autumnal);

    In honor of a friend and ex-colleague, I call this style "The Swedish".

    swedish :=
     background      = ColorTools:-Color("RGB", [0/255, 107/255, 168/255])
    ,color           = [ ColorTools:-Color("RGB",[169/255, 158/255, 112/255])
                        ,ColorTools:-Color("RGB",[126/255,  24/255,   9/255])
                        ,ColorTools:-Color("RGB",[254/255, 205/255,   0/255])]
    ,axes            = frame
    ,axis            = [gridlines = [10, color = ColorTools:-Color("RGB",[134/255,134/255,134/255])]]
    ,font            = [Arial]
    ,size            = [400*1.78, 400]
    ,labeldirections = [horizontal, vertical]
    ,filled          = false
    ,thickness       = 10:
    plot([data_1, data_2, data_3], swedish);

    This looks like a plot from a journal article.

    experimental_data_mono :=
    background       = white
    ,color           = black
    ,axes            = box
    ,axis            = [gridlines = [linestyle = dot, color = ColorTools:-Color("RGB",[0.5, 0.5, 0.5])]]
    ,font            = [Arial, 11]
    ,legendstyle     = [font = [Arial, 11]]
    ,size            = [400, 400]
    ,labeldirections = [horizontal, vertical]
    ,style           = point
    ,symbol          = [solidcircle, solidbox, soliddiamond]
    ,symbolsize      = [15,15,20]:
    plot([data_1, data_2, data_3], experimental_data_mono, legend = ["Annihilation", "Authority", "Acceptance"]);

    If you're willing to tinker a little bit, you can add some real character and personality to your visualizations. Try it!

    I'd also be very interested to learn what you find attractive in a plot - please do let me know.

    I'm really fed up with installing the Physics package!

    My install directory of Maple is C:\Maths\Maple

    I don't want to do it manually each time a new Physics Package is installed.

    "C:\Maths\Maple\lib", "C:\Users\jm\maple\toolbox\GRTensorIII\lib"

    The "Physics Updates" version "717" is installed but is not

       active. The active version of Physics is within the library

       C:\Maths\Maple\lib\maple.mla, created 2020, March 5, 2:36 hours

    Please software developpers at Maplesoft do something...Before rel 713 everything was OK.

    Thank you and kind regards to each and every one.



    I would like to share this work I've done. 
    No big math here, just a demonstrator of Maple's capabilities in 3D visualization.

    All the plots in the file have been discarded to reduce the size of this post. Here is a screen capture to give you an idea of what is inside the file.



    In a recent post  (Monte Carlo Integration) Radaar shared its work about the numerical integration, with the Monte Carlo method, of a function defined in polar coordinates.
    Radaar used a raw strategy based on a sampling in cartesian coordinates plus an ad hoc transformation.
    Radaar obtained reasonably good results, but I posted a comment to show how Monte Carlo summation in polar coordinates can be done in a much simpler way. Behind this is the choice of a "good" sampling distribution which makes the integration problem as simple as Monte Carlo integration over a 2D rectangle with sides parallel to the co-ordinate axis.

    This comment I sent pushed me to share the present work on Monte Carlo integration over simple polygons ("simple" means that two sides do not intersect).
    Here again one can use raw Monte Carlo integration on the rectangle this polygon is inscribed in. But as in Radaar's post, a specific sampling distribution can be used that makes the summation method more elegant.

    This work relies on three main ingredients:

    1. The Dirichlet distribution, whose one form enables sampling the 2D simplex in a uniform way.
    2. The construction of a 1-to-1 mapping from this simplex into any non degenerated triangle (a mapping whose jacobian is a constant equal to the ratio of the areas of the two triangles).
    3. A tesselation into triangles of the polygon to integrate over.

    This work has been carried out in Maple 2015, which required the development of a module to do the tesselation. Maybe more recent Maple's versions contain internal procedures to do that.


    This is an animation of the spread of the COVID-19 over the U.S. in the first 150 days.  It was created in Maple 2020, making extensive use of DataFrames. 



    The animation of 150 Day history includes COVID-19 data published by the NY Times and geographic data assembled from other sources. Each cylinder represents a county or in two special cases New York City and Kansas City. The cross-sectional area of each cylinder is the area in square miles of the corresponding county. The height of each cylinder is on a logarithmic scale (in particular it is 100*log base 2 of the case count for the county. The argument of the logarithm function is the number of cases per county divided by the are in square miles-so an areal density.  Using a logarithmic scale facilitates showing super high density areas (e.g., NYC) along with lower density areas.  The heights are scaled by a prefactor of 100 for visualization.

    Hi. My name is Eugenio and I’m a Professor at the Departamento de Didáctica de las Ciencias Experimentales, Sociales y Matemáticas at the Facultad de Educación of the Universidad Complutense de Madrid (UCM) and a member of the Instituto de Matemática Interdisciplinar (IMI) of the UCM.

    I have a 14-year-old son. In the beginning of the pandemic, a confinement was ordered in Spain. It is not easy to make a kid understand that we shouldn't meet our friends and relatives for some time and that we should all stay at home in the first stage. So, I developed a simplified explanation of virus propagation for kids, firstly in Scratch and later in Maple, the latter using an implementation of turtle geometry that we developed long ago and has a much better graphic resolution (E. Roanes-Lozano and E. Roanes-Macías. An Implementation of “Turtle Graphics” in Maple V. MapleTech. Special Issue, 1994, 82-85). A video (in Spanish) of the Scratch version is available from the Instituto de Matemática Interdisciplinar (IMI) web page:


    Surely you are uncomfortable being locked up at home, so I will try to justify that, although we are all looking forward going out, it is good not to meet your friends and family with whom you do not live.

    I firstly need to mention a fractal is. A fractal is a geometric object whose structure is repeated at any scale. An example in nature is Romanesco broccoli, that you perhaps have eaten (you can search for images on the Internet). You can find a simple fractal in the following image (drawn with Maple):

    Notice that each branch is divided into two branches, always forming the same angle and decreasing in size in the same proportion.

    We can say that the tree in the previous image is of “depth 7” because there are 7 levels of branches.

    It is quite easy to create this kind of drawing with the so called “turtle geometry” (with a recursive procedure, that is, a procedure that calls itself). Perhaps you have used Scratch programming language at school (or Logo, if you are older), which graphics are based in turtle geometry.

    All drawings along these pages have been created with Maple. We can easily reform the code that generated the previous tree so that it has three, four, five,… branches at each level, instead of two.

    But let’s begin with a tale that explains in a much simplified way how the spread of a disease works.

    - o O o -

    Let's suppose that a cat returns sick to Catland suffering from a very contagious disease and he meets his friends and family, since he has missed them so much.

    We do not know very well how many cats each sick cat infects in average (before the order to STAY AT HOME arrives, as cats in Catland are very obedient and obey right away). Therefore, we’ll analyze different scenarios:

    1. Each sick cat infects two other cats.
    2. Each sick cat infects three other cats.
    3. Each sick cat infects five other cats


    1. Each Sick Cat Infects Two Cats

    In all the figures that follow, the cat initially sick is in the center of the image. The infected cats are represented by a red square.

    · Before everyone gets confined at home, it only takes time for that first sick cat to see his friends, but then confinement is ordered (depth 1)

    As you can see, with the cat meeting his friends and family, we already have 3 sick cats.

    · Before all cats confine themselves at home, the first cat meets his friends, and these in turn have time to meet their friends (depth 2)

    In this case, the number of sick cats is 7.

    · Before every cat is confined at home, there is time for the initially sick cat to meet his friends, for these to meet their friends, and for the latter (friends of the friends of the first sick cat) to meet their friends (depth 3).

    There are already 15 sick cats...

    · Depth 4: 31 sick cats.

    · Depth 5: 63 sick cats.

    Next we’ll see what would happen if each sick cat infected three cats, instead of two.


    2. Every Sick Cat Infects Three Cats

    · Now we speed up, as you’ve got the idea.

    The first sick cat has infected three friends or family before confining himself at home. So there are 4 infected cats.

    · If each of the recently infected cats in the previous figure have in turn contact with their friends and family, we move on to the following situation, with 13 sick cats:

    · And if each of these 13 infected cats lives a normal life, the disease spreads even more, and we already have 40!

    · At the next step we have 121 sick cats:

    · And, if they keep seeing friends and family, there will be 364 sick cats (the image reminds of what is called a Sierpinski triangle):


    4. Every Sick Cat Infects Five Cats

    · In this case already at depth 2 we already have 31 sick cats.


    5. Conclusion

    This is an example of exponential growth. And the higher the number of cats infected by each sick cat, the worse the situation is.

    Therefore, avoiding meeting friends and relatives that do not live with you is hard, but good for stopping the infection. So, it is hard, but I stay at home at the first stage too!

    We have just released an update to Maple, Maple 2020.1.

    Maple 2020.1 includes corrections and improvement to the mathematics engine, export to PDF, MATLAB connectivity, support for Ubuntu 20.04, and more.  We recommend that all Maple 2020 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 2020.1 download page, where you can also find more details.

    In particular, please note that this update includes a fix to the SMTLIB problem reported on MaplePrimes. Thanks for the feedback!

    Monte Carlo integration uses random sampling unlike classical techniques like the trapezoidal or Simpson's rule in evaluating the integration numerically.

    restart; ff := proc (rho, phi) return exp(rho*cos(phi))*rho end proc; aa := 0; bb := 1; cc := 0; dd := 2*Pi; alfa := 5; nrun := 15000; sum1 := 0; sum2 := 0; X := Statistics:-RandomVariable(Uniform(0, 1)); SX := Statistics:-Sample(X); for ii to nrun do u1 := SX(1)[1]; u2 := SX(1)[1]; xx1 := aa+(bb-aa)*u1; xx2 := cc+(dd-cc)*u2; xx3 := (bb-aa)*(1-u1); xx4 := (dd-cc)*(1-u2); sum1 := sum1+evalf(ff(xx1, xx2)); sum2 := sum2+evalf(ff(xx1, xx2))+evalf(ff(xx1, xx4))+evalf(ff(xx3, xx2))+evalf(ff(xx3, xx4)) end do; area1 := (bb-aa)*(dd-cc)*sum1/nrun; area2 := (bb-aa)*(dd-cc)*sum2/(4*nrun); area2



    evalf(Int(exp(rho*cos(phi))*rho, rho = 0 .. 1, phi = 0 .. 2*Pi))








    The purpose of this document is:

    a) to correct the physics that was used in the document "Minimal Road Radius for Highway Superelevation" recently submitted to the Maple Applications Center;

    b) to confirm the values found in the manual for the American Association of State Highway and Transportation Officials (AASHTO) that engineers use to design and build these banked curves are physically sound. 

    c) to highlight the pedagogical value inherent in the Maple language to distinguish between assignment ( := )  and equivalence (  =  );

    d) but most importantly, to demonstrate the pedagogical value Maple has in thinking about solving a problem involving a physical process. Given Maple's symbolic mathematics capabilities, one can implement a top-down approach to the physics and the mathematics, working from the general principle to the specific example. This allows one to avoid the types of errors that occur when translating the problem into a bottom up approach, from specific values of the example to the general principle, an approach that is required by most other computational systems.

    I hope that others are willing to continue to engage in discussions related to the pedagogical value of Maple beyond mathematics.

    I was asked to post this document to both here and the Maple Applications Center

    [Document edited for typos.]


    We’re excited to announce a new version of MapleSim! The MapleSim 2020 family of products lets you build and test models faster than ever, including faster simulations, powerful new tools for machine builders, and expanded modeling capabilities. Improvements include:

    • Faster results, with more efficient models, faster simulations, and more powerful design tools.
    • Powerful new features for machine builders, with new components, improved visualizations, and automation-focused connectivity tools that make it faster than ever to build and test digital twins.
    • Improved modeling capabilities, with an extensive collection of updates to components, libraries, and analysis tools.
    • More realistic machine visualizations with an expanded Kinematic Cam Generation App.
    • New product: MapleSim Insight, giving machine builders powerful, simulation-based debugging and 3-D visualization capabilities that connect directly to common automation platforms.
    • New add-on library:  MapleSim Ropes and Pulleys Library for the easy creation of winch and pulley systems as part of your machine development. 

    See What’s New in MapleSim 2020 for more information about these and other improvements!

    Maple's pdsolve() is quite capable of solving the PDE that describes the motion of a single-span Euler beam.  As far as I have been able to ascertain, there is no obvious way of applying pdsolve() to solve multi-span beams.  The worksheet attached to this post provides tools for solving multi-span Euler beams.  Shown below are a few demos.  The worksheet contains more demos.


    A module for solving the Euler beam with the method of lines



    The beamsolve proc solves a (possibly multi-span) Euler beam equation:``

    "rho ((&PartialD;)^2u)/((&PartialD;)^( )t^2)+ ((&PartialD;)^2)/((&PartialD;)^( )x^2)(EI ((&PartialD;)^(2)u)/((&PartialD;)^( )x^(2)))=f"

    subject to initial and boundary conditions.  The solution u = u(x, t) is the

    transverse deflection of the beam at point x at time t, subject to the load
    density (i.e., load per unit length) given by f = f(x, t). The coefficient rho 

    is the beam's mass density (mass per unit length), E is the Young's modulus of

    the beam's material, and I is the beam's cross-sectional moment of inertia

    about the neutral axis.  The figure below illustrates a 3-span beam (drawn in green)
    supported on four supports, and loaded by a variable density load (drawn in gray)
    which may vary in time.  The objective is to determine the deformed shape of the
    beam as a function of time.

    The number of spans, their lengths, and the nature of the supports may be specified

    as arguments to beamsolve.


    In this worksheet we assume that rho, E, I are constants for simplicity. Since only
    the product of the coefficients E and I enters the calculations, we lump the two

    together into a single variable which we indicate with the two-letter symbol EI.

    Commonly, EI is referred to as the beam's rigidity.


    The PDE needs to be supplied with boundary conditions, two at each end, each

    condition prescribing a value (possibly time-dependent) for one of u, u__x, u__xx, u__xxx 
    (that's 36 possible combinations!) where I have used subscripts to indicate

    derivatives.  Thus, for a single-span beam of length L, the following is an admissible

    set of boundary conditions:
    u(0, t) = 0, u__xx(0, t) = 0, u__xx(L, t) = 0, u__xxx(t) = sin*t.   (Oops, coorection, that last
    condition was meant to be uxxx(L,t) = sin t.)

    Additionally, the PDE needs to be supplied with initial conditions which express

    the initial displacement and the initial velocity:
    "u(x,0)=phi(x),   `u__t`(x,0)=psi(x)."


    The PDE is solved through the Method of Lines.  Thus, each span is subdivided into

    subintervals and the PDE's spatial derivatives are approximated through finite differences.

    The time derivatives, however, are not discretized.  This reduces the PDE into a set of

    ODEs which are solved with Maple's dsolve().  


    Calling sequence:

            beamsolve(L, n, options)



            L:  List of span lengths, in order from left to right, as in [L__1, L__2 .. (), `L__&nu;`].

            n The number of subintervals in the shortest span (for the finite difference approximation)




    It is assumed that the spans are laid back-to-back along the x axis, with the left end
    of the overall beam at x = 0.


    The interior supports, that is, those supports where any two spans meet, are assumed
    to be of the so-called simple type.  A simple support is immobile and it doesn't exert
    a bending moment on the beam.  Supports at the far left and far right of the beam can
    be of general type; see the BC_left and BC_right options below.


    If the beam consists of a single span, then the argument L may be entered as a number
    rather than as a list. That is, L__1 is equivalent to [L__1].



            All options are of the form option_name=value, and have built-in default values.

            Only options that are other than the defaults need be specified.


            rho: the beam's (constant) mass density per unit length (default: rho = 1)

            EI: the beam's (constant) rigidity (default: EI = 1)

            T: solve the PDE over the time interval 0 < t and t < T (default: T = 1)

            F: an expression in x and t that describes the applied force f(x, t)  (default: F = 0)
            IC: the list [u(x, 0), u__t(x, 0)]of the initial displacement and velocity,  as
                    expressions in x (default: IC = [0,0])

            BC_left: a list consisting of a pair of boundary conditions at the left end of
                    the overall (possibly multi-span beam.  These can be any two of
                    u = alpha(t), u_x = beta(t), u_xx = gamma(t), u_xxx = delta(t). The right-hand sides of these equations

                    can be any expression in t.  The left-hand sides should be entered literally as indicated.

                    If a right-hand side is omitted, it is taken to be zero.   (default: BC_left = [u, u_xx] which

                    corresponds to a simple support).

            BC_right: like BC_left, but for the right end of the overall beam (default: BC_right = "[u,u_xx])"


    The returned module:

            A call to beamsolve returns a module which presents three methods.  The methods are:


            plot (t, refine=val, options)

                    plots the solution u(x, t) at time t.  If the discretization in the x direction

                    is too coarse and the graph looks non-smooth, the refine option

                    (default: refine=1) may be specified to smooth out the graph by introducing

                    val number of intermediate points within each discretized subinterval.

                    All other options are assumed to be plot options and are passed to plots:-display.


            plot3d (m=val, options)

                    plots the surface u(x, t).  The optional m = val specification requests

                    a grid consisting of val subintervals in the time direction (default: "m=25)"

                    Note that this grid is for plotting purposes only; the solution is computed

                    as a continuous (not discrete) function of time. All other options are assumed

                    to be plot3d options and are passed to plots:-display.


            animate (frames=val, refine=val, options)

                    produces an animation of the beam's motion.  The frames option (default = 50)

                    specifies the number of animation frames.  The refine option is passed to plot
                    (see the description above. All other options are assumed to be plot options and
                    are passed to plots:-display.


            In specifying the boundary conditions, the following reminder can be helpful.  If the beam

            is considered to be horizontal, then u is the vertical displacement, `u__x ` is the slope,  EI*u__xx

            is the bending moment, and EI*u__xxx is the transverse shear force.


    A single-span simply-supported beam with initial velocity


    The function u(x, t) = sin(Pi*x)*sin(Pi^2*t) is an exact solution of a simply supported beam with

    "u(x,0)=0,   `u__t`(x,0)=Pi^(2)sin(Pi x)."  The solution is periodic in time with period 2/Pi.

    sol := beamsolve(1, 25, 'T'=2/Pi, 'IC'=[0, Pi^2*sin(Pi*x)]):

    The initial condition u(x, 0) = 0, u__t(x, 0) = 1  does not lead to a separable form, and

    therefore the motion is more complex.

    sol := beamsolve(1, 25, 'T'=2/Pi, 'IC'=[0, 1]):
    sol:-animate(frames=200, size=[600,250]);


    A single-span cantilever beam


    A cantilever beam with initial condition "u(x,0)=g(x),  `u__t`(x,0)=0," where g(x) is the
    first eigenmode of its free vibration (calculated in another spreadsheet).  The motion is
    periodic in time, with period "1.787018777."

    g := 0.5*cos(1.875104069*x) - 0.5*cosh(1.875104069*x) - 0.3670477570*sin(1.875104069*x) + 0.3670477570*sinh(1.875104069*x):
    sol := beamsolve(1, 25, 'T'=1.787018777, 'BC_left'=[u,u_x], 'BC_right'=[u_xx,u_xxx], 'IC'=[g, 0]):

    If the initial condition is not an eigenmode, then the solution is rather chaotic.

    sol := beamsolve(1, 25, 'T'=3.57, 'BC_left'=[u,u_x], 'BC_right'=[u_xx,u_xxx], 'IC'=[-x^2, 0]):
    sol:-animate(size=[600,250], frames=100);


    A single-span cantilever beam with a weight hanging from its free end


    sol := beamsolve(1, 25, 'T'=3.57, 'BC_left'=[u,u_x], 'BC_right'=[u_xx,u_xxx=1]):
    sol:-animate(size=[600,250], frames=100);


    A single-span cantilever beam with oscillating support


    sol := beamsolve(1, 25, 'T'=Pi, 'BC_left'=[u=0.1*sin(10*t),u_x], 'BC_right'=[u_xx,u_xxx]):
    sol:-animate(size=[600,250], frames=100);


    A dual-span simply-supported beam with moving load


    Load moves across a dual-span beam.

    The beam continues oscillating after the load leaves.

    d := 0.4:  T := 4:  nframes := 100:
    myload := - max(0, -6*(x - t)*(d + x - t)/d^3):
    sol := beamsolve([1,1], 20, 'T'=T, 'F'=myload):
    plots:-animate(plot, [2e-3*myload(x,t), x=0..2, thickness=1, filled=[color="Green"]], t=0..T, frames=nframes):
    plots:-display([%%,%], size=[600,250]);


    A triple-span simply-supported beam with moving load


    Load moves across a triple-span beam.

    The beam continues oscillating after the load leaves.

    d := 0.4:  T := 6: nframes := 100:
    myload := - max(0, -6*(x - t)*(d + x - t)/d^3):
    sol := beamsolve([1,1,1], 20, 'T'=T, 'F'=myload):
    plots:-animate(plot, [2e-3*myload(x,t), x=0..3, thickness=1, filled=[color="Green"]], t=0..T, frames=nframes):
    plots:-display([%%,%], size=[600,250]);



    A triple-span beam, moving load falling off the cantilever end


    In this demo the load move across a multi-span beam with a cantilever section at the right.

    As it skips past the cantilever end, the beam snaps back violently.

    d := 0.4:  T := 8: nframes := 200:
    myload := - max(0, -6*(x - t/2)*(d + x - t/2)/d^3):
    sol := beamsolve([1,1,1/2], 10, 'T'=T, 'F'=myload, BC_right=[u_xx, u_xxx]):
    plots:-animate(plot, [1e-2*myload(x,t), x=0..3, thickness=1, filled=[color="Green"]], t=0..T, frames=nframes):
    plots:-display([%%,%], size=[600,250]);


    Download worksheet:


    4 5 6 7 8 9 10 Last Page 6 of 278