Download Bug-solve_ineqs.mw

This Maplesoft guest blog post is from Prof. Dr. Johannes Blümlein from Deutsches Elektronen-Synchrotron (DESY), one of the world’s leading particle accelerator centres used by thousands of researchers from around the world to advance our knowledge of the microcosm. Prof. Dr. Blümlein is a senior researcher in the Theory Group at DESY, where he and his team make significant use of Maple in their investigations of high energy physics, as do other groups working in Quantum Field Theory. In addition, he has been involved in EU programs that give PhD students opportunities to develop their Maple programming skills to support their own research and even expand Maple’s support for theoretical physics.

The use of Maple in solving frontier problems in theoretical high energy physics

For several decades, progress in theoretical high energy physics relies on the use of efficient computer-algebra calculations. This applies both to so-called perturbative calculations, but also to non-perturbative computations in lattice field theory. In the former case, large classes of Feynman diagrams are calculated analytically and are represented in terms of classes of special functions. In early approaches started during the 1960s, packages like Reduce [1] and Schoonship [2] were used. In the late 1980s FORM [3] followed and later on more general packages like Maple and Mathematica became more and more important in the solution of these problems. Various of these problems are related to data amounts in computer-algebra of O(Tbyte) and computation times of several CPU years currently, cf. [4].

Initially one has to deal with huge amounts of integrals. An overwhelming part of them is related by Gauss’ divergence theorem up to a much smaller set of the so-called master integrals (MIs). One performs first the reduction to the MIs which are special multiple integrals. No general analytic integration procedures for these integrals exist. There are, however, several specific function spaces, which span these integrals. These are harmonic polylogarithms, generalized harmonic polylogarithms, root-valued iterated integrals and others. For physics problems having solutions in these function spaces codes were designed to compute the corresponding integrals. For generalized harmonic polylogarithms there is a Maple code HyperInt [5] and other codes [6], which have been applied in the solution of several large problems requiring storage of up to 30 Gbyte and running times of several days. In the systematic calculation of special numbers occurring in quantum field theory such as the so-called β-functions and anomalous dimensions to higher loop order, e.g. 7–loop order in Φ⁴ theory, the Maple package HyperLogProcedures [7] has been designed. Here the largest problems solved require storage of O(1 Tbyte) and run times of up to 8 months. Both these packages are available in Maple only.

A very central method to evaluate master integrals is the method of ordinary differential equations. In the case of first-order differential operators leading up to root-valued iterative integrals their solution is implemented in Maple in [8] taking advantage of the very efficient differential equation solvers provided by Maple. Furthermore, the Maple methods to deal with generating functions as e.g. gfun, has been most useful here. For non-first order factorizing differential equation systems one first would like to factorize the corresponding differential operators [9]. Here the most efficient algorithms are implemented in Maple only. A rather wide class of solutions is related to 2^nd order differential equations with more than three singularities. Also here Maple is the only software package which provides to compute the so-called ₂F₁ solutions, cf. [10], which play a central role in many massive 3-loop calculations

The Maple-package is intensely used also in other branches of particle physics, such as in the computation of next-to-next-to leading order jet cross sections at the Large Hadron Collider (LHC) with the package NNLOJET and double-parton distribution functions. NNLOJET uses Maple extensively to build the numerical code. There are several routines to first build the driver with automatic links to the matrix elements and subtraction terms, generating all of the partonic subprocesses with the correct factors. To build the antenna subtraction terms, a meta-language has been developed that is read by Maple and converted into calls to numerical routines for the momentum mappings, calls to antenna and to routines with experimental cuts and plotting routines, cf. [11].

In lattice gauge calculations there is a wide use of Maple too. An important example concerns the perturbative predictions in the renormalization of different quantities. Within different European training networks, PhD students out of theoretical high energy physics and mathematics took the opportunity to take internships at Maplesoft for several months to work on parts of the Maple package and to improve their programming skills. In some cases also new software solutions could be obtained. Here Maplesoft acted as industrial partner in these academic networks.

References

[1] A.C. Hearn, Applications of Symbol Manipulation in Theoretical Physics, Commun. ACM 14 No. 8, 1971.

[2] M. Veltman, Schoonship (1963), a program for symbolic handling, documentation, 1991, edited by D.N. Williams.

[3] J.A.M. Vermaseren, New features of FORM, math-ph/0010025.

[4] J. Blümlein and C. Schneider, Analytic computing methods for precision calculations in quantum field theory, Int. J. Mod. Phys. A 33 (2018) no.17, 1830015 [arXiv:1809.02889 [hep-ph]].

