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
  • Over the last few days, I’ve been creating worksheets on oscillators to support my class’s understanding of these fundamental physics concepts. I wanted to share one of these worksheets that I found particularly useful for illustrating energy exchange and motion dynamics.

    A simple pendulum is a classic physics example that exhibits periodic motion. It consists of a mass m (called a bob) attached to a string or rod of length L, which swings back and forth under the influence of gravity. When the bob is displaced from its equilibrium position and released, it swings back and forth under the influence of gravity.
    To derive the equation of motion, we can examine the forces acting on the pendulum bob and use Newton’s second law.

    Period of a Pendulum:

    • 

    Frequency (f) "-" the number of cycles the pendulum completes in one second. Measued in hertz ("Hz)."

    f = 1/T

    • 

    Period ("T) -" the time it takes the pendulum to complete one cycle. Measued in seconds (s).

    T = 2*Pi*sqrt(L/g)

    This period depends only on the length Land gravitational acceleration "g,"meaning it is independent of the amplitude for small oscillations.

    What is the period and the frequency of a single pendulum that is 70 cm long on the earth and on the moon?

    L := .7; g__earth := 9.8; g__moon := 1.6

    .7

     

    9.8

     

    1.6

    (1)

    T__earth := 2*Pi*sqrt(L/g__earth); f__earth := L/T__earth

    1.679251909

     

    .4168522878

    (2)

    T__moon := 2*Pi*sqrt(L/g__moon); f__moon := L/T__moon

    4.155936442

     

    .1684337597

    (3)

    The above image is taken from https://www.researchgate.net/publication/365297210_Scientific_counterfactuals_as_make-believe

    1. 

    Forces on the Pendulum Bob:

    The main forces acting on the bob are:

    • 

    The gravitational force"`f__g`=mg, "acting vertically downward.

    • 

    The tension Tauin the string, acting along the string toward the pivot point.

    2. 

    Components of the Gravitational Force:

    Since the pendulum swings in an arc, it’s helpful to resolve the gravitational force into two components:

    • 

    Radial Component (along the string): This component, "`f__y`=mgcostheta ," is countered by the tension in the string and does not contribute to the pendulum’s motion.

    • 

    Tangential Component (perpendicular to the string): This component, f__z = -`mgsinθ`(restoring force), acts along the arc of the pendulum’s swing and is responsible for its motion.

    3. 

    Applying Newton's Second Law
    Since the tangential component of the gravitational force causes the pendulum’s motion, we can apply Newton's second law in the tangential direction:``

    f__X = ma__x

    Substituting for f__x and the tangential acceleration a__xNULL

    m*(diff(s(t), t, t)) = -`mgsinθ`

    where diff(s(t), t, t) = a__x and a__x = diff(x, t, t)

    Now, we want to write everything in terms of θ

    s = `Lθ`

    we obtain:

    diff(theta(t), t, t) = -g*`sinθ`/L

    Small-Angle Approximation

    For small angles (typically) theta  , the approximation

    `≈`(sin(theta), theta)(where theta is in radians) can be used. This simplifies the equation:

    diff(theta(t), t, t) = -g*theta/L

    This equation is now in the form of a simple harmonic oscillator

    diff(theta(t), t, t) = -omega^2*theta

    where omega = sqrt(g/L)is the angular frequency of the pendulum.

    restart; with(plots); with(DEtools)

    L := 1; m := .2; g := 9.8

    1

     

    .2

     

    9.8

    (4)

    T := 2*Pi*sqrt(L/g)

    2.007089924

    (5)

    omega := sqrt(g/L)

    3.130495168

    (6)

    ODE__1 := diff(theta(t), t, t)+omega^2*theta = 0; IC := theta(0) = A, (D(theta))(0) = 0

    diff(diff(theta(t), t), t)+9.799999997*theta(t) = 0

     

    theta(0) = A, (D(theta))(0) = 0

    (7)

    sol := dsolve({IC, ODE__1}, theta(t))

    theta(t) = A*cos((1/100000)*97999999970^(1/2)*t)

    (8)

    plot_1 := subs(A = 0.873e-1, sol); plotsresult := plot([rhs(plot_1)], t = 0 .. 2, color = [red])

     

    `θ_t` := rhs(subs(A = 0.873e-1, sol)); v_t := diff(`θ_t`, t)

    0.873e-1*cos((1/100000)*97999999970^(1/2)*t)

     

    -0.8730000000e-6*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

    (9)

    T := (1/2)*m*L^2*v_t^2; V := m*g*L*(1-cos(`θ_t`)); H := T+V

    0.7468864200e-2*sin((1/100000)*97999999970^(1/2)*t)^2

     

    1.96-1.96*cos(0.873e-1*cos((1/100000)*97999999970^(1/2)*t))

     

    0.7468864200e-2*sin((1/100000)*97999999970^(1/2)*t)^2+1.96-1.96*cos(0.873e-1*cos((1/100000)*97999999970^(1/2)*t))

    (10)

    energy_plot := plot([eval(T), eval(V), eval(H)], t = 0 .. 5, color = [red, blue, green], legend = ["Kinetic Energy", "Potential Energy", "Total Energy"], title = "Energy Exchange in Simple Pendulum", labels = ["Time (s)", "Energy (Joules)"])

     

    directionfield := DEplot([diff(theta(t), t) = v(t), diff(v(t), t) = -omega^2*theta(t)], [theta(t), v(t)], t = 0 .. 20, theta = -20 .. 20, v = -40 .. 40, arrows = medium, title = "Direction Field for Simple Harmonic Oscillator", axes = boxed, color = navy)

    sol1 := dsolve({ODE__1, theta(0) = 3, (D(theta))(0) = 0}, theta(t)); sol2 := dsolve({ODE__1, theta(0) = 6.5, (D(theta))(0) = 0}, theta(t)); sol3 := dsolve({ODE__1, theta(0) = -8, (D(theta))(0) = 0}, theta(t)); sol4 := dsolve({ODE__1, theta(0) = 9.7, (D(theta))(0) = 2.5}, theta(t))

    theta(t) = 3*cos((1/100000)*97999999970^(1/2)*t)

     

    theta(t) = (13/2)*cos((1/100000)*97999999970^(1/2)*t)

     

    theta(t) = -8*cos((1/100000)*97999999970^(1/2)*t)

     

    theta(t) = (25000/9799999997)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)+(97/10)*cos((1/100000)*97999999970^(1/2)*t)

    (11)

    theta1 := rhs(sol1); theta2 := rhs(sol2); theta3 := rhs(sol3); theta4 := rhs(sol4)

    3*cos((1/100000)*97999999970^(1/2)*t)

     

    (13/2)*cos((1/100000)*97999999970^(1/2)*t)

     

    -8*cos((1/100000)*97999999970^(1/2)*t)

     

    (25000/9799999997)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)+(97/10)*cos((1/100000)*97999999970^(1/2)*t)

    (12)

    v1 := diff(theta1, t); v2 := diff(theta2, t); v3 := diff(theta3, t); v4 := diff(theta4, t)

    -(3/100000)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

     

    -(13/200000)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

     

    (1/12500)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

     

    (5/2)*cos((1/100000)*97999999970^(1/2)*t)-(97/1000000)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

    (13)

    phase_plot := plot([[eval(theta1, t = tval), eval(v1, t = tval), tval = 0 .. 20], [eval(theta2, t = tval), eval(v2, t = tval), tval = 0 .. 20], [eval(theta3, t = tval), eval(v3, t = tval), tval = 0 .. 20], [eval(theta4, t = tval), eval(v4, t = tval), tval = 0 .. 20]], style = line, title = "Phase Portrait for Simple Harmonic Oscillator", labels = ["x (Displacement)", "v (Velocity)"], color = ["red", "blue", "green", "orange"], axes = boxed)

    display(directionfield, phase_plot)

     

    theta_at_1_sec := subs(t = 1, A = 0.873e-1, rhs(sol)); evalf(theta_at_1_sec)

    0.873e-1*cos((1/100000)*97999999970^(1/2))

     

    -0.8729462437e-1

    (14)
     

    NULL

    Download Simple_Pendulum.mw

    Maple Transactions has just published the Autumn 2024 issue at mapletransactions.org

    From the header:

    This Autumn Issue contains a "Puzzles" section, with some recherché questions, which we hope you will find to be fun to think about.  The Borwein integral (not the Borwein integral of XKCD fame, another one) set out in that section is, so far as we know, open: we "know" the value of the integral because how could the identity be true for thousands of digits but yet not be really true? Even if there is no proof.  But, Jon and Peter Borwein had this wonderful paper on Strange Series and High Precision Fraud showing examples of just that kind of trickery.  So, we don't know.  Maybe you will be the one to prove it! (Or prove it false.)

    We also have some historical papers (one by a student, discussing the work of his great grandfather), and another paper describing what I think is a fun use of Maple not only to compute integrals (and to compute them very rapidly) but which actually required us to make an improvement to a well-known tool in asymptotic evaluation of integrals, namely Watson's Lemma, just to explain why Maple is so successful here.

    Finally, we have an important paper on rational interpolation, which tells you how to deal well with interpolation points that are not so well distributed.

    Enjoy the issue, and keep your contributions coming.

    We have just released updates to Maple and MapleSim.

    Maple 2024.2 includes ability to tear away tabs into new windows, improvements to scrollable matrices, corrections to PDF export, small improvements throughout the math engine, and moreWe recommend that all Maple 2024 users install this update.

    This update also include a fix to the problem with the simplify extension mechanism, as first reported on MaplePrimes. Thanks, as always, for helping us make Maple better.

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

    At the same time, we have also released an update to MapleSim, which contains a variety of improvements to MapleSim and its add-ons. You can find more information on the MapleSim 2024.2 download page.

    Consider the following two systems over the reals:

    sys__1:=(r*(387*r+52)+2<r*(226*q+121*s)+9*q*(q*(2*q-5)-3*s+2)+6*s,4*q**3+r*(27*r+4)+s**2=q*(q+18*r),q>=0,r>=0):
    sys__2:=((392-1739*q)*r+4*(2-9*q)**2+2151*r**2<75*r*s,4*q**3+r*(27*r+4)+s**2=q*(q+18*r),q>=0,r>=0):

    The following results indicate that both and are satisfiable 

    QuantifierElimination:-QuantifierEliminate(exists([s,q,r],And(sys__1)));
                                  true
    
    QuantifierElimination:-QuantifierEliminate(exists([s,q,r],And(sys__1)));
                                  true
    

    but RealDomain:-solve simply returns an empty list (that is, no solution exists) in both cases

    RealDomain:-solve(And(sys__1),[q,s,r]); # thus sys1 cannot be satisfied
                                   []
    
    RealDomain:-solve(And(sys__2),[q,s,r]); # thus sys2 cannot be satisfied
                                   []
    

    Evidently, the above two results conflict with each other, which suggests that there must be a bug in one of these two functions. 
    As discussed in the previous problem, the use of RealDomain:-solve is unsafe, Nevertheless, it appears that using QuantifierElimination:-QuantifierEliminate is still not always safe (since the unsatisfiability is also checked by another software). 

    Updated. The “inconsistency” problem still exists in Maple 2025: 
     

    restart;

    kernelopts('version');

    `Maple 2025.1, X86 64 WINDOWS, Jun 12 2025, Build ID 1932578`

    (1)

    `#msub(mi("sys"),mi("1"))` := r*(387*r+52)+2 < r*(226*q+121*s)+9*q*(q*(2*q-5)-3*s+2)+6*s, 4*q^3+r*(27*r+4)+s^2 = q*(q+18*r), q >= 0, r >= 0

    `#msub(mi("sys"),mi("2"))` := (392-1739*q)*r+4*(2-9*q)^2+2151*r^2 < 75*r*s, 4*q^3+r*(27*r+4)+s^2 = q*(q+18*r), q >= 0, r >= 0

    RegularChains:-SemiAlgebraicSetTools:-QuantifierElimination(`&E`([q, r]), `&and`(sys__1))

    RegularChains:-SemiAlgebraicSetTools:-QuantifierElimination(`&E`([q, r]), `&and`(sys__2))

    false

     

    false

    (2)

    QuantifierElimination:-QuantifierEliminate(exists([q, r], And(sys__1)))

    QuantifierElimination:-QuantifierEliminate(exists([q, r], And(sys__2)))

    s-RootOf(27*_Z^3-216*_Z^2+9*_Z+1, 17498784809929/2199023255552 .. 139990278479459/17592186044416) = 0

     

    s-RootOf(13500*_Z^2-210681*_Z-9826, 144368837035598577819/9223372036854775808 .. 288737674071197155665/18446744073709551616) = 0

    (3)

    NULL


     

    Download 239277-Which-Function-Is-More-Credible.mw

    With the 2024 Maple Conference coming up this week, I imagine one of two of you have noticed something missing. We chose not to have a Conference Art Gallery this year, because we have been working to launch new section of MaplePrimes:  The MaplePrimes Art Gallery. This new section of MaplePrimes is designed for showing off your Maple related images, in a gallery format that puts the images up front, like Instagram but for Math.

    To get the ball rolling on the gallery, I have populated it with many of the works from past years' Maple Art Galleries, so you can now browse them all in one place, and maybe give "Thumbs Ups" to art pieces that you felt should have won the coveted "People's Choice Award".

    Once you are inspired by past entries, we welcome you to submit your new artwork!  Just like the conference galleries, we are looking for all sorts of mathematical art ranging from computer graphics and animations, to photos of needlework, geometrical sculptures, or almost anything!  Art Gallery posts are very similar to regular MaplePrimes posts except that you are asked to supply an Image File and a Caption that will displayed on the main Gallery Page, the post itself can include a longer description, Maple commands, additional images, and whatever else you like.  See for example this Art Gallery post describing Paul DeMarco's sculpture from the 2021 Maple Conference Gallery - it includes an embedded worksheet as well as additional images.

    I can't wait to see what new works of art our MaplePrimes community creates!

     

    I recently prepared a worksheet to teach vector fundamentals in one of my classes, and I wanted to share it with you all. It's nothing special, but I found Maple really helpful in demonstrating the concepts visually. Below is a breakdown of what the worksheet covers, with some Maple code examples included.

    Feel free to take a look and use it if you find it useful! Any feedback or suggestions on how to improve it would be appreciated.

    restart

    NULL

    v := `<|>`(`<,>`(2, 3)); w := `<|>`(`<,>`(4, 1))

    Matrix(%id = 36893488152076804092)

    (1)

    Basic Vector Operations

    • 

    Addition and Subtraction

     

    We  can add and subtract vectors easily if they are of the same dimension.

    NULL

    u_add := v+w; u_sub := v-w

    Matrix(%id = 36893488152076803596)

    (2)

    NULL

    NULL

    Typesetting[delayDotProduct]((((Triangle(L)*a*w*o*f*A*d*d)*i*t*i)*o*n*o)*f, Vector(s), true)

    "The famous triangle law can be used for the addition of vectors and this method is also called the head-to-tail method,As per this law,two vectors can be added together by placing them together in such away that the first vector's head joins the tail of the second vector. Thus,by joining the first vector's tail to the head of the second vector,we can obtain the resultant vector sum."

    NULL

    with(plots); display(arrow([0, 0], [2, 3], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([2, 3], [4, 1], color = green, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [6, 4], difference = true, color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), scaling = constrained, labels = [x, y], title = "
Triangle Law of Addition of Vectors", titlefont = [times, 20, bold])

     

    NULL

    NULL

    Parallelogram Law of Addition of Vectors

    "An other law that can be used for the addition of vectors is the parallelogram law of the addition of vectors*Let's take two vectors v and u,as shown below*They form the two adjacent sides of aparallelogram in their magnitude and direction*The sum v+u is represented in magnitude and direction by the diagonal of the parallelogram through thei rcommon point."

    NULL

    display(arrow([0, 0], [2, 3], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [6, 4], color = green, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [4, 1], difference = true, color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), scaling = constrained, labels = [x, y], title = "

Parallelogram Law of Addition of Vectors", titlefont = [times, 20, bold])

     

    NULL

    NULL

    NULL

    NULL

    NULL

    NULL

    NULL

    • 

    Scalar Multiplication

    We can multiply a vector by a scalar. To multiply a vector by a scalar (a constant), multiply each of its components by the constant.

    Suppose we let the letter  λ represent a real number and  u⃗ = (x,y) be the vector

    v_scaled := 3*v

    Matrix(%id = 36893488152152005076)

    (3)

    NULL

    NULL

    • 

    -Define*the*opposite*vector*v

    Error, (in LinearAlgebra:-Multiply) invalid arguments

     

    Two vectors are opposite if they have the same magnitude but opposite direction.

    v_opposite := -v; vec1 := arrow([0, 0], [2, 3], color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1); vec2 := arrow([0, 0], [-2, -3], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1); display([vec1, vec2], scaling = constrained, title = "Original Vector and its Opposite (-v)")

     
    • 

    Dot Product

    Sometimes the dot product is called the scalar product. The dot product is also an example of an inner product and so on occasion you may hear it called an inner product.

    Given the two vectors  `#mover(mi("a"),mo("&rarr;"))` = (x__1, y__1) and `#mover(mi("b"),mo("&rarr;"))` = (x__2, y__2)
     the dot product is, `#mover(mi("a"),mo("&rarr;"))`.`#mover(mi("b"),mo("&rarr;"))` = x__1*x__2+y__1*y__2

    `#mover(mi("a"),mo("&rarr;"))`.`#mover(mi("b"),mo("&rarr;"))` = x__1*x__2+y__1*y__2

    (4)

    `#mover(mi("a"),mo("&rarr;"))` := `<|>`(`<,>`(5, -8))

    Matrix(%id = 36893488152255219820)

    (5)

    `#mover(mi("b"),mo("&rarr;"))` := `<|>`(`<,>`(1, 2))

    Matrix(%id = 36893488152255215244)

    (6)


    display(arrow([0, 0], [5, -8], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [1, 2], difference = true, color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), scaling = constrained, labels = [x, y])

     

     

     

     

     

     

     

     

     

     

     

      with(LinearAlgebra)

    dot_product := DotProduct(`#mover(mi("a"),mo("&rarr;"))`, `#mover(mi("b"),mo("&rarr;"))`)

    -11

    (7)
    • 

    Vector Norm (Magnitude)

    To find the magnitude (or length) of a vector, use Norm.

    norm_a := Norm(`#mover(mi("a"),mo("&rarr;"))`, 2); norm_b := Norm(`#mover(mi("b"),mo("&rarr;"))`, 2)

    5^(1/2)

    (8)
    • 

    Calculate the Cosine Between Two Vectors

    cos_theta := dot_product/(norm_a*norm_b)

    -(11/445)*89^(1/2)*5^(1/2)

    (9)

    angle_radians := arccos(cos_theta)

    Pi-arccos((11/445)*89^(1/2)*5^(1/2))

    (10)

    angle_degrees := evalf(convert(angle_radians, degrees))

    121.4295656*degrees

    (11)

    NULL

    We can determine whether two vectors are parallel by using the scalar multiple method or the determinant (area of the parallelogram formed by the vectors) method.

    ``

    • 

    Scalar Multiple Method

    a := Vector([5, -8]); b := Vector([10, -16]); k := a[1]/b[1]; is_parallel := a[2]/b[2] = k

    1/2 = 1/2

    (12)
    • 

    Determinant Method

    determinant := a[1]*b[2]-a[2]*b[1]; result := determinant = 0

    0 = 0

    (13)

     

     

    display(arrow([0, 0], [5, -8], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [10, -16], difference = true, color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), scaling = constrained, labels = [x, y], title = "Parallel Vectors")

     

     

     

     

    Two vectors a and b are perpendicular (vertical)

    NULL

    • 

    if and only if their dot product is zero

    a1 := Vector([1, 2]); b1 := Vector([-2, 1]); dot_product := DotProduct(a1, b1); is_perpendicular := dot_product = 0

    0 = 0

    (14)

    display(arrow([0, 0], [-2, 1], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [1, 2], difference = true, color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), scaling = constrained, labels = [x, y], title = "Perpendicular Vectors")

     
    • 

    Slope Method

    slope_a := a1[2]/a1[1]; slope_b := b1[2]/b1[1]; is_perpendicular := slope_a*slope_b = -1

    -1 = -1

    (15)

    NULL

    • 

    Special case: If one vector is vertical (undefined slope) and the other is horizontal (zero slope), they are perpendicular.

    a2 := Vector([0, 3]); b2 := Vector([4, 0]); is_perpendicular := a2[1] = 0 and b2[2] = 0

    true

    (16)

    vec1 := arrow([0, 0], [4, 0], color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1); vec2 := arrow([0, 0], [0, 3], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1); display([vec1, vec2], scaling = constrained, title = "Perpendicular Vectors")

     

    NULL

    Download vectors.mw

    The tab key, and the mouse can be used to trigger completion selection in Maple 2024.1 like previous versions. 

    For interface resposiveness, a new option (disabled by default) has been added to the interface tab of Maple 2024.1 Options (available from the Tools menu). Users that prefer to use enter to trigger completion selection can enable the option. 

    Change the option here and apply to the session or globally if using enter to trigger completion selection is preferred.


    This is an answer to the last question ASHAN recently asked.
    As this answer goes beyond the simple example supplied by ASHAN in that it extends the possibilities of Statistics:-ColumnGraph, I thought it would be appropriate to share it with the entire community of Maple users.

    Possible improvements:

    • add controls for parameter types
    • add consistency controls (numer of colors = number of data rows for instance)
    • add some custom parameters
    • improve legend position
       

    restart

    CustomColumnGraph := proc(
                           data
                           , BarColors
                           , GroupTitles
                           , BarWidth         
                           , GroupOffset      
                           , LegendBarHeight
                           , VerticalShiftFactor
                           , BarFont
                           , GroupFont
                           , LegendFont
                           , BarTitleFormat
                         )



      local Ni, Nj, W, DX, LH, LW, SF, L, DL:

      local SingleBar := (w, h, c) -> rectangle([0, 0], [w, h], color=c):
      local BarGraph, Legends:

      uses plots, plottools:

      Ni := numelems(data[..,1]):
      Nj := numelems(data[1]):
      W  := BarWidth:
      DX := GroupOffset:
      LH := LegendBarHeight:
      SF := VerticalShiftFactor;
      L  := DX*(Ni-1) + W*(Nj-1) + DX/2 + DX/4:
      LW := L / Ni / Nj / 2;                          # LegendBarWidth  
      DL := solve(dl*(Nj+1) + Nj*(3*LW) = L, dl);



      BarGraph := display(
                    seq(
                      seq(
                        translate(SingleBar(W, data[i, j], BarColors[j]), DX/4 + DX*(i-1) + W*(j-1), 0)
                        , i=1..Ni
                      )
                      , j=1..Nj
                    )
                    , axis[1] = [gridlines = [], tickmarks = []]
                    , axis[2] = [gridlines = [11, linestyle = dot, thickness = 0, color = "LimeGreen"]]
                  ):

      Legends := display(
                   seq(
                     translate(SingleBar(LW, LH, BarColors[j]), DL/2 + LW + (3*LW+DL)*(j-1), min(data)*(2*SF-1))
                     , j=1..Nj
                   )
                   ,
                   seq(
                     translate(
                       textplot([0, LH/2, GroupTitles[j]], align=right, font = LegendFont)
                       , DL/2 + 2*LW + (3*LW+DL)*(j-1), min(data)*1.4)
                       , j=1..Nj
                   )
                 ):

      display(
        BarGraph  
        , Legends

        #------------------------------- X base line
        , plot(0, x=0..DX*(Ni-1)+W*(Nj-1)+DX/2 + DX/4)
     
        #------------------------------- Left bar group boundaries
        ,
        seq(
          line(
            [DX/4 + DX*(i-1) - W/2, min(data)*SF],
            [DX/4 + DX*(i-1) - W/2, max(data)*SF]
            , color="LightGray"
          )
          , i=1..Ni
        )


        #------------------------------- Right bar group boundaries
        ,
        seq(
          line(
            [DX/4 + DX*(i-1) + W*Nj + W/2, min(data)*SF],
            [DX/4 + DX*(i-1) + W*Nj + W/2, max(data)*SF]
            , color="LightGray"
          )
          , i=1..Ni
        )
      
        #------------------------------- Top bar group boundaries
        ,
        seq(
          line(
            [DX/4 + DX*(i-1) - W/2, max(data)*SF],
            [DX/4 + DX*(i-1) + W*Nj + W/2, max(data)*SF]
            , color="LightGray"
          )
          , i=1..Ni
        )
     
        #------------------------------- Bar titles
        ,
        seq(
          seq(
            textplot(
              # [DX/4 + DX*(i-1) + W*(j-1) +W/2, data[i, j], evalf[3](data[i, j])]
              [DX/4 + DX*(i-1) + W*(j-1) +W/2, data[i, j], sprintf(BarTitleFormat, data[i, j])]
              , font=BarFont
              , align=`if`(data[i, j] > 0, above, below)
            )
            , i=1..Ni
          )
          , j=1..Nj
        )
      
        #------------------------------- Bar group titles
        ,
        seq(
          textplot(
            [DX/4 + DX*(i-1) + W*(Nj/2), max(data)*SF, cat("B", i)]
            , font=GroupFont
            , align=above
          )
          , i=1..Ni
        )
      
        , labels=["", ""]
        , size=[800, 400]
        , view=[0..DX*(Ni-1) + W*(Nj-1) + DX/2 + DX/4, default]
      );
    end proc:

     

    Data from  AHSAN

     

    A := [0.1, 0.5, 0.9]:
    B := [0.5, 2.5, 4.5]:


    Nu_A := [-0.013697, 0.002005, 0.017707, -0.013697, 0.002005, 0.017707, -0.013697, 0.002005, 0.017707]:
    Nu_B := [0.010827, 0.011741, 0.012655, 0.013967, 0.014881, 0.015795, 0.017107, 0.018021, 0.018935]:


    i:=4:
    data := Matrix([[Nu_A[i], Nu_B[i]], [Nu_A[i + 1], Nu_B[i + 1]], [Nu_A[i + 2], Nu_B[i + 2]]]):
     

     

    Column Graph

     

    Ni, Nj := LinearAlgebra:-Dimensions(data):

    CustomColumnGraph(
                       data                                        # data
                       , ["Salmon", "LightGreen"]                  # BarColors
                       , ["Sensitivity to A", "Sensitivity to B"]  # GroupTitles
                       , 1/(numelems(data[1])+2)                   # BarWidth  
                       , 1                                         # GroupOffset  
                       , max(abs~(data))/10                        # LegendBarHeight
                       , 1.2                                       # VerticalShiftFactor
                       , [Tahoma, bold, 10]                        # BarFont
                       , [Tahoma, bold, 14]                        # GroupFont
                       , [Tahoma, bold, 12]                        # LegendFont
                       , "%1.5f"                                   # BarTitleFormat
                     )

     

     

    Column Graph of an extended set of data

     

     

    ExtentedData := < data | Statistics:-Mean(data^+)^+ >:
    ExtentedData := < ExtentedData, Statistics:-Mean(ExtentedData) >:

    interface(displayprecision=6):
    data, ExtentedData;

    Ni, Nj := LinearAlgebra:-Dimensions(ExtentedData):


    CustomColumnGraph(
                       ExtentedData                                  
                       , ["Salmon", "LightGreen", "Moccasin"]                 
                       , ["Sensitivity to A", "Sensitivity to B", "Mean Sensitivity to A and B"]
                       , 1/(numelems(data[1])+2)                   # BarWidth  
                       , 1                                         # GroupOffset  
                       , max(abs~(data))/10                        # LegendBarHeight
                       , 1.2                                       # VerticalShiftFactor
                       , [Tahoma, 10]                              # BarFont
                       , [Tahoma, bold, 14]                        # GroupFont
                       , [Tahoma, 12]                              # LegendFont
                       , "%1.3f"                                   # BarTitleFormat
                     )

    Matrix(3, 2, {(1, 1) = -0.13697e-1, (1, 2) = 0.13967e-1, (2, 1) = 0.2005e-2, (2, 2) = 0.14881e-1, (3, 1) = 0.17707e-1, (3, 2) = 0.15795e-1}), Matrix(4, 3, {(1, 1) = -0.13697e-1, (1, 2) = 0.13967e-1, (1, 3) = HFloat(1.349999999999997e-4), (2, 1) = 0.2005e-2, (2, 2) = 0.14881e-1, (2, 3) = HFloat(0.008442999999999999), (3, 1) = 0.17707e-1, (3, 2) = 0.15795e-1, (3, 3) = HFloat(0.016751000000000002), (4, 1) = HFloat(0.002005), (4, 2) = HFloat(0.014881), (4, 3) = HFloat(0.008443)})

     

     

     

     


     

    Download CustomColumnGraph.mw

     

    This post is about the visualization of a gyroscopic phenomenon of a rotating body. MapleSim models and a description for those who do not have MapleSim are provided for their own analysis. Implementation with other tools like Maple might give further insight into the phenomenon.

    With appropriate initial conditions, a ball thrown into a tube can pop out of the tube. This can be reproduced with a MapleSim model

    Throwing_a_ball_into_a_tube_A.msim

    To hit a perfect shot without trial and error, time reversal was applied for the model (reversed calculation results of a ball exiting the tube are used as initial conditions for the shot). This worked straight away and shows that this model is sufficiently conservative.

    This phenomenon has recently attracted attention on YouTube. For example, Steve Mold demonstrates the effect and provides an intuitive explanation which he considers incomplete because the resulting vertical oscillation of the ball does not match theory and his experiments. He suspects that the assumption of a constant axis of rotation of the ball is responsible for this discrepancy.

    However, he cannot demonstrate a change of the axis of rotation. In general, the visualization of the rotation axis of a ball is difficult to achieve in an experiment. On the contrary, visualization is much easier in a simulation experiment with this model:

    Throwing_a_ball_into_a_tube_B.msim

    The following can be observed for a trajectroy that does not exit the tube:

    At the apex (the top) of the trajectory, the vector of rotation (red bold in the following images) points downwards and is essentially parallel to the axis of the cylinder. The graph to the left shows the vertical (in green) position and one horizontal position (in red). The model applies gravity in negative y direction.

    Ein Bild, das Text, Diagramm, Screenshot, Reihe enthält.

Automatisch generierte Beschreibung 

    On the way down, the axis of rotation points away from the direction of travel (the ball orbits counterclockwise in the top view).

    Ein Bild, das Text, Diagramm, Screenshot, Reihe enthält.

Automatisch generierte Beschreibung

    At the bottom, the vector of rotation points towards the axis of the cylinder.

    Ein Bild, das Text, Diagramm, Screenshot, Reihe enthält.

Automatisch generierte Beschreibung

    On the way up, the axis of rotation points in the direction of travel.

    Ein Bild, das Text, Diagramm, Screenshot, Reihe enthält.

Automatisch generierte Beschreibung

    These observations confirm that the assumption of a constant axis of rotation is too simplified. Effectively the ball performs a precession movement know from gyroscopes. More specifically, the precession movement of the rotation axis rotates in the opposite direction of the rotation of the ball.

    However, the knowledge and the visualization of this precession movement do not provide more insight for a better intuitive explanation of the effect. As the ball acts like a gyroscope, a second attempt is to visualize forces that perturb the motion of the ball. Besides gravity, there are contact forces exerted by the tube. The normal force at the contact as well as the gravitational force cannot generate a perturbing momentum since they point to the center of the ball. Only frictional forces at the contact can cause a perturbing momentum.

    Contrary to the visualization of the axis of rotation, visualization of contact forces is not straight forward in MapleSim, because neither the contact point nor the contact forces are directly provided by components of the MapleSim library. Only for a single contact point, a work-around is possible by measuring the reactive forces on the tube and then displaying these forces in a moving reference frame at the contact point. The location and the orientation of this frame are calculated with built-in mathematical components. To illustrate the additional effort, the image below highlights in yellow the components only needed for the visualization of the above images, all other components were required to visualize the contact forces and frictional moments.