Hoping to get your very own copy of Maple for Christmas? 

No, I'm not about to announce a sale. No crass commercialism allowed on MaplePrimes!  But we have created this handy-dandy letter for Santa that students can use to try to convince him to put Maple under the tree this year*.

Happy holidays from Maplesoft!

*Results not guaranteed, but we hope you will agree we made a valiant attempt. :-)

On the convex part of the surface we place a curve (not necessarily flat, as in this case). We divide this curve into segments of equal length (in the text Ls [i]) and divide the path that our surface will roll (in the text L [i]) into segments of the same length as segments of curve. Take the next segment of the trajectory L [i] and the corresponding segment on the curve Ls [i], calculate the angles between them. After that, we perform well-known transformations that place the curve in the space so that the segment Ls [i] coincides with the segment L [i]. At the same time, we perform exactly the same transformations with the equation of surface.

For example, the ellipsoid rolls on the oX1 axis, and each position of the ellipsoid in space corresponds to the equation in the figure.
Rolling_of_surface.mw

And similar examples:

Download PDE_and_BC_during_2018.mw

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

Error when running 'm_model.m' to generate simulink function block when using 'Group all parameters into nested structure'. No errors when using 'Do not group parameters'.

Error using m_fullSwing (line 125)
Invalid setting in 'MapleSim_fullSwing/MapleSim_fullSwing/MapleSimParameters' for parameter 'Value'.
Caused by:
    Error using m_fullSwing (line 125)
    Expression '[Param.Address.ArmPlane, Param.Address.ClubRelHand,  ... ]' for
    parameter 'Value' in 'MapleSim_fullSwing/MapleSim_fullSwing/MapleSimParameters' cannot be evaluated.
        Error using m_fullSwing (line 125)
        Error: Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for
        mismatched delimiters.
 

Yesterday, I accidentally discovered a nasty bug in a fairly simple example:

restart;
Expr:=a*sin(x)+b*cos(x);
maximize(Expr, x=0..2*Pi);
minimize(Expr, x=0..2*Pi);
                                    

I am sure the correct answers are  sqrt(a^2+b^2)  and  -sqrt(a^2+b^2)  for any real values  a  and  b .  It is easy to prove in many ways. The simplest method does not require any calculations and can be done in the mind. We will consider  Expr  as the scalar product (or the dot product) of two vectors  <a, b>  and  <sin(x), cos(x)>, one of which is a unit vector. Then it is obvious that the maximum of this scalar product is reached if the vectors are codirectional and equals to the length of the first vector, that is, sqrt(a^2+b^2).

Bugs in these commands were noted by users and earlier (see search by keywords bug, maximize, minimize) but unfortunately are still not fixed. 

The attached worksheet develops a procedure for extrapolating boundary data from a square to its interior.  Specifically, let's consider the square [−1,1]×[−1,1] and the continuous functions u₁, u₂, u₃, u₄ defined on its edges. The procedure constructs a continuous function u(x,y) in the interior of the square which matches the boundary data.

The function u(x,y) is necessarily discontinuous at a corner if the boundary data of the edges meeting at that corner are inconsistent.  However, if the boundary data are consistent at all corners, then u(x,y) is continuous everywhere including the boundary.

Here is what the extrapolating function looks like:

Download putnam2018.mw

Edit.
Maple can be also very useful in solving the difficult problem B6; it can be reduced to compute the (huge) coefficient of x^1842 in the polynomial

g := (1 + x + x^2 + x^3 + x^4 + x^5 + x^9)^2018;

The computation is very fast:

coeff(g, x, 1842):   evalf(%);
         0.8048091229e1147

While generating a 3D plot of the solution of an ODE with a parameter, I noticed that better performance could be obtained by calling the plot3d command with a procedure argument, done in a special manner.

I don't recall this being discussed before, so I'll share it. (It it has been discussed, and this is widely known, then I apologize.)

I tweaked the initial conditions of the original ODE system, to obtain a non-trivial solution. I don't think that the particular nature of the solution has a bearing on this note.

Download ode_param_plot.mw

We have just released updates to Maple and MapleSim, which provide corrections and improvements to product functionality.

As usual, the Maple update is available through Tools>Check for Updates in Maple, and is also available from our website on the Maple 2018.2.1 download page, where you can also find more details.

For MapleSim users, use Help>Check for Updates or visit the MapleSim 2018.2.1 download page.

Just a simple use of DataFrames.  Looking for a new flashlight, I decided to compile a small selection of flashlights for comparison.
 

Download flashlight3.mw