[5] E. Panzer, Algorithms for the symbolic integration of hyperlogarithms with applications to Feynman integrals, Comput. Phys. Commun. 188 (2015) 148–166 [arXiv:1403.3385 [hep-th]].

[6] J. Ablinger, J. Blümlein, C .Raab, C. Schneider and F. Wissbrock, Calculating Massive 3-loop Graphs for Operator Matrix Elements by the Method of Hyperlogarithms, Nucl. Phys. 885 (2014) 409-447 [arXiv:1403.1137 [hep-ph]].

[7] O. Schnetz, φ⁴ theory at seven loops, Phys. Rev. D 107 (2023) no.3, 036002 [arXiv: 2212.03663 [hep-th]].

[8] J. Ablinger, J. Blümlein, C. G. Raab and C. Schneider, Iterated Binomial Sums and their Associated Iterated Integrals, J. Math. Phys. 55 (2014) 112301 [arXiv:1407.1822 [hep-th]].

[9] M. van Hoeij, Factorization of Differential Operators with Rational Functions Coefficients, Journal of Symbolic Computation, 24 (1997) 537–561.

[10] J. Ablinger, J. Blümlein, A. De Freitas, M. van Hoeij, E. Imamoglu, C. G. Raab, C. S. Radu and C. Schneider, Iterated Elliptic and Hypergeometric Integrals for Feynman Diagrams, J. Math. Phys. 59 (2018) no.6, 062305 [arXiv:1706.01299 [hep-th]].

[11] A. Gehrmann-De Ridder, T. Gehrmann, E.W.N. Glover, A. Huss and T.A. Morgan, Precise QCD predictions for the production of a Z boson in association with a hadronic jet, Phys. Rev. Lett. 117 (2016) no.2, 022001 [arXiv:1507.02850 [hep-ph]].

     Happy Easter to all those who celebrate! One common tradition this time of year is decorating Easter eggs. So, we’ve decided to take this opportunity to create some egg-related math content in Maple Learn. This year, a blog post by Tony Finch inspired us to create a walkthrough exploring the four-point egg. The four-point egg is a method to construct an egg-shaped graph using just a compass and a ruler, or in this case, Maple Learn. Here's the final product: 

     The Maple Learn document, found here, walks through the steps. In general, each part of the egg is an arc corresponding to part of a circle centred around one of the points generated in this construction. 

     For instance, starting with the unit circle and the three red points in the image below, the blue circle is centred at the bottom point such that it intersects with the top of the unit circle, at (0,1). The perpendicular lines were constructed using the three red points, such that they intersect at the bottom point and pass through opposite side points, either (-1,0) or (1,0). Then, the base of the egg is constructed by tracing an arc along the bottom of the blue circle, between the perpendicular lines, shown in red below.

     Check out the rest of the steps in the Maple Learn Document. Also, be sure to check out other egg-related Maple Learn documents including John May’s Egg Formulas, illustrating other ways to represent egg-shaped curves with mathematics, and Paige Stone’s Easter Egg Art, to design your own Easter egg in Maple Learn. So, if you’ve had your fill of chocolate eggs, consider exploring some egg-related geometry - Happy Easter!  

Here are the Maple sources for the paper "Maximum Gap among Integers Having a Common Divisor with an Odd semi-prime" published on Journal of Advances in Mathematics and Computer Science 39 (10):51-61, at :https://doi.org/10.9734/jamcs/2024/v39i101934.

gap := proc(a, b) return abs(a - b) - 1; end proc

HostsNdivisors := proc(N)

local i, j, g, d, L, s, t, m, p, q, P, Q, np, nq;

m := floor(1/2*N - 1/2);

L := evalf(sqrt(N));

P := Array();

Q := Array();

s := 1; t := 1;

for i from 3 to m do

   d := gcd(i, N);

    if 1 < d and d < L then P(s) := i; s++;

    elif L < d then Q(t) := i; t++; end if;

end do;

  np := s - 1;

  nq := t - 1;

 for i to np do printf("%3d,", P(i)); end do;

  printf("\n");

  for i to nq do printf("%3d,", Q(i)); end do;

  printf("\n gaps: \n");

  for i to np do

     for j to nq do

      p := P(i); q := Q(j);

      g := gap(p, q);

      printf("%4d,", g);

  end do;

    printf("\n");

end do;

end proc

HostOfpq := proc(p, q)

local alpha, s, t, g, r, S, T, i, j;

   S := 1/2*q - 1/2;

   T := 1/2*p - 1/2;

   alpha := floor(q/p);

    r := q - alpha*p;

   for s to S do

     for t to T do

       g := abs((t*alpha - s)*p + t*r) - 1;

        printf("%4d,", g);

      end do;

     printf("\n");

 end do;

end proc

MapleSource.pdf

