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
  • fyi, there is new video showing Maple's 2025 new interface

    It seems oriented to document mode which I do not use. May be they also improved worksheet mode.

    I am still getting my Maple desktop getting shuffled few times each day where I have to close Maple and reopen it to clear it. I hope they fixed this in the new interface.,

    Hello everyone,

    I have created a Maple worksheet titled "ΕΜΒΑΔΟΝ ΕΠΙΠΕΔΟΥ ΧΩΡΙΟΥ", designed to help my students prepare for their final exams as they qualify for university. This worksheet focuses on area calculations in plane geometry, using Maple to visualize and solve problems efficiently.

    This worksheet is aimed at high school students preparing for university entrance exams, as well as teachers who want to integrate Maple into their teaching.

    I would love to hear your thoughts and feedback!

    Have you used Maple for similar exam preparation?
    εμβαδόν_χωρίου.mw

    MaplePrimes offers spell check and correction

    for the function "Contact Author"

     

    Suggestion:
    For the sake of message clarity (and to save time) it would be desireable to have spell checking and correction in other MaplePrimes message bodies of as well.

     

    Recently, @zenterix asked a question related to solving the quantum mechanics of the finite-wall-height particle-in-a-box problem. In the last two decades or so I've been teaching a quantum mechanics course every second year, in which I sketch the method of solution of this problem for the class, I but do not get into the mathematical details. A few times I've started to work on it in Maple, but abandoned it because of lack of time. So I decided to persist this time.

    This post is more about some of the challenges and tradeoffs in making it work, hence the title. The title is also a nod to the fact that this problem appears in Engel's textbook [1] in a chapter entitled "Applying the particle in a box model to real-world topics". You don't need to know much about quantum mechanics to understand it. From the mathematical point of view you are solving three ordinary differential equations (in three regions) and stitching the solutions together so that the solution, a wavefunction, is continuous with continuous derivative and tends to zero at plus infinity and minus infinity.

    The three regions and the wavefunctions for three eigenvalues (energies) are shown here in the figure, which is the final output of the worksheet. Next is the worksheet and I'll comment further on it below.

    (It's a quirk of quantum mechanics that we plot two things with different units (wavefunctions and energy) on the same axis, and worse, we offset one by the other.)

    Bound states for particle in finite height box of wall height
    V(x) = piecewise(x < -(1/2)*a, V__0, -(1/2)*a <= x and x <= (1/2)*a, 0, (1/2)*a < x, V__0)
    I'll call these regions left, box and right respectively.

    restart

    Schroedinger equation. Here as elsewhere "=0" is implied.

    Schroedinger := -`&hbar;`^2*(diff(psi(x), x, x))/(2*m)+V*psi(x)-E*psi(x)

    -(1/2)*`&hbar;`^2*(diff(diff(psi(x), x), x))/m+V*psi(x)-E*psi(x)

    (1)

    Nondimensionalize length as units of box length, X = x/a. Since psi(x)^2*dx is unitless probability, psi(x) has dimensions 1/length^(1/2)so we define a dimensionless Phi = a^(1/2)*psi

    tr := {x = a*X, psi(x) = Phi(X)/sqrt(a)}; NDSchroedinger := PDETools:-dchange(tr, Schroedinger, [X, Phi(X)], params = [`&hbar;`, m, a], simplify)

    {x = a*X, psi(x) = Phi(X)/a^(1/2)}

     

    -((1/2)*`&hbar;`^2*(diff(diff(Phi(X), X), X))+a^2*m*Phi(X)*(E-V))/(a^(5/2)*m)

    (2)

    This will also change the normalization integrals, e.g., over the box region:

    normint := Int(psi(x)^2, x = -(1/2)*a .. (1/2)*a); NDnormint := PDETools:-dchange(tr, normint, [X, Phi(X)], params = [`&hbar;`, m, a], simplify)

    Int(psi(x)^2, x = -(1/2)*a .. (1/2)*a)

     

    Int(Phi(X)^2, X = -1/2 .. 1/2)

    (3)

    Outside the box we have 0 <= E and E < V with V = V__0. Define a dimensionless positive constant kappa:

    `&kappa;__eqn` := kappa = a*sqrt(2*m*(V__0-E)/`&hbar;`^2); de_outside := simplify(sqrt(a)*kappa^2*(eval(NDSchroedinger, {isolate(`&kappa;__eqn`, m), V = V__0}))/(E-V__0)); outside := rhs(dsolve(de_outside))

    kappa = a*2^(1/2)*(m*(V__0-E)/`&hbar;`^2)^(1/2)

     

    diff(diff(Phi(X), X), X)-kappa^2*Phi(X)

     

    c__1*exp(-kappa*X)+c__2*exp(kappa*X)

    (4)

    Inside the box we have V = 0 and E > 0. Define a dimensionless positive constant k.

    k__eqn := k = a*sqrt(2*m*E/`&hbar;`^2); de_box := simplify(-sqrt(a)*k^2*(eval(NDSchroedinger, {isolate(k__eqn, m), V = 0}))/E); box := rhs(eval(dsolve(de_box), {c__1 = A, c__2 = B}))

    k = a*2^(1/2)*(m*E/`&hbar;`^2)^(1/2)

     

    diff(diff(Phi(X), X), X)+k^2*Phi(X)

     

    A*sin(k*X)+B*cos(k*X)

    (5)

    We have seven unknowns {A, B, E, c__1(L), c__1(R), c__2(L), c__2(R)}, where for example c__1(L) means Maple's c__1 for the left region. We need seven equations: goes to 0 at x = -infinityand at x = infinity(or Phi(X) wouldn't be square integrable), continuity at the two box boundaries, slopes continuous at the two box boundaries, and normalization. We deal with the ones at `&+-`(infinity)first.

    At X = -infinity the wavefunction will be unbounded unless one of the constants goes to zero. We'll rename the other one A__L or A__R. The code here is just to handle the fact that c__1 and c__2 can be the the opposite way around, depending on the session.

    `assuming`(['limit(outside, X = -infinity)' = limit(outside, X = -infinity)], [kappa > 0]); left0 := eval(outside, `~`[`=`](`minus`(indets(rhs(%), name), {infinity}), 0)); left1 := eval(left0, `~`[`=`](indets(left0, suffixed(c)), A__L)); `assuming`(['limit(outside, X = infinity)' = limit(outside, X = infinity)], [kappa > 0]); right0 := eval(outside, `~`[`=`](`minus`(indets(rhs(%), name), {infinity}), 0)); right1 := eval(right0, `~`[`=`](indets(right0, suffixed(c)), A__R))

    limit(c__1*exp(-kappa*X)+c__2*exp(kappa*X), X = -infinity) = signum(c__1)*infinity

     

    A__L*exp(kappa*X)

     

    limit(c__1*exp(-kappa*X)+c__2*exp(kappa*X), X = infinity) = signum(c__2)*infinity

     

    A__R*exp(-kappa*X)

    (6)

    Now we have five unknowns {A, A__L, A__R, B, E}. The strategy will be to eliminate A, A__L, A__R and then find the eigenvalues E. Then we will have expressions for all three regions containing the unknown B, which we will find by the normalization condition (The normalization condition fixes the scale of the wavefunction so that the probability of finding the particle somewhere is one.)

    Continuity of the wavefunction and its derivative at the left boundary are below
    We eliminate A__L and A using this boundary

    bcleft := {eval(left1-box, X = -1/2), eval(diff(left1, X)-(diff(box, X)), X = -1/2)}; sol := solve(bcleft, {A, A__L}); left2 := eval(left1, sol); box2 := eval(box, sol)

    {A__L*exp(-(1/2)*kappa)+A*sin((1/2)*k)-B*cos((1/2)*k), A__L*kappa*exp(-(1/2)*kappa)-A*k*cos((1/2)*k)-B*k*sin((1/2)*k)}

     

    {A = -B*(sin((1/2)*k)*k-cos((1/2)*k)*kappa)/(sin((1/2)*k)*kappa+cos((1/2)*k)*k), A__L = B*k*(sin((1/2)*k)^2+cos((1/2)*k)^2)/(exp(-(1/2)*kappa)*(sin((1/2)*k)*kappa+cos((1/2)*k)*k))}

     

    B*k*(sin((1/2)*k)^2+cos((1/2)*k)^2)*exp(kappa*X)/(exp(-(1/2)*kappa)*(sin((1/2)*k)*kappa+cos((1/2)*k)*k))

     

    -B*(sin((1/2)*k)*k-cos((1/2)*k)*kappa)*sin(k*X)/(sin((1/2)*k)*kappa+cos((1/2)*k)*k)+B*cos(k*X)

    (7)

    Eliminate A__R at the right box boundary

    bcright := {eval(box2-right1, X = 1/2), eval(diff(box2, X)-(diff(right1, X)), X = 1/2)}; elim := eliminate(bcright, A__R); right2 := eval(right1, elim[1]); energy_eqn := elim[2][]/B

    {-B*(sin((1/2)*k)*k-cos((1/2)*k)*kappa)*sin((1/2)*k)/(sin((1/2)*k)*kappa+cos((1/2)*k)*k)+B*cos((1/2)*k)-A__R*exp(-(1/2)*kappa), -B*(sin((1/2)*k)*k-cos((1/2)*k)*kappa)*k*cos((1/2)*k)/(sin((1/2)*k)*kappa+cos((1/2)*k)*k)-B*k*sin((1/2)*k)+A__R*kappa*exp(-(1/2)*kappa)}

     

    [{A__R = B*(2*sin((1/2)*k)*cos((1/2)*k)*kappa+2*cos((1/2)*k)^2*k-k)/(exp(-(1/2)*kappa)*(sin((1/2)*k)*kappa+cos((1/2)*k)*k))}, {-2*(sin((1/2)*k)*cos((1/2)*k)*k^2-sin((1/2)*k)*cos((1/2)*k)*kappa^2-2*cos((1/2)*k)^2*k*kappa+k*kappa)*B}]

     

    B*(2*sin((1/2)*k)*cos((1/2)*k)*kappa+2*cos((1/2)*k)^2*k-k)*exp(-kappa*X)/(exp(-(1/2)*kappa)*(sin((1/2)*k)*kappa+cos((1/2)*k)*k))

     

    -2*sin((1/2)*k)*cos((1/2)*k)*k^2+2*sin((1/2)*k)*cos((1/2)*k)*kappa^2+4*cos((1/2)*k)^2*k*kappa-2*k*kappa

    (8)

    Now we will get the quantized energies by finding which values will make energy_eqn zero. We recall the relationships of `&kappa;__` and k to energy:

    k__eqn; `&kappa;__eqn;`

    k = a*2^(1/2)*(m*E/`&hbar;`^2)^(1/2)

     

    kappa = a*2^(1/2)*(m*(V__0-E)/`&hbar;`^2)^(1/2)

    (9)

    The energy scale is set by V__0, so can get a non-dimensionalized energy e = E/V__0 and a non-dimensional parameter b = a*sqrt(2*m*V__0/`&hbar;`^2).

    b__eqn := b = a*sqrt(2*m*V__0/`&hbar;`^2); e__eqn := e = E/V__0; eqns := `assuming`([{simplify(eval(eval(k__eqn, isolate(e__eqn, E)), isolate(b__eqn, V__0))), simplify(eval(eval(`&kappa;__eqn`, isolate(e__eqn, E)), isolate(b__eqn, V__0)))}], [a > 0])

    b = a*2^(1/2)*(m*V__0/`&hbar;`^2)^(1/2)

     

    e = E/V__0

     

    {k = (e*b^2)^(1/2), kappa = (-b^2*(e-1))^(1/2)}

    (10)

    So now for any values of the box length and height, one finds the parameter b describing the problem and solves for the possible e values.

    energy_eqn2 := `assuming`([simplify((eval(energy_eqn, eqns))/b^2)], [b > 0, e > 0, e < 1])

    2*e^(1/2)*(-e+1)^(1/2)*cos(b*e^(1/2))-2*e*sin(b*e^(1/2))+sin(b*e^(1/2))

    (11)

    Maple cannot find an analytical solution, so we will find the eigenvalues numerically. The zero solution is unphysical.

    solve(energy_eqn2, e)

    0, RootOf(-2*(-(_Z^2-b^2)/b^2)^(1/2)*(_Z^2/b^2)^(1/2)*b^2+2*tan(_Z)*_Z^2-tan(_Z)*b^2)^2/b^2

    (12)

    For any b, read off the possible energies off the vertical line for that b, e.g. there are 3 bound states for b = 9. For other values, choose bval and nsols on the next line appropriately.

    bval := 9; nsols := 3; plots:-implicitplot([b = bval, energy_eqn2], b = 0 .. 15, e = 0 .. 1)

    9

     

    3

     

     

    evals := [fsolve(eval(energy_eqn2, b = bval), e = 0 .. 1, maxsols = nsols)]

    [0.8115199822e-1, .3189598393, .6884466500]

    (13)

    In terms of the parameters b and e and normalization constant B, the non-dimensionalized wavefunctions of the 3 regions are as below. The scale factor gives a simpler form and prevents 0/0 problems later in the calculation.

    scale := `assuming`([denom(simplify(eval(left2, eqns)))], [positive]); left3 := `assuming`([simplify(eval(scale*left2, eqns))], [positive]); box3 := `assuming`([simplify(eval(scale*box2, eqns))], [positive]); right3 := `assuming`([simplify(eval(scale*right2, eqns))], [positive])

    sin((1/2)*b*e^(1/2))*(-e+1)^(1/2)+cos((1/2)*b*e^(1/2))*e^(1/2)

     

    B*e^(1/2)*exp((1/2)*b*(-e+1)^(1/2)*(2*X+1))

     

    -(sin((1/2)*b*e^(1/2))*(e^(1/2)*sin(b*e^(1/2)*X)-(-e+1)^(1/2)*cos(b*e^(1/2)*X))-cos((1/2)*b*e^(1/2))*(cos(b*e^(1/2)*X)*e^(1/2)+sin(b*e^(1/2)*X)*(-e+1)^(1/2)))*B

     

    B*exp(-(1/2)*b*(-e+1)^(1/2)*(-1+2*X))*(e^(1/2)*cos(b*e^(1/2))+(-e+1)^(1/2)*sin(b*e^(1/2)))

    (14)

    It remains to find the normalization constant B. The three parts of the normalization integral are:

    P__L := `assuming`([int(left3^2, X = -infinity .. -1/2)], [positive]); P__box := `assuming`([int(box3^2, X = -1/2 .. 1/2)], [positive]); P__R := `assuming`([int(right3^2, X = 1/2 .. infinity)], [positive])

    (1/2)*B^2*e/(b*(-e+1)^(1/2))

     

    -(1/4)*B^2*(2*e^(1/2)*(-e+1)^(1/2)*cos(2*b*e^(1/2))-2*e^(1/2)*(-e+1)^(1/2)-2*b*e^(1/2)-2*e*sin(2*b*e^(1/2))+sin(2*b*e^(1/2)))/(b*e^(1/2))

     

    (1/2)*B^2*(e^(1/2)*cos(b*e^(1/2))+(-e+1)^(1/2)*sin(b*e^(1/2)))^2/(b*(-e+1)^(1/2))

    (15)

    Solve for B. There are two (messy) solutions for B of opposite sign; it is conventional to choose the sign such that the ground state has positive values.

    Bval := sort([solve(P__L+P__box+P__R = 1, B)])[2]

    Assemble it into a single piecewise function.

    `&Phi;__all` := eval(piecewise(X < -1/2, left3, `and`(X >= -1/2, X <= 1/2), box3, X > 1/2, right3), B = Bval)

    Plot. Plots offset vertically by the energy as usual

    Xmax := 1.5; psiscale := .15; colors := [red, blue, magenta]; wavefns := plot([seq(eval(psiscale*`&Phi;__all`, {b = bval, e = evals[i]})+evals[i], i = 1 .. nsols)], X = -Xmax .. Xmax, color = colors); walls := plot([-1/2, y, y = 0 .. 1], color = black), plot([1/2, y, y = 0 .. 1], color = black), plot(1, X = -Xmax .. -1/2, color = black), plot(0, X = -1/2 .. 1/2, color = black), plot(1, X = 1/2 .. Xmax, color = black); evalplot := plot(evals, X = -Xmax .. Xmax, linestyle = dash, color = colors); plots:-display(wavefns, walls, evalplot, axes = boxed, labels = [X, psiscale*Phi+E/V__0])

     

    NULL

    Download finite_box_non-dim2.mw

    Comments
    1. The first general comment is about the issue of how much Maple can do. In principle one should be able to give dsolve the three ODEs and the boundary conditions and it should output the result. We are still far away from that. Even for one region, the boundary conditions at infinity are not handled by dsolve. The other extreme is to solve the ODEs for their general solutions, and follow the algebraic manipulations for implementing the boundary conditions as one might do it on paper. Engel for example has a comment "At this point we notice that by dividing the equations in each pair, the coefficients can be eliminated to give [...]". Maple will not easily do that trick or the manipulations that led to the special form on which the trick is applied. Levine's textbook [2] gives a different non-intuitive set of steps.

    But Maple can solve complicated equations, so I wanted a more general strategy that avoids as much as possible the requirement to manipulate into special forms. On the other hand, it is tempting (as @zenterix tried) to just give the three general solutions with 6 arbitrary constants and the boundary conditions to Maple's solve, and solve for the six constants. One finds that all six are zero. This us a consequence of the fact that the ODEs are linear, and misses the crucial point that it is an eigenvalue problem. So there must be some sort of step-by-step guidance for Maple and there is a general but non-obvious strategy to solving such problems.

    2. The second general comment is that there is a trade-off between presenting the solution steps without any arcane Maple code, and getting to the solution efficiently and programmatically, without too many manual manipulations. One of my suggestions to @zenterix involved manually dividing both sides of an equations by a fairly complicated expression, to which @acer expressed a preference for programmatic solutions. At the time, I mainly did that to keep close the to existing solution, and was thinking that I usually set things up so that is not necessary. But I found here that I do it a lot, and it is a natural consequence of the interactive way I use Maple. Here's two examples:

    In the second line of label (4), I originally got a more complicated form of de_outside, then manually figured out the factor to put inside simplify to get the form you see. I do this frequently with Maple, especially when I use dchange.

    In the last line of label (8), I manually divided by B to remove it. Had B been in the denominator, I could have removed it programmatically with numer (an operation I use frequently, though it somewhat recklessly ignores denominator zeroes).

    Programmatically doing things can be relatively unreadable and disrupt the readability for the non-Maple reader. For example in label (6) there is "extraneous" code to find which coefficient to set to zero so that the solution goes to zero at +/- infinity. This was forced on me since after generating the nice figure I re-ran the worksheet only to have multiple errors. I originally used eval(outside, {c__1 = B__L, c__2 = A__L}); to change Maple's dsolve constants to read the way I wanted. But I realized that dsolve does not reproducibly assign c__1 and c__2 to the same term each time, raising the general issue of: 

    3. Session to session reproducibility. If I am only going to use a worksheet myself, I don't usually worry too much about this, since I can usually debug this fairly easily. On the other hand if the worksheet is for others, it needs to be foolproof. I have taken to always sorting the output of solve, e.g., the assignment to bval in the execution group after label (15). The code added for the dsolve constants works, but is non-trivial Maple (suffixed type testing). (Solving this issue for Eigenvectors is yet more difficult.)

    4. Some efforts to make the worksheet more readable were hard to do.

    - For some reason, Engel chose to call the two constants in the left region B and B'. Entering B*exp(-k*x) + B'*exp(k*x) in 2D leads to B(x)*exp(-k*x) + diff(B(x), x)*exp(k*x). No doubt this can be fixed, but I didn't pursue this.

    - I would normally use upper case Psi for the non-dimensionalized psi, but Psi is the name of a special function in Maple. But wait, now we can use "local Psi;" to solve this problem. However diff(Psi(x),x,x) displays as Psi(2,x), so I had to settle for Phi. (SCR submitted.)

    - I build up the left wavefunctions incrementally as expressions assigned to left1, left2, etc, which is rather ugly. Why not use arrow procedures so we can use psi(x), D(psi)(x) etc. Aside from some awkwardness in updating such procedures by redefining them as needed, I ran into the following problem:

    limit(A*exp(-k*x) + B*exp(k*x), x = infinity) assuming k > 0; gives as expected
    signum(B)*infinity;

    but

    psi := x -> A*exp(-k*x) + B*exp(k*x);
    limit(psi(x), x = infinity) assuming k > 0;
    returns unevaluated. What happened to full evaluation?

    - I'm used to entering hbar as hbar in 1D; it displays as h with the bar through. In 2-D entering hbar^2 is displayed nicely, but internally it is `&hbar;`. We can have the following puzzle:

    (hbar/m)^2-`&hbar;`^2/m^2

    hbar^2/m^2-`&hbar;`^2/m^2

    NULL

    - I originally wanted a vertical axis label Psi, E/V__0, but labels = [X, Psi, E/V__0] obviously won't work. It took some time to figure out the answer is to make "Psi, E/V__0" an atomic variable. I finally settled on the more honest, if cryptic, label shown.

    5. solve gives an "invalid" solution. Originally I worked with solutions left2 etc (see label (6)) after simplify(eval(left2, eqns)). Every second wavefunction was about 10^9 times higher than the others, which was difficult to debug. It turns out the denominator of those wavefunctions was exactly zero but numerically 10^(-10), meaning they plotted as large-scaled versions of the right shape, without any division-by-zero error. The fix is to remove the denominators by scaling (see label (14)).

    This is just a consequence of Maple giving generic solutions, and the only way around it is to be aware of it.

    Summary

    To solve such problems with Maple requires one to play around in an interactive way. After hacking to a solution, it can take some effort to make it presentable and usable to others. But the final result is pleasing, and is certainly much easier than working through the problem on paper. The ability to manipulate complicated expressions means circumventing tricks that might be needed in manual solutions. 

    References
    [1] T. Engel, Quantum Chemistry and spectroscopy. Pearson 2019, 4th ed. Sec. 5.1, Prob. 5.3. 
    [2] I. N. Levine, Quantum Chemistry. Pearson 2009, 6th ed. Sec. 2.4. 

    To Maplesoft,

    Please consider changing the name from Maple to something else.

    It is almost impossible to search for anything related to maple, since google keeps giving results about trees called Maple in Canada and about some maple syrup products which I have no interest in at all.

    One has to go through pages and pages of links looking for a real Maple software hit.

    I know the name Maple has been around for long time, but a new unique name will make searching easier and people will get used to the new name very quickly (may be in 1-2 years).

    Some examples:

     

    I know when Maple was created almost 40 years ago, the inernet itself was not even here (I forgot when VP. Gore created the internet but I think that was in early 90's), and search was not thought about then.

    But these days, the ability to search for something and to easily find it is very important for companies and having a unique name for Maple will make it also much more popular and easier to find things about it instead of finding  information about Maple syrup and Maple trees all the time.

    This is posted  under product suggestions.

    Hi Maple lovers, and others,

    I wrote up some simple Maple code that 
    checks if an expression (a^3+2) is a prime number.

    The output fits on one page of paper.
    n_cubed_plus_2.mw

    n_cubed_plus_2.pdf

    Also, I have an account at mersenneforum.org

    Kindest Regards,

    Matt

    PS Does this help?  The code works.

    I'm happy to announce the publication of Volume 5, Issue #1 of Maple Transactions.  You can find it at

    mapletransactions.org
     

    We have a survey paper by Veselin Jungic and Naomi Borwein on teaching Experimental Mathematics courses as our Featured Contribution.  Many of you will find it interesting and useful.

    In the refereed paper section we have a paper on Metaprogramming with Maple and C by Ilias Kotsireas; a paper on fast transposed Vandermonde solving by Hyukho Kwon & Michael Monagan; a paper by David Ulgenes (an honours student in Oslo) on Gamma, Pseudogamma, and Inverse Gamma functions; and a paper by John Campbell on applications of Gosper's nonlocal derangement identity (which, if you don't know that the word "derangement" has a technical mathematical meaning, may give you the wrong impression!).

    As usual I've also written something, and I hope you like it: it's about Chladni figures and standing modes in an elliptical drum, and visualizing such in Maple.  It uses Mathieu functions in Maple and noodles a bit about zerofinding (but winds up using fsolve because that's so convenient).
     

    Keep the papers coming.  This is the 12th issue of Maple Transactions, and I remind you that it has a "Diamond Class" designation, which means there are no page charges to authors, and the articles are free to read for everyone.  This means that there's some volunteer labour needed, of course: you have to write the articles, and what we want is that you write articles that people in the Maple community actually want to read.

    I'd also like to thank the copyeditor, Michelle Hatzel, for her very hard work on this issue.  She's really made a difference, and I think you will be able to see it.   

    In the book by W.G. Chinn, N.E. Steenrod "First Concepts of Topology" the another remarkable theorem was proved: any two flat bounded regions can be cut by a single straight line so that each of these regions is divided into two regions of equal area (the second  pancake problem). This is an existence theorem which does not provide any way to find this cut. In this post I made an attempt to find such cut for any 2 convex regions on the plane bounded by a piecewise smooth self-non-intersecting curves.
    The Each_Into_2_Equal_Areas procedure returns a list of coordinates of 4 endpoints of the cutting segments. This procedure significantly uses my old procedures  Area  and  Picture , which can be found in detail at the link  https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-  . The formal arguments of the Each_Into_2_Equal_Areas procedure are the lists  L1  and  L2 specifying the boundaries of the regions to be cut. When specifying  L1  and  L2 , the boundary can be passed clockwise or counterclockwise, but it is necessary that the parameter t (when specifying each link) should go in ascending order. This can always be achieved by replacing  t  with  -t  if necessary. The Pic procedure draws a picture of the source regions and cutting segments. For ease of use, the code for the  Area  and  Picture   procedure is also provided. It is also worth noting that the procedure also works for "not too" non-convex regions (see examples below).

    restart;
    Area := proc(L) 
    local i, var, e, e1, e2, P; 
    for i to nops(L) do 
    if type(L[i], listlist(algebraic)) then 
    P[i] := (1/2)*add(L[i, j, 1]*L[i, j+1, 2]-L[i, j, 2]*L[i, j+1, 1], j = 1 .. nops(L[i])-1) else 
    var := lhs(L[i, 2]); 
    if type(L[i, 1], algebraic) then e := L[i, 1]; 
    if nops(L[i]) = 3 then P[i] := (1/2)*(int(e^2, L[i, 2])) else 
    if var = y then P[i] := (1/2)*simplify(int(e-var*(diff(e, var)), L[i, 2])) else 
    P[i] := (1/2)*simplify(int(var*(diff(e, var))-e, L[i, 2])) end if end if else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; 
    P[i] := (1/2)*simplify(int(e1*(diff(e2, var))-e2*(diff(e1, var)), L[i, 2])) end if end if end do; 
    abs(add(P[i], i = 1 .. nops(L))); 
    end proc:
    
    Picture := proc(L, C, N::posint := 100, Boundary::list := [linestyle = 1]) 
    local i, var, var1, var2, e, e1, e2, P, Q, h; 
    global Border;
    uses plottools; 
    for i to nops(L) do 
    if type(L[i], listlist(algebraic)) then P[i] := op(L[i]) else 
    var := lhs(L[i, 2]); var1 := lhs(rhs(L[i, 2])); var2 := rhs(rhs(L[i, 2])); h := (var2-var1)/N; 
    if type(L[i, 1], algebraic) then e := L[i, 1]; 
    if nops(L[i]) = 3 then P[i] := seq(subs(var = var1+h*i, [e*cos(var), e*sin(var)]), i = 0 .. N) else 
    P[i] := seq([var1+h*i, subs(var = var1+h*i, e)], i = 0 .. N) fi else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; P[i] := seq(subs(var = var1+h*i, [e1, e2]), i = 0 .. N) fi; fi; od; 
    Q := [seq(P[i], i = 1 .. nops(L))]; Border := curve([op(Q), Q[1]], op(Boundary)); [polygon(Q, C), Border] 
    end proc:
    
    Each_Into_2_Equal_Areas:=proc(L1::list, L2::list)
    local D, n, m, L10, L20, S1,S2, f, L11, L21, i, j, k, s, A, B, C , sol;
    
    f:=(X,Y)->expand((y-X[2])*(Y[1]-X[1])-(x-X[1])*(Y[2]-X[2]));
    L10:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L1);
    L20:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L2);
    S1:=Area(L1); S2:=Area(L2);  
    n:=nops(L1); m:=nops(L2);
    
    for i from 1 to n do
    for j from i to n do
    
    for k from 1 to m do
    for s from k to m do
    
    if not ((nops({i,j})=1 and type(L1[i],listlist)) or (nops({k,s})=1 and type(L2[k],listlist))) then
    
    A:=eval(L10[i,1],t=t1): 
    B:=eval(L10[j,1],t=t2):
    C:=eval(L20[k,1],t=t3): 
    D:=eval(L20[s,1],t=t4):
    
    L11:=`if`(j=i,[subsop([2,2]=t1..t2,L10[i]),[B,A]],`if`(j=i+1,[subsop([2,2]=t1..op([2,2,2],L10[i]),L10[i]),subsop([2,2]=op([2,2,1],L10[j])..t2,L10[j]),[B,A]], [subsop([2,2]=t1..op([2,2,2],L10[i]),L10[i]),op(L10[i+1..j-1]),subsop([2,2]=op([2,2,1],L10[j])..t2,L10[j]),[B,A]])):
    
    L21:=`if`(s=k,[subsop([2,2]=t3..t4,L20[k]),[D,C]],`if`(s=k+1,[subsop([2,2]=t3..op([2,2,2],L20[k]),L20[k]),subsop([2,2]=op([2,2,1],L20[s])..t4,L20[s]),[D,C]], [subsop([2,2]=t3..op([2,2,2],L20[k]),L20[k]),op(L20[k+1..s-1]),subsop([2,2]=op([2,2,1],L20[s])..t4,L20[s]),[D,C]])):
    
    sol:=fsolve(simplify({Area(L11)-S1/2,Area(L21)-S2/2,eval(f(A,B),[x=C[1],y=C[2]]),eval(f(A,B),[x=D[1],y=D[2]])}),{t1=op([2,2,1],L10[i])..op([2,2,2],L10[i]),t2=op([2,2,1],L10[j])..op([2,2,2],L10[j]),t3=op([2,2,1],L20[k])..op([2,2,2],L20[k]),t4=op([2,2,1],L20[s])..op([2,2,2],L20[s])});
    
    if type(sol,set(`=`)) then  return eval([A,B,C,D],sol) fi;
    
    fi;
    od: od: od: od:
    end proc:
    
    Pic := proc(L1,L2,col1,col2,Size:=[800,400])
    local P1, P2, P3, T, P;
    uses plots, plottools;
    P1, P2 := Picture(L1, color=col1), Picture(L2, color=col2):
    P3 := line(Sol[1..2][],color=red,thickness=3), line(Sol[3..4][],color=red,thickness=3), line(Sol[1],Sol[4],linestyle=2,thickness=3,color=red):
    T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"]], font=[times,18], align=[left,above]);
    P:=pointplot(Sol, symbol=solidcircle, color=red, symbolsize=10);
    display(P1,P2,P3,T,P, scaling=constrained, size=Size, axes=none);
    end proc: 
    


    Examples of use.

    local D:
    L1:=[[[8,0],[6,7]],[[6,7],[2,5]],[[2,5],[0,2]],[[0,2],[0,0]],[[0,0],[8,0]]]:
    L2:=[[[5*cos(t)+16,5*sin(t)],t=Pi/2..Pi],[[5*cos(t)+16,5*sin(t)/2],t=Pi..2*Pi],[[21,0],[16,5]]]:
    Sol:=Each_Into_2_Equal_Areas(L1,L2): Points:=[A,B,C,D]:
    seq(Points[i]=Sol[i], i=1..4);
    Pic(L1,L2,"Yellow","LightBlue",[900,400]);
    

       

    The specified regions may overlap:

    L1:=[[[8,0],[6,7]],[[6,7],[2,5]],[[2,5],[0,2]],[[0,2],[0,0]],[[0,0],[8,0]]]:
    L2:=[[[5*cos(t)+9,5*sin(t)],t=Pi/2..Pi],[[5*cos(t)+9,5*sin(t)/2],t=Pi..2*Pi],[[14,0],[9,5]]]:
    Sol:=Each_Into_2_Equal_Areas(L1,L2):  Points:=[A,B,C,D]:
    seq(Points[i]=Sol[i], i=1..4);
    Pic(L1,L2,"Yellow","LightBlue");
    

       


    If there is a solution for which the cutting segments intersect the boundary of each of the regions at 2 points, then the procedure also works for such non-convex regions:

    L1:=[[[cos(t),sin(t)],t=Pi/3..2*Pi-Pi/3],[[cos(-t)+1,sin(-t)],t=-Pi-Pi/3..-Pi+Pi/3]]:
    L2:=[[[cos(t)+2,sin(t)],t=-Pi/6..Pi+Pi/6],[[cos(-t)+2,sin(-t)-1],t=-5*Pi/6..-Pi/6]]:
    Sol:=Each_Into_2_Equal_Areas(L1,L2): Points:=[A,B,C,D]:
    seq(Points[i]=Sol[i], i=1..4);
    Pic(L1,L2,"Yellow","LightBlue");
    

       


    A number of other interesting examples can be found in the attached file.

    Each_Into_2_Equal_Areas1.mw

    After a long chat with ChatGPT I finally received a fully working code for a proc for a generalized Woodbury Identity for the inversion of the sum of two or more positive definite matrices.
    I was inspired by a family member who is a trained professional programmer, who told me that in his professional work he uses ChatGTP for an initial draft of his program.  
    I did find out that ChatGTP makes errors: from simple ones like writing 'Simplify' instead of 'simplify' to serious conceptual errors, for example in recursive loops. However, ChatGTP seems to 'understand' the error after given specific feedback. Although, this does not mean that the next proposal does not contain the same logical error. But after a long chat I received a nice proc that seems to work. 
    My second surprise was that Gemini suggested a formula for the generalized Woodbury lemma that was unknown to me, and I was unable to find on Scholar Google or https://math.stackexchange.com. Based on a special case of that formula, I was able to write the second proc myself. 
    In conclusion, to start working it can be helpful to collaborate with AI friend, a little patience may help, AI may not be astute as someone on Mapleprimes wrote, but neither am I. I am now retired, and it is fun to play with Maple and AI. 
    By the way, the search term Woodbury did not give a single hit on Mapleprimes.With_a_little_help_from_my_friends.mw
    kind regards,Harry 

     

    Greetings, dear Maple developers!

    I, Yegor Volovodenko, together with my supervisor, Igor Zinoviev, Associate Professor, Head of the Department of General Mathematics at ZNU, am conducting research on ‘Solving equations and inequalities of elementary mathematics by means of computer mathematics: opportunities, problems and ways to solve them’. I express my sincere gratitude for your work on the development and popularisation of mathematics, for the constant improvement and enhancement of CAS in particular.

    Unfortunately, during the study, we found some, in our opinion, shortcomings in the work of Maxima algorithms for solving elementary equations and inequalities

    1. Solving equations with parameters;



      In these two cases, the solution obtained is formal, without analysing the cases when the coefficient of a variable is zero. Thus, not all solutions of equations with a parameter are obtained. I believe that if the user is not familiar with the methodology for solving equations with parameters, he or she may lose some solutions that may be important for further work. So the solution can be considered incomplete.

     

    1. Solving equations with two variables.
       

    In the equation with two variables, the answer was given in complex numbers, which is incomprehensible to a person who is not familiar with complex numbers. There are also comments on the course of the solution, in this equation it was necessary to select the square of the difference, and then solve the equation x^2+(y-2)^2=0, getting x=0, y=2.

    I hope that the results I have obtained (the identified shortcomings) will help to correct the work of the algorithms and improve the work of the Maple system.

    Sincerely, Yegor Volovodenko, Igor Zinoviev.

     

     

    For Maple 17, in 2013, we introduced the GroupTheory package. It has seen a lot of improvements since its introduction, and I figured I would write something about how you can use it to teach a group theory course.

     

    First of all, I think it would be a great idea to have your students just play with the GroupTheory package in Maple, and get some intuition for what makes group theory tick by getting their hands dirty. But that is not what I want to talk about today. Instead, I want to discuss how you might use it for giving your lectures - to show some visualizations while you show your beamer presentation or write on the black- or whiteboard.

     

    When starting out with the definition of a group, you might want to compare the Cayley table of a group with that of a binary operation that doesn't define a group. A set with any binary operation is called a magma; to find good examples, we might want to use one group; one magma that has an identity and is associative but not invertible; and one magma that is has an identity and is invertible but not associative. Maybe we also want the group to not be cyclic (so the group order will have to not be prime), because those Cayley tables look a little boring. For this we can use the Magma package.

     

    We note that a magma that is invertible and has an identity element is called a loop. We can enumerate all magmas of a given order and with certain properties using the command Magma:-Enumerate, finding either the number of such objects (by default) or a list of the multiplication tables (with the output = list option). Can we find interesting examples of order 4?

     

    with(Magma)

    Enumerate(4, group), Enumerate(4, loop), Enumerate(4, associative, identity)

    2, 2, 35

    (1)

     

    No: all loops of order 4 are groups. Let's try order 6.

     

    Enumerate(6, group), Enumerate(6, loop), Enumerate(6, associative, identity)

    2, 109, 2237

    (2)

     

    Yes, we can find examples here. We find the Cayley tables of three suitable magmas; we call the non-cyclic group one m1, the non-associative loop one m2, and the non-invertible magma m3.

     

    Enumerate(6, group, output = list)

    (3)

    "m1 := ?[2]:"

    Enumerate(6, loop, output = list)[1 .. 5]

    (4)

     

    "m2 := ?[2]:"

    Enumerate(6, associative, identity, output = list)[1 .. 5]

    (5)

     

    Most of these appear to have many more of one value than of another (e.g., many more threes). Let's find one where that is not so extreme. For example, we might want to minimize the standard deviation of the frequencies of the values occurring in the Cayley table.

     

    lst := Enumerate(6, associative, identity, output = list)

    lst := remove(IsLeftInvertible, lst)

    numelems(lst)

    2235

    (6)

    frequencies := map(proc (m) options operator, arrow; map2(numboccur, m, [seq(1 .. 6)]) end proc, lst)``

    with(Statistics)

    sds := map(StandardDeviation, convert(frequencies, 'set'))

    {HFloat(2.449489742783178), HFloat(3.0983866769659336), HFloat(3.1622776601683795), HFloat(3.22490309931942), HFloat(3.286335345030997), HFloat(3.3466401061363027), HFloat(3.4058772731852804), HFloat(3.4641016151377544), HFloat(3.464101615137755), HFloat(3.521363372331802), HFloat(3.5213633723318023), HFloat(3.5777087639996634), HFloat(3.6331804249169903), HFloat(3.687817782917155), HFloat(3.7416573867739413), HFloat(3.794733192202055), HFloat(3.8470768123342687), HFloat(3.847076812334269), HFloat(3.898717737923586), HFloat(3.9496835316262997), HFloat(3.9496835316263), HFloat(3.9999999999999996), HFloat(4.0), HFloat(4.049691346263318), HFloat(4.09878030638384), HFloat(4.147288270665544), HFloat(4.195235392680606), HFloat(4.1952353926806065), HFloat(4.2895221179054435), HFloat(4.33589667773576), HFloat(4.3817804600413295), HFloat(4.427188724235731), HFloat(4.47213595499958), HFloat(4.5166359162544865), HFloat(4.560701700396552), HFloat(4.604345773288535), HFloat(4.604345773288536), HFloat(4.6475800154489), HFloat(4.69041575982343), HFloat(4.732863826479693), HFloat(4.774934554525329), HFloat(4.77493455452533), HFloat(4.8166378315169185), HFloat(4.857983120596447), HFloat(4.898979485566356), HFloat(4.9396356140913875), HFloat(4.939635614091388), HFloat(4.979959839195493), HFloat(5.019960159204453), HFloat(5.059644256269407), HFloat(5.0990195135927845), HFloat(5.138093031466052), HFloat(5.176871642217914), HFloat(5.215361924162119), HFloat(5.253570214625479), HFloat(5.291502622129181), HFloat(5.329165037789691), HFloat(5.366563145999495), HFloat(5.403702434442518), HFloat(5.440588203494177), HFloat(5.477225575051661), HFloat(5.513619500836089), HFloat(5.549774770204643), HFloat(5.585696017507577), HFloat(5.621387729022079), HFloat(5.656854249492381), HFloat(5.692099788303083), HFloat(5.727128425310542), HFloat(5.761944116355173), HFloat(5.796550698475776), HFloat(5.830951894845301), HFloat(5.865151319446072), HFloat(5.8991524815010505), HFloat(5.93295878967653), HFloat(5.932958789676531), HFloat(5.966573556070519), HFloat(6.0), HFloat(6.0332412515993425), HFloat(6.066300355241241), HFloat(6.066300355241242), HFloat(6.099180272790763), HFloat(6.131883886702357), HFloat(6.164414002968976), HFloat(6.196773353931867), HFloat(6.228964600958975), HFloat(6.260990336999411), HFloat(6.29285308902091), HFloat(6.324555320336759), HFloat(6.356099432828282), HFloat(6.418722614352485), HFloat(6.48074069840786), HFloat(6.511528238439883), HFloat(6.542170893518451), HFloat(6.572670690061994), HFloat(6.603029607687671), HFloat(6.603029607687672), HFloat(6.6332495807108), HFloat(6.663332499583073), HFloat(6.6932802122726045), HFloat(6.693280212272605), HFloat(6.723094525588644), HFloat(6.723094525588645), HFloat(6.752777206453653), HFloat(6.782329983125268), HFloat(6.811754546370561), HFloat(6.928203230275509), HFloat(6.957010852370435), HFloat(6.985699678629192), HFloat(7.014271166700073), HFloat(7.042726744663604), HFloat(7.0710678118654755), HFloat(7.127411872482185), HFloat(7.155417527999327), HFloat(7.183313998427188), HFloat(7.18331399842719), HFloat(7.293833011524188), HFloat(7.429670248402684), HFloat(7.4565407529228995), HFloat(7.4565407529229), HFloat(7.483314773547883), HFloat(7.536577472566709), HFloat(7.53657747256671), HFloat(7.563068160475615), HFloat(7.64198926981712), HFloat(7.77174369109018), HFloat(7.8993670632526), HFloat(7.92464510246358), HFloat(7.9498427657407165), HFloat(8.024961059095553), HFloat(8.12403840463596), HFloat(8.366600265340756), HFloat(8.390470785361211), HFloat(8.390470785361213), HFloat(8.414273587185052), HFloat(8.438009243891594), HFloat(8.508818954473059), HFloat(8.854377448471462), HFloat(8.87693640846886), HFloat(8.921883209278185), HFloat(9.338094023943), HFloat(9.338094023943002), HFloat(9.359487165438072), HFloat(9.81835016690686), HFloat(9.818350166906862), HFloat(10.295630140987)}

    (7)

     

    Aha, the smallest standard deviation is a bit under two and a half, and the next possible value is greater than 3. Let's examine the Cayley tables that show this standard deviation.

    NULL

    small_sd_index := select(proc (i) options operator, arrow; StandardDeviation(frequencies[i]) < 3 end proc, [seq(1 .. numelems(frequencies))])

    [1636, 1638, 2234, 2235]

    (8)

    small_sd := lst[small_sd_index]

    (9)

     

    The last of these looks interesting.

     

    m3 := small_sd[-1]

    NULL

    Now we can display these three Cayley tables:

    CayleyColorTable(m1); CayleyColorTable(m2); CayleyColorTable(m3)

     

    It's clear that the multiplication shown in the third Cayley table is not invertible (because multiplying by any element other than 1, on either side, is not a permutation of the elements: rows and columns beyond the first don't have all colors). It's less visually obvious that the second Cayley table doesn't show an associative multiplication, but a quick loop showed that 5*(3*3) = 5*5 and 5*5 = 4 whereas 3*(3*5) = 3 and 3 = 3. I'm not sure if there is something real underpinning this, but to my mind this is a nice parallel to the fact that when you're proving a set with a given operation is a group, it's easier to show that there are inverses than that the operation is associative.

     

    A second nice visualization, also usable in the first week or two of a group theory course, is also named after Cayley: Cayley graphs. Here is a Cayley graph for the dihedral group D[3] of order 6:

     

    with(GroupTheory)

    D3 := DihedralGroup(3)

    _m132248572811840

    (10)

    with(GraphTheory)

    DrawGraph(CayleyGraph(D3), layout = spring)

     

    NULL

    Here is a well-known presentation of the alternating group on 4 letters. I find "grabbing and turning" the 3D plot shows the structure much better than simply an embedding of the graph in the plane.

     

    a4 := `<|>`(`<,>`(a, b), `<,>`(a^3, b^2, a.b.a = b.(a^2).b))

    _m132249105458016

    (11)

    DrawGraph(CayleyGraph(a4), layout = spring, dimension = 3, size = [800, 800])

     

     

    Another great opportunity for a visualization comes a couple of weeks later, when you talk about subgroup lattices. Here are a few examples:

     

    with(GroupTheory)``

    qd := QuasiDihedralGroup(2)

    _m132249751531168

    (12)

    GroupOrder(qd)

    16

    (13)

    DrawSubgroupLattice(qd, normal = false, center = false)

     

     

    This is the plain subgroup lattice. Note how the subgroups in the first and second nontrivial "layer" from the bottom are grouped by conjugacy class. The default visualization of a subgroup lattice highlights the centre in light blue and the (in this case six) other normal subgroups in green.

     

    DrawSubgroupLattice(qd)

     

     

    We see that every conjugacy class of size one is a normal subgroup. We can also highlight other subgroups, if we want. Maybe we want to show which are the (cyclic) subgroups generated by one of the chosen generators:

     

    generators := Generators(qd)

    [_m132249751506336, _m132249751509600]

    (14)

    DrawSubgroupLattice(qd, highlight = [seq({g}, `in`(g, generators))])

     

     

    Of course you can also show more advanced topics in this way. For example, here is the Frattini series of this group.

     

    DrawSubgroupLattice(qd, highlight = FrattiniSeries(qd))

     

    NULL NULL


     

    Download blogpost-algebra.mw

    Recently, a local Math olympiad included the following problem:

     

     

     Compute:Limit(n*(Int(x^n/(x^n+x+2024), x = 0 .. 1)), n = infinity).

     

     

    I was almost sure that Maple can compute it, but unfortunately it needed help!

    restart;

    J:=n->n*int(x^n/(x^n+x+2024),x=0..1);

    proc (n) options operator, arrow; n*(int(x^n/(x^n+x+2024), x = 0 .. 1)) end proc

    (1)

    limit(J(n), n=infinity); # Maple gives up.

    limit(n*(int(x^n/(x^n+x+2024), x = 0 .. 1)), n = infinity)

    (2)

    simplify(IntegrationTools:-Change(J(n), x=t^(1/n), t)) assuming n>1;

    int(t^(1/n)/(t+t^(1/n)+2024), t = 0 .. 1)

    (3)

    ans:=limit(%, n=infinity);

    ln(2)+ln(1013)-4*ln(3)-2*ln(5)

    (4)

    #map(combine,ans, ln);

    combine(ans,ln, integer);

    4*ln((1/15)*2026^(1/4)*5^(1/2))

    (5)

    ln(exp(%)); # should be obtainable with convert

    ln(2026/2025)

    (6)

    evalf(% = J(100)); # sanity check;

    0.4937051076e-3 = 0.488824e-3+0.*I

    (7)


    Download lim-int-vv.mw

    Try it yourself, you will understand what I mean. (MF2024.2.1)

     

    Referring to the screenshots, "J" can be converted to "N m" in MF2024.1, but not in M2024.2.
    Is this some sort of bug in M2024.2?

     

     

    Just around a month after the first release I am glad to announce the second public release of this project.

    Changes from last release:

    • added angled cuts to beam ends
    • fixed bug for bolt connections with thick steel plates
    • rewrite check if fasteners are placed within beams
    • removed some obsolete procedures in NODEFunctions
    • minor changes in XML file headers

    For more information see https://github.com/Anthrazit68/NODEMaple.

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