From time to time, people ask me about visualizing knots in Maple. There's no formal "Knot Theory" package in Maple per se, but it is certainly possible to generate many different knots using a couple of simple commands. The following shows various examples of knots visualized using the plots:-tubeplot and algcurves:-plot_knot commands.

The unknot can be defined by the following parametric equations:

x=sin(t)

y=cos(t)

z=0

The trefoil knot can be defined by the following parametric equations:

x = sin(t) + 2*sin(2*t)

y = cos(t) + 2*sin(2*t)

z = sin(3*t)

The figure-eight can be defined by the following parametric equations:

x = (2 + cos(2*t)) * cos(3*t)

y = (2 + cos(2*t)) * sin(3*t)

z = sin(4*t)

The Lissajous knot can be defined by the following parametric equations:

Where n[x], n[y], and n[z] are integers and the phase shifts phi[x], phi[y], and phi[z] are any real numbers.
The 8 21 knot ( n[x] = 3, n[y] = 4, and n[z] = 7) appears as follows:
 

A star knot can be defined by using the following polynomial:
 

By switching x and y, different visualizations can be generated:

Please feel free to add more of your favourite knot visualizations in the comments below!

You can interact with the examples or download a copy of these examples from the MapleCloud here: https://maple.cloud/app/5654426890010624/Examples+of+Knots

In this post an interesting geometric problem is solved: for an arbitrary convex polygon, find a straight line that divides both the area and the perimeter in half. The following results on this problem are well known:
1. For any convex polygon there is such a straight line.
2. For any convex polygon into which a circle can be inscribed, in particular for any triangle, the desired line must pass through the center of the inscribed circle.
3. For a triangle, the number of solutions can be 1, 2, or 3.
4. If the polygon has symmetry with respect to a point, then any straight line passing through this point is a solution.

The procedure called  InHalf  (the code below) symbolically solves this problem. The formal parameter of the procedure is the list of coordinates of the vertices of a convex polygon (vertices must be passed opposite or clockwise). The procedure returns all solutions in the form of a list of pairs of points, lying on the perimeter of the polygon, that are the ends of segments that implement the desired dividing.

Download In_Half.mw

We have just released an update to Maple, Maple 2018.2. This release includes improvements in a variety of areas, including code edit regions, Workbooks, and Physics, as well as support for macOS 10.14.

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

For MapleSim users, the update includes optimizations for handling large models, improvements to model import and export, updates to the hydraulics and pneumatics libraries, and more. For more details and download instructions, visit the MapleSim 2018.2 download page.

The procedure presented here does independence tests of a contingency table by four methods:

1. Pearson's chi-squared (equivalent to Statistics:-ChiSquareIndependenceTest),
2. Yates's continuity correction to Pearson's,
3. G-chi-squared,
4. Fisher's exact.

(All of these have Wikipedia pages. Links are given in the code below.) All computations are done in exact arithmetic. The coup de grace is Fisher's. The first three tests are relatively easy computations and give approximations to the p-value (the probability that the categories are independent), but Fisher's exact test, as its name says, computes it exactly. This requires the generation of all matrices of nonnegative integers that have the same row and column sums as the input matrix, and for each of these matrices computing the product of the factorials of its entries. So, there are relatively few implementations of it, and perhaps none that do it exactly. (Could some with access check Mathematica please?)

Our own Joe Riel's amazing and fast Iterator package makes this computation considerably easier and faster than it otherwise would've been, and I also found inspiration in his example of recursively counting contingency tables found at ?Iterator,BoundedComposition. 

ContingencyTableTests:= proc(
   O::Matrix(nonnegint), #contingency table of observed counts 
   {method::identical(Pearson, Yates, G, Fisher):= 'Pearson'}
)
description 
   "Returns p-value for Pearson's (w/ or w/o Yates's continuity correction)" 
   " or G chi-squared or Fisher's exact test."
   " All computations are done in exact arithmetic."
;
option
   author= "Carl Love <carl.j.love@gmail.com>, 27-Oct-2018",
   references= (                                                           #Ref #s:
      "https://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test",         #*1
      "https://en.wikipedia.org/wiki/Yates%27s_correction_for_continuity",  #*2
      "https://en.wikipedia.org/wiki/G-test",                               #*3
      "https://en.wikipedia.org/wiki/Fisher%27s_exact_test",                #*4
      "Eric W Weisstein \"Fisher's Exact Test\" _MathWorld_--A Wolfram web resource:"
      " http://mathworld.wolfram.com/FishersExactTest.html"                 #*5
   )