MapleSource.mw

(EDITED 2024/03/11  GMT 17H)

In a recent Question, @cq mentionned in its last reply "In fact, I wanted to do parameter sensitivity analysis and get the functional relationship between the [...] function and [parameters]. Later, i will study how the uncertainty of [the parameters] affects the [...] function".
I did not keep exchanging further on with @cq, simply replying that I could provide it more help if needed.

In a few words the initial problem was this one:

• Let X_1 and X_2 two random variables and G the random variable defined by  G = 1 - (X_1 - 1)^2/9 - (X_2 - 1)^3/16.
   
•  X_1 and X_2 are assumed to be gaussian random variables with respective mean and standard deviation equal to (theta_1, theta_3) and (theta_2, theta_4).
   
• The four theta parameters are themselves assumed to be realizations of four mutually independent uniform random variables Theta_1, ..., Theta_4 whose parameters are constants.
   
• Let QOI  (Quantity Of Interest) denote some scalar statistic of G (for instance its Mean, Variance, Skweness, ...).
  For instance, if QOI = Mean(G), then  QOI expresses itself as a function of the four parameters theta_1, ..., theta_4.
  The goal of @cq is to understand which of those parameters have the greatest influence on QOI.

For a quick survey of Sensitivity Analysys (SA) and a presentation of some of the most common strategies see Wiki-Overview


The simplest SA is the Local SA (LSA) we are all taught at school: having chosen some reference point P in the [theta_1, ..., theta_4] space the 1st order partial derivatives d[n] = diff(QOI, theta_n) expressed at point P give a "measure" (maybe after some normalization) of the sensitivity, at point P, of QOI regardibg each parameter theta_n.


A more interesting situation occurs when the parameters can take values in a neighorhood of  P which is not infinitesimal, or more generally in some domain without reference to any specific point P.
That is where Global SA (GSA) comes into the picture.
While the notion of local variation at some point P is well established, GSA raises the fundamental question of how to define how to measure the variation of a function over an arbitrary domain?
Let us take a very simple example while trying to answer this question "What is the variation of sin(x) over [0, 2*Pi]?"

1. If we focus on the global trend of sin(x)  mean there is no variation at all.
2. If we consider peak-to-peak amplitude the variation is equal to 2.
3. At last, if we consider L2 norm the variation is equal to Pi.
   (but the constant function x -> A/sqrt(2) has the same L2 norm but it is flat, and in some sense les fluctuating). 

Statisticians are accustomed to use the concept of variance as a measure to quantify the dispersion of a random variable. At the end of the sixties  one of them, Ilya Meyerovich Sobol’,  introduced the notion of Variance-Based GSA as the key tool to define the global variation of a function. This notion naturally led to that of Sobol' indices as a measure of the sensitivity of a funcion regarding one of its parameters or, which most important, regarding any combination (on says interaction) of its parameters.

The aim of this post is to show how Sobol' indices can be computed when the function under study has an analytic expression.
 

The Sobol' analysis is based on an additive decomposition of this function in terms of 2^P mutually orthogonal fiunctions where P is the number of its random parameters.
This decomposition and the ensuing integrations whose values will represent the Sobol' indices can be done analytically in some situation. When it is no longer the case specific numerical estimation methods have to be used;


The attached file contains a quite generic procedure to compute exact Sobol' indices and total Sobol' indices for a function whose parameters have any arbitrary statistical distribution.
Let's immediately put this into perspective by saying that these calculations are only possible if Maple is capable to find closed form expressions of some integrals, which is of course not always the case.

A few examples are also provided, including the one corresponding to @cq's original question.
At last two numerical estimation methods are presented.

SOBOL.mw

A checkered figure is a connected flat figure consisting of unit squares. The problem is to cut this figure into several equal parts (in area and shape). Cuts can only be made on the sides of the cells. In mathematics, such figures are called polyominoes, and the problem is called the tiling of a certain polyomino with copies of a single polyomino. See https://en.wikipedia.org/wiki/Polyomino

Below are 3 examples of such figures:

We will define such figures by the coordinates of the centers of the squares of which it consists. These points must lie in the first quarter, and points of this figure must lie on each of the coordinate axes.

Below are the codes for two procedures named  CutEqualParts  and  Picture . Required formal parameters of the first procedure: set  S  specifies the initial figure, r is the initial cell for generating subfigures, m is the number of parts into which the original figure needs to be divided. The optional parameter  s  equals (by default onesolution) or allsolutions. The starting cell  r  should be the corner cell of the figure. Then the set of possible subshapes for partitioning will be smaller. If there are no solutions, then the empty set will be returned. The second procedure  Picture  returns a picture of the obtained result as one partition (for a single solution) or in the form of a matrix if there are several solutions. In the second case, the optional parameter  d  specifies the number of rows and columns of this matrix.

restart;
CutEqualParts:=proc(S::set(list),r::list,m::posint, s:=onesolution)
local OneStep, n, i1, i2, j1, j2, R, v0, Tran, Rot, Ref, OneStep1, M, MM, MM1, T, T0, h, N, L;
n:=nops(S)/m;
if irem(nops(S), m)<>0 then error "Should be (nops(S)/m)::integer" fi;
if not (r in S) then error "Should be r in S" fi;
if m=1 then return {S} fi;
if m=nops(S) then return map(t->{t}, S) fi;
i1:=min(map(t->t[1],select(t->t[2]=0,S)));
i2:=max(map(t->t[1],select(t->t[2]=0,S)));
j1:=min(map(t->t[2],select(t->t[1]=0,S)));
j2:=max(map(t->t[2],select(t->t[1]=0,S)));
OneStep:=proc(R)
local n1, R1, P, NoHoles;
R1:=R;
n1:=nops(R1);
R1:={seq(seq(seq(`if`(r1 in S and not (r1 in R1[i]) , subsop(i={R1[i][],r1}, R1)[],NULL),r1=[[R1[i,j][1],R1[i,j][2]-1],[R1[i,j][1]+1,R1[i,j][2]],[R1[i,j][1],R1[i,j][2]+1],[R1[i,j][1]-1,R1[i,j][2]]]), j=1..nops(R1[i])), i=1..n1)};
NoHoles:=proc(s)
local m1, m2, M1, M2, M;
m1:=map(t->t[1],s)[1]; M1:=map(t->t[1],s)[-1];
m2:=map(t->t[2],s)[1]; M2:=map(t->t[2],s)[-1];
M:={seq(seq([i,j],i=m1..M1),j=m2..M2)}; 
if ormap(s1->not (s1 in s) and `and`(seq(s1+t in s, t=[[1,0],[-1,0],[0,1],[0,-1]])), M) then return false fi;
true;
end proc:
P:=proc(t)
if `and`(seq(seq(seq(([i,0] in t) and ([j,0] in t) and not ([k,0] in t) implies not ([k,0] in S), k=i+1..j-1), j=i+2..i2-1), i=i1..i2-2)) and `and`(seq(seq(seq(([0,i] in t) and ([0,j] in t) and not ([0,k] in t) implies not ([0,k] in S), k=i+1..j-1), j=i+2..j2-1), i=j1..j2-2))  then true else false fi;
end proc:
select(t->nops(t)=nops(R[1])+1 and NoHoles(t) and P(t) , R1);
end proc:
R:={{r}}:
R:=(OneStep@@(n-1))(R):
v0:=[floor(max(map(t->t[1], S))/2),floor(max(map(t->t[2], S))/2)]:
h:=max(v0);
Tran:=proc(L,v) L+v; end proc:
Rot:=proc(L, alpha,v0) <cos(alpha),-sin(alpha); sin(alpha),cos(alpha)>.convert(L-v0,Vector)+convert(v0,Vector); convert(%,list); end proc:
Ref:=proc(T) map(t->[t[2],t[1]], T); end proc:
OneStep1:=proc(T)
local T1, n2, R1;
T1:=T; n2:=nops(T1);
T1:={seq(seq(`if`(r1 intersect `union`(T1[i][])={}, subsop(i={T1[i][],r1}, T1), NULL)[], r1=MM1 minus T1[i]), i=1..n2)};
end proc:
N:=0; 
for M in R do
MM:={seq(seq(seq(map(t->Tran(Rot(t,Pi*k/2,v0),[i,j]),M),i=-h-1..h+1),j=-h-1..h+1),k=0..3),seq(seq(seq(map(t->Tran(Rot(t,Pi*k/2,v0),[i,j]),Ref(M)),i=-h-1..h+1),j=-h-1..h+1), k=0..3)}:
MM1:=select(t->(t intersect S)=t, MM):
T:={{M}}:
T:=(OneStep1@@(m-1))(T):
T0:=select(t->nops(t)=m, T):
if T0<>{} then if s=onesolution then return T0[1] else N:=N+1;
 L[N]:=T0[] fi; fi; 
od:
L:=convert(L,list);
if L[]::symbol then return {} else L fi;
end proc:

Picture:=proc(L::{list,set},Colors::list,d:=NULL)
local r;
uses plots, plottools;
if L::set or (L::list and nops(L)=1) or d=NULL then return
display( seq(polygon~(map(t->[[t[1]-1/2,t[2]-1/2],[t[1]+1/2,t[2]-1/2],[t[1]+1/2,t[2]+1/2],[t[1]-1/2,t[2]+1/2]] ,`if`(L::set,L[j],L[1][j])), color=Colors[j]),j=1..nops(Colors)) , scaling=constrained, size=[800,600]) fi;
if d::list then r:=irem(nops(L),d[2]);
if r=0 then return
display(Matrix(d[],[seq(display(seq(polygon~(map(t->[[t[1]-1/2,t[2]-1/2],[t[1]+1/2,t[2]-1/2],[t[1]+1/2,t[2]+1/2],[t[1]-1/2,t[2]+1/2]]  ,L[i,j]), color=Colors[j]),j=1..nops(Colors)), scaling=constrained, size=[400,300], axes=none), i=1..nops(L))])) else
display(Matrix(d[],[seq(display(seq(polygon~(map(t->[[t[1]-1/2,t[2]-1/2],[t[1]+1/2,t[2]-1/2],[t[1]+1/2,t[2]+1/2],[t[1]-1/2,t[2]+1/2]]  ,L[i,j]), color=Colors[j]),j=1..nops(Colors)), scaling=constrained, size=[400,300], axes=none), i=1..nops(L)), seq(plot([[0,0]], axes=none, size=[10,10]),k=1..d[2]-r)]))  fi; fi; 
end proc:

Examples of use for figures 1, 2, 3
In the first example for Fig.1 we get 4 solutions for m=4:

S:=({seq(seq([i,j], i=0..4), j=0..2)} union {[2,3],[3,3],[3,4]}) minus {[0,0],[0,1]}:
L:=CutEqualParts(S,[0,2],4,allsolutions);
C:=["Cyan","Red","Yellow","Green"]:
nops(L);
Picture(L,C,[2,2]);

In the second example for Fig.2 for m=2, we get 60 solutions (the first 16 are shown in the figure):

S:={seq(seq([i,j], i=0..4), j=0..4)} minus {[2,2]}:
L:=CutEqualParts(S,[0,0],2,allsolutions):
nops(L);
 C:=["Cyan","Red"]:
Picture(L[1..16],C,[4,4]);

In the third example for Fig.3 and  m=2  there will be a unique solution:

S:={seq(seq([i,j], i=0..5), j=0..3)}  minus {[5,0],[4,2]} union {[1,4],[2,4]}:
L:=CutEqualParts(S,[0,0],2):
 C:=["Cyan","Red"]:
Picture(L,C);

Addition. It is proven that the problem of tiling a certain polyomino with several copies of a single polyomino is NP-complete. Therefore, it is recommended to use the CutEqualParts procedure when the numbers  nops(S)  and  nops(S)/m  are relatively small (nops(S)<=24  and nops(S)/m<=12), otherwise the execution time may be unacceptably long.

Cutting_equal_parts.mw

A ball on a turntable can move in circles instead of falling off the edge (provided the initial conditions are appropriate). The effect was demonstrated in a video and can be simulated with MapleSim. The amination below shows a simulation of a frictionless case (falling off the table) and the case with a coefficient of friction of one.

Also demonstrated in the video: Tilting the table leads to a sideward (not a downhill) movement of the ball.

The presenter of the video noted that in the untilted state, the ball eventually drifts off the table, attributing this to slippage. This drift is also observable in the animation above, where the ball starts moving in a spiral, whereas in a Maple simulation (below to the left), the ball follows a perfect circle. Only after optimizing contact and initial conditions, MapleSim produced a result (using the same parameters) that exhibits a similar circle, with a slight difference in angular orientation after completing two revolutions about the center of the circle.

Some observations on the MapleSim model:

• The results only slightly depend on the solvers. Numerical inaccuracies do not seem to be the reason for the difference in orientation. (Edit: see update below for the reason).
• The ball bounces up and down in the MapleSim simulation (0.0025 of the balls radius). The bouncing is caused by the fact that the initial position of the ball is above the elastic equilibrium position with respect to the table (the elastic contact makes the ball sink in a bit). Adding damping in the settings of the contact element attenuates this effect and reduces the drift.
• Drift is not observable anymore if in the contact element setting for "kmu" (smoothness coefficient of sliding friction) is set to larger values (above 10 in this example). This is an indication that sliding friction occurs during the simulation.
• Choosing the equilibrium position as initial condition for the ball does not prevent initial bouncing because MapleSim sets an initial velocity for the ball that is directed away from the table. I did not manage to enforce strictly zero velocity. Maybe someone can tell why that is and how to set MapleSim to start the simulation without vertical velocity. (Edit: see update below for the reason).
• Assuming an initial velocity towards the contact to cancel the initial vertical velocity set by MapleSim substantially reduces bouncing and increases the diameter of the circle. This finally leads to a diameter that matches the Maple simulation. Therefore the initial bouncing combined with slippage seems to dissipate energy which leads to smaller circles. Depending on the contact settings and initial conditions for vertical movement the diameter of the circle varied moderately by about 15%.