;
uses AT= ArrayTools, St= Statistics, It= Iterator;
local
   #column & row sums: 
   C:= AT:-AddAlongDimension(O,1), R:= AT:-AddAlongDimension(O,2),
   r:= numelems(R), c:= numelems(C), #counts of rows & columns
   n:= add(R), #number of observations
   #matrix of expected values under null hypothesis (independence):
   #(A 0 entry would mean a 0 row or column in the original, which is not allowed.)
   E:= Matrix((r,c), (i,j)-> R[i]*C[j], datatype= 'positive') / n,
   #Pearson's, Yates's, and G all use a chi-sq statistic, each computed by 
   #slightly different formulae.
   Chi2:= add@~table([
       'Pearson'= (O-> (O-E)^~2 /~ E),                     #see *1
       'Yates'= (O-> (abs~(O - E) -~ 1/2)^~2 /~ E),        #see *2
       'G'= (O-> 2*O*~map(x-> `if`(x=0, 0, ln(x)), O/~E))  #see *3
   ]), 
   row, #alternative rows generated for Fisher's
   Cutoff:= mul(O!~), #cut-off value for less likely matrices
   #Generate recursively all contingency tables whose row and column sums match O.
   #Compute their probabilities under independence. Sum probabilities of all those
   #at most as probable as O. (see *5, *4)
   #Parameters: 
   #   C = column sums remaining to be filled; 
   #   F = product of factorials of entries of contingency table being built;
   #   i = row to be chosen this iteration
   AllCTs:= (C, F, i)->
      if i = r then #Recursion ends; last row is C, the unused portion of column sums. 
         (P-> `if`(P >= Cutoff, 1/P, 0))(F*mul(C!~))
      else
         add(
            thisproc(C - row[], F*mul(row[]!~), i+1), 
            row= It:-BoundedComposition(C, R[i])
         )
      fi      
;
   userinfo(1, ContingencyTableTests, "Table of expected values:", print(E));
   if method = 'Fisher' then AllCTs(C, 1, 1)*mul(R!~)*mul(C!~)/n!
   else 1 - St:-CDF(ChiSquare((r-1)*(c-1)), Chi2[method](O)) 
   fi   
end proc:

The worksheet below contains the code above and one problem solved by the 4 methods

Download FishersExact.mw

Problem:

Suppose you have a bunch of 2D data points which:

1. May include points with the same x-value but different y-values; and
2. May be unevenly-spaced with respect to the x-values.

How do you clean up the data so that, for instance, you are free to construct a connected data plot, or perform a Discrete Fourier Transform? Please note that Curve Fitting and the Lomb–Scargle Method, respectively, are more effective techniques for these particular applications. Let's start with a simple example for illustration. Consider this Matrix:

A := < 2, 5; 5, 8; 2, 1; 7, 8; 10, 10; 5, 7 >;

Consolidate:

First, sort the rows of the Matrix by the first column, and extract the sorted columns separately:

P := sort( A[..,1], output=permutation ); # permutation to sort rows by the values in the first column
U := A[P,1]; # sorted column 1
V := A[P,2]; # sorted column 2

We can regard the sorted bunches of distinct values in U as a step in a stair case, and the goal is replace each step with the average of the y-values in V located on each step.

Second, determine the indices for the first occurrences of values in U, by selecting the indices which give a jump in x-value:

m := numelems( U );
K := [ 1, op( select( i -> ( U[i] > U[i-1] ), [ seq( j, j=2..m ) ] ) ), m+1 ];
n := numelems( K );

The element m+1 is appended for later convenience. Here, we can quickly define the first column of the consolidated Matrix:

X1 := U[K[1..-2]];

Finally, to define the second column of the consolidated Matrix, we take the average of the values in each step, using the indices in K to tell us the ranges of values to consider:

Y1 := Vector[column]( n-1, i -> add( V[ K[i]..K[i+1]-1 ] ) / ( K[i+1] - K[i] ) );

Thus, the consolidated Matrix is given by:

B := < X1 | Y1 >;

Spread Evenly:

To spread-out the x-values, we can use a sequence with fixed step size:

X2 := evalf( Vector[column]( [ seq( X1[1]..X1[-1], (X1[-1]-X1[1])/(m-1) ) ] ) );

For the y-values, we will interpolate:

Y2 := CurveFitting:-ArrayInterpolation( X1, Y1, X2, method=linear );

This gives us a new Matrix, which has both evenly-spaced x-values and consolidated y-values:

C := < X2 | Y2 >;

Plot:

plots:-display( Array( [
        plots:-pointplot( A, view=[0..10,0..10], color=green, symbol=solidcircle, symbolsize=15, title="Original Data", font=[Verdana,15] ),
        plots:-pointplot( B, view=[0..10,0..10], color=red, symbol=solidcircle, symbolsize=15, title="Consolidated Data", font=[Verdana,15] ),
        plots:-pointplot( C, view=[0..10,0..10], color=blue, symbol=solidcircle, symbolsize=15, title="Spread-Out Data", font=[Verdana,15] )
] ) );

Sample Data with Noise:

For another example, let’s take data points from a logistic curve, and add some noise:

# Noise generators
f := 0.5 * rand( -1..1 ):
g := ( 100 - rand( -15..15 ) ) / 100:

# Actual x-values
X := [ seq( i/2, i=1..20 ) ];

# Actual y-values
Y := evalf( map( x -> 4 / ( 1 + 3 * exp(-x) ), X ) );

# Matrix of points with noise
A := Matrix( zip( (x,y) -> [x,y], map( x -> x + f(), X ), map( y -> g() * y, Y ) ) );

Using the method outlined above, and the general procedures defined below, define:

B := ConsolidatedMatrix( A );
C := EquallySpaced( B, 21, method=linear );

Visually:

plots:-display( Array( [
    plots:-pointplot( A, view=[0..10,0..5], symbol=solidcircle, symbolsize=15, color=green, title="Original Data", font=[Verdana,15] ),
    plots:-pointplot( B, view=[0..10,0..5], symbol=solidcircle, symbolsize=15, color=red, title="Consolidated Data", font=[Verdana,15]  ),
    plots:-pointplot( C, view=[0..10,0..5], symbol=solidcircle, symbolsize=15, color=blue, title="Spread-Out Data", font=[Verdana,15] )
] ) );

Generalization:

Below are more generalized custom procedures, which are used in the above example. These also account for special cases.

# Takes a matrix with two columns, and returns a new matrix where the new x-values are unique and sorted,
# and each new y-value is the average of the old y-values corresponding to the x-value.
ConsolidatedMatrix := proc( A :: 'Matrix'(..,2), $ )

        local i, j, K, m, n, P, r, U, V, X, Y:
  
        # The number of rows in the original matrix.
        r := LinearAlgebra:-RowDimension( A ):

        # Return the original matrix should it only have one row.
        if r = 1 then
               return A:
        end if:

        # Permutation to sort first column of A.
        P := sort( A[..,1], ':-output'=permutation ):       

        # Sorted first column of A.
        U := A[P,1]:

        # Corresponding new second column of A.
        V := A[P,2]:

        # Return the sorted matrix should all the x-values be distinct.
        if numelems( convert( U, ':-set' ) ) = r then
               return < U | V >:
        end if:

        # Indices of first occurrences for values in U. The element m+1 is appended for convenience.
        m := numelems( U ):
        K := [ 1, op( select( i -> ( U[i] > U[i-1] ), [ seq( j, j=2..m ) ] ) ), m+1 ]:
        n := numelems( K ):

        # Consolidated first column.
        X := U[K[1..-2]]:

        # Determine the consolidated second column, using the average y-value.
        Y := Vector[':-column']( n-1, i -> add( V[ K[i]..K[i+1]-1 ] ) / ( K[i+1] - K[i] ) ):

        return < X | Y >:

end proc:

# Procedure which takes a matrix with two columns, and returns a new matrix of specified number of rows
# with equally-spaced x-values, and interpolated y-values.
# It accepts options that can be passed to ArrayInterpolation().
EquallySpaced := proc( M :: 'Matrix'(..,2), m :: posint )

        local A, i, r, U, V, X, Y:

        # Consolidated matrix, the corresponding number of rows, and the columns.
        A := ConsolidatedMatrix( M ):
        r := LinearAlgebra:-RowDimension( A ):
        U, V := evalf( A[..,1] ), evalf( A[..,2] ):

        # If the consolidated matrix has only one row, return it.
        if r = 1 then
               return A:
        end if:

        # If m = 1, i.e. only one equally-spaced point is requested, then return a matrix of the averages.
        if m = 1 then
               return 1/r * Matrix( [ [ add( U ), add( V ) ] ] ):
        end if:

        # Equally-spaced x-values.
        X := Vector[':-column']( [ seq( U[1]..U[-1], (U[-1]-U[1])/(m-1), i=1..m ) ] ):

        # Interpolated y-values.
        Y := CurveFitting:-ArrayInterpolation( U, V, X, _rest ):    

        return < X | Y >:

end proc:

Worth Checking Out:

• The CurveFitting:-ArrayInterpolation() command
• The dataplot() command
• The plots:-pointplot() command
• The Statistics:-ScatterPlot() command
• The Interpolation package (new in Maple 2018)