In summary, MapleSim can be parametrized to simulate an ideal case without slippage. From there it should be possible to tune contact parameters to closely reproduce experiments where parameters are often not well known in advance.  

Some thoughts for future enhancement of MapleSim:

• In the model presented here, a patterned ball would have been helpful to visualize the tumbling movement of the ball. Marking two opposite sides of the ball with colored smaller spheres is a fiddly exercise that does not look nice.
• A sensor component that calculates the energy of a moving rigid body would have helped identifying energy loss of the system. Ideally the component could have two ports for the rotational and translational energy components. I see professional interest in such calculations and not only educational value for toy experiments.
• A port for slippage on the contact elements would have helped understanding where slippage occurs. Where slippage is, there is wear and this is also of interest for industrial applications.

Turntable_Paradox.msim (contains parameter sets for the above observations)

February 2nd is Groundhog Day in North America. A day when we look to small marmots to prognosticate the weather. If the groundhog sees its shadow when it emerges from its burrow, then it predicts 6 more weeks of winter, and if not, then spring begins today! Unfortunately there are many official weather predicting groundhogs. Fortunately, the excellent website Groundhog-Day.com tracks each of their predictions. Unfortunately, it doesn't tell you which groundhog to trust. Fortunately, it has an API and we can take the data and map it in Maple:

This map assumes that each groundhog's prediction is valid at it's exact geographic coordinates, but that it's predictive powers fall off in inverse proportion to the distance away.  So, exactly halfway between a groundhog predicting early spring, and one predicting 6 more weeks of winter, we expect 3 more weeks of winter.  I handle that with Maple's Interpolation:-InverseDistanceWeightedInterpolation command with a radius of 1500 miles.  I plot a contourplot of that interpolating function, and then display it with the world map in DataSets:-Builtin:-WorldMap to generate the image above.

All the code to do that can be found in the following worksheet which also using the URL packaget to fetch the most recent groundhog data possible from the website.

Groundhog-Map.mw

I've commented out a few lines that you might use to explore other possible maps.  You can filter to file to just include real living groundhogs and not all the other precitions (some from puppets, some from other animals) if you find that more trust worthy. You might also prefer to change the interpolation command, one of my collegues suggests that Interpolation:-NaturalNeighborInterpolation might be a better choice.

[Right-click on image and open in new tab to see larger]

To Scan Math with the Maple Calculator and show solution steps in Maple:

1. first scan some math with the calculator

2. Maple calculator immediately shows the solution if that is what you are looking for:

3. Calculator gives options to show the solution steps in the calculator itself ( footprint button in top-right) 

4. Or to upload the math to the MapleCloud (cloud icon with up arrow)

5. Once the math is uploaded, MapleCloud can be loaded on a desktop computer and the file opened from your account's Maple Calculator group of files:

6. Again, the solution and some more details are visible on Maple Cloud:

7. To open this math in Maple, click the blue button to Download the file.

The downloaded file can then be loaded in Maple:

8. The Maple commands to solve this math are shown, and the result. 

To show steps in Maple at this point, convert the math to inert form, then run the Student:-Calculus1:-ShowSolution() command on it:

Download MapleCalculatorMathCloudUpload.mw

I'm an engineer and when showing results of calculations, some values will display as fractions, and I would prefer that instead floating numbers are displayed.  Also, there is kind of a quirk where if the multiplier of a unit is 1, the result displays as a unit only. I would prefer to see 1*A rather than A.

I wrote this simple proc to convert a value with or without attached unit to a floating point number if it is a fraction or if it has a unit and the coefficient is 1.

Let me know if there is a more elegant way to do this or you have any suggestions or questions.

   unrationalize := proc(x)
        local 
            returnval,
            localcoeff,
            localunit
        ;
        description
            "Converts a fractional number to float",
            "Units are supported"
        ;
        if type(x, fraction) or type(x, with_unit(fraction)) then 
            returnval := evalf(x)
        elif type(x, with_unit(1, 'localcoeff', 'localunit')) then
            returnval :=  evalf(localcoeff)*localunit
        else 
            returnval := x;
        end if;
        return returnval;
    end  proc;
# Testing the proc
list1 := [1/2, 1/2*Unit(('A')/('V')), 1, Unit('A')];
listDescription := ["Fraction", "Fraction with Unit", "Unity", "Unity with Unit"];
for i, myValue in list1 do
    [listDescription[i], "evalf:", evalf(myValue), "unrationalize:", unrationalize(myValue)];
end do;

Curve sketching is an important skill for all calculus students to learn. In an era where technology is increasingly relied upon to perform mathematical computations and representations, maintaining fundamental skills such as curve sketching is more important than ever.

The new “Curve Sketching” collection is now available on Maple Learn. This collection provides background information on the process of curve sketching and opportunities to put this knowledge into practice. By starting with the “Curve Sketching Guide” and “Relationships Between Derivatives” documents, students are exposed to observational and computational strategies for drawing a function and its 1st and 2nd derivatives.

After looking through these documents, students can begin to practice sketching by observing and interpreting graphical properties with the “Sketch Derivative From Function Graph”, “Sketch Second Derivative From Function Graph”, and “Sketch Function From Derivative Graphs” activities:

Once a student has mastered extracting sketching information by graphical observation, they are ready for the next step: extracting information from a function’s definition. At this point, the student is ready to try sketching from a blank canvas with the “Sketch Curve From Function Definition” activity:

This collection also has activities for students below the calculus level. For example, the “Curve Sketching Quadratics Activity”, can be completed using only factoring strategies:

Whether you are a quadratics rookie or a calculus pro, this collection has an interactive activity to challenge your knowledge. Have fun sketching!

My friend and colleague Nic Fillion and I are writing another book, this one on perturbation methods using backward error analysis (and Maple).  We have decided to make the supporting materials available by means of Jupyter notebooks with a Maple kernel (there are some Maple worksheets and workbooks already, but going forward we will use Jupyter).

The presentation style is meant to aid reproducibility, and to allow others to solve related problems by changing the scripts as needed.

The first one is up at 

https://github.com/rcorless/Perturbation-Methods-in-Maple

Comments very welcome.  This particular method is a bit advanced in theory (but it's very simple in practice, for weakly nonlinear oscillators).  I haven't coded for efficiency and there may be some improvements possible ("may" he says, sheesh).  Comments on that are also welcome.

-r

In the most recent issue of Maple Transactions, I published (with David Jeffrey, and with a student named Johan Joby) a paper that used Jupyter Notebook with a Maple kernel as the main vehicle.  Have a look, and let me know what you think.

Two-cycles in the infinite exponential tower

Download bugs-irrationals.mw

Why this post
This work was intended to be a simple reply to a question asked a few days ago.
At some point, I realised that the approach I was using could have a more general interest which, in my opinion, was worth a post.
In a few words, this post is about solving an algebra problem using a method originally designed to tackle statistical problems.

The Context
Recently @raj2018 submitted a question I'm going to resume this way:

Let S(phi ;  beta, f) a function of phi parameterized by beta and f.
Here is the graph of S(phi ;  0.449, 0.19)  @raj2018 provided

@raj2018 then asked how we can find other values (A, B)  of values for (beta, f) such that the graph of S(phi, A, B) has the same aspect of the graph above.
More precisely, let phi_0 the largest strictly negative value of phi such that  S(phi_0, A, B) = 0.
Then  S(phi, A, B) must be negative (strictly negative?) in the open interval (phi_0, 0), and must have exactly 3 extrema in this range.
I will said the point  (A, B) is admissible is S(phi, A, B) verifies thess conditions

The expression of S(phi, A, B) is that complex that it is likely impossible to find an (several?, all?) admissible point using analytic developments.

The approach

When I began thinking to this problem I early thought to find the entire domain of admissible points: was it something possible, at least with some reasonable accuracy? 

Quite rapidly I draw an analogy with an other type of problems whose solution is part of my job: the approximate construction of the probability density function (PDF) of multivariate random variables (obviously this implies that no analytical expression of this PDF is available). This is a very classical problem in Bayesian Statistics, for instance when we have to construt an approximation of a posterior PDF.

To stick with this example and put aside the rare situations where this PDF can be derived analytically, building a posterior PDF is largely based on specific numerical methods. 
The iconic one is known under the generic name MCMC  which stands for Markov Chain Monte Carlo.

Why am I speaking about MCMC or PDF or even random variables?
Let us consider some multivariate random variable R whose PDF as a constant on some bounded domain D and is equal to 0 elsewhere. R is then a uniform random variable with support supp(R) = D.
Assuming the domain Adm of admissible (beta, f) is bounded, we may  think of it as the support of some uniform random variable. Following this analogy we may expect to use some MCMC method to "build the PDF of the bivariate random variable (beta, f)", otherwise stated "to capture​​​​​​ the boundary of​ Adm".

The simplest MCMC method is the Metropolis-Hastings algorithm (MH).
In a few simple words MH builds a Markov chain this way:

1. Let us assume that the chain already contains elements e₁, ..., e_n.
   Let  f  some suitable "fitness" function (whose nature is of no importance right now).
2. A potential new element c is randomly picked in some neighborhood or e_n.
3. If the ratio f (c) / f (e_n) is larger than 1, we decide to put c into the chain (thus e_n+1 = c) otherwise we leave it to chance to decide whether or not c iis put into the chain.
   If chance decides the contrary,  then e_n is duclicated (thus e_n+1 = e_n).

MH is not the most efficient MCMC algorithm but it is efficient enough for what we want to achieve.
The main difficulty here is that there is no natural way to build the fitness function  f , mainly because the equivalent random variable I talked about is a purely abstract construction.

A preliminary observation is that if S(phi, beta, f) < 0 whatever phi in (phi_0, 0), then S has an odd number of extrema in (phi_0, 0). The simplest way to find these extrema is to search for the zeros of the derivative S' of S with respect to phi, while discardinq those where the second derivative can reveal "false" extrema where both S'' of S is null (I emphasize this last point because I didn't account for it in attached file).
The algorithm designed in this file probably misses a few points for not checking if S''=0, but it is important to keep in mind that we don't want a complete identification of  Adm but just the capture of its boundary.
Unless we are extremely unlucky there is only a very small chance that omitting to check if S''=0 will deeply modify this boundary.

How to define function f  ?
What we want is that  f (c) / f (e_n) represents the probability to decide wether c is an admissible point or not. In a Markov chain this  ratio represents how better or worse c is relatively to e_n, and this is essential for the chain to be a true Markov chain.
But as our aim is not to build a true Markov chain but simply a chain which looks like a Markov chain, we we can take some liberties and replace  f (c) / f (e_n) by some function  g(c) which quantifies the propability for c to be an admissible couple. So we want that  g(c) = 1 if  S(phi, c) has exactly M=3 negative extrema and  g(c) < 1 if M <> 3.
The previous algorihm transforms into:

1. Let us assume that the chain already contains elements e₁, ..., e_n.
   Let  g  a function which the propability that element is admissible
2. A potential new element c is randomly picked in some neighborhood or e_n.
3. If the ratio g(c) is larger than 1, we decide to put c into the chain (thus e_n+1 = c) otherwise we leave it to chance to decide whether or not c iis put into the chain.
   If chance decides the contrary,  then e_n is duclicated (thus e_n+1 = e_n).

This algorithm can also be seen as a kind of genetic algorithm.

A possible choice is  g(c)= exp(-|3-M|).
In the attached file I use instead the expression g(c) = (M + 1) / 4 fo several reasons:

• It is less sharp at M=3 and thus enables more often to put c into the chain, which increases its exploratory capabilities.
• The case M > 3, which no preliminary investigation was able to uncover, is by construction eliminated in the procedure Extrema which use an early stopping strategy (if as soon as more than M=3 negative extrema are found the procedure stops).

The algorithm I designed basically relies upon two stages:

1. The first one is aimed to construct a "long" Markov-like chain ("long" and not long because Markov chains are usually much longer than those I use).
   There are two goals here:
   1. Check if Adm is or not simply-connected or not (if it has holes or not).
   2. Find a first set of admissible points that can be used as starting points for subsequent chains.
       
2. Run several independent Markov-like chains from a reduced set of admissible points.
   The way this reduced set is constructed depends on the goal to achieve:
   1. One may think of adding points among those already known in order to assess the connectivity of Adm,
   2. or refinining the boundary of Adm.

These two concurent objectives are mixed in an ad hoc way depending on the observation of the results already in hand.

We point here an important feature of MCMC methods: behind their apparent algorithmic simplicity, it is common that high-quality results can only be obtained efficiently at the cost of problem-dependent tuning.

A last word to say that after several trials and failures I found it simpler to reparameterize the problems in terms of (phi_0, f) instead of (beta, f).

Codes and results

Choice g(c) = (M + 1) / 4 
The code : Extrema_and_MCMC.mw

To access the full results I got load this m file (do not bother its extension, Mapleprimes doesn't enable uploading m files) MCMC_20231209_160249.mw (save it and change it's extension in to m instead mw)

EDITED: choice  g(c)= exp(-|3-M|)
Here are the files contzining the code and the results:
Extrema_and_MCMC_g2.mw
MCMC_20231211_084053.mw

To ease the comparison of the two sets of results I used the same random seeds inn both codes.
Comparing the results got around the first admissible point is straightforward.
It's more complex for @raj2018's solution because the first step of the algorithim (drawing of a sibgle chain of length 1000) finds six times more admissible point with g(c)= exp(-|3-M|) than with g(c) = (M + 1) / 4.