Applications, Examples and Libraries

Share your work here

 

Minimize the number of tensor components according to its symmetries
(and relabel, redefine or count the number of independent tensor components)

 

 

The nice development described below is work in collaboration with Pascal Szriftgiser from Laboratoire PhLAM, Université Lille 1, France, used in the Mapleprimes post Magnetic traps in cold-atom physics

 

A new keyword in Define  and Setup : minimizetensorcomponents, allows for automatically minimizing the number of tensor components taking into account the tensor symmetries. For example, if a tensor with two indices in a 4D spacetime is defined as antisymmetric using Define with this new keyword, the number of different tensor components will be exactly 6, and the elements of the diagonal are automatically set equal to 0. After setting this keyword to true with Setup , all subsequent definitions of tensors automatically minimize the number of components while using this keyword with Define  makes this minimization only happen with the tensors being defined in the call to Define .

 

Related to this new functionality, 4 new Library routines were added: MinimizeTensorComponents, NumberOfIndependentTensorComponents, RelabelTensorComponents and RedefineTensorComponents

 

Example:

restart; with(Physics)

 

Define an antisymmetric tensor with two indices

Define(F[mu, nu], antisymmetric)

`Defined objects with tensor properties`

 

{Physics:-Dgamma[mu], F[mu, nu], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu]}

(1.1)

Although the system knows that F[mu, nu] is antisymmetric, you need to use Simplify to apply the (anti)symmetry

F[mu, nu]+F[nu, mu]

F[mu, nu]+F[nu, mu]

(1.2)

 

Simplify(F[mu, nu]+F[nu, mu])

0

(1.3)

so by default the components of F[mu, nu] do not automatically reflect the (anti)symmetry; likewise

F[1, 2]+F[2, 1]

F[1, 2]+F[2, 1]

(1.4)

Simplify(F[1, 2]+F[2, 1])

0

(1.5)

and computing the array form of F[mu, nu]we do not see the elements of the diagonal equal to zero nor the lower-left triangle equal to the upper-right triangle but for a different sign:

TensorArray(F[mu, nu])

Matrix(%id = 18446744078270093062)

(1.6)

 

On the other hand, this new functionality, here called minimizetensorcomponents, makes the symmetries of the tensor be explicitly reflected in its components.

 

There are three ways to use it. First, one can minimize the number of tensor components of a tensor previously defined. For example

 

Library:-MinimizeTensorComponents(F)

Matrix(%id = 18446744078270064630)

(1.7)

After this, both (1.2) and (1.3) are automatically equal to 0 without having to use Simplify

F[mu, nu]+F[nu, mu]

0

(1.8)

0

0

(1.9)

And the output of TensorArray  in (1.6) becomes equal to (1.7).

 

NOTE: in addition, after using minimizetensorcomponents in the definition of a tensor, say F, all the keywords implemented for Physics tensors are available for F:

 

F[]

F[mu, nu] = Matrix(%id = 18446744078247910206)

(1.10)

F[trace]

0

(1.11)

F[nonzero]

F[mu, nu] = {(1, 2) = F[1, 2], (1, 3) = F[1, 3], (1, 4) = F[1, 4], (2, 1) = -F[1, 2], (2, 3) = F[2, 3], (2, 4) = F[2, 4], (3, 1) = -F[1, 3], (3, 2) = -F[2, 3], (3, 4) = F[3, 4], (4, 1) = -F[1, 4], (4, 2) = -F[2, 4], (4, 3) = -F[3, 4]}

(1.12)

"F[~1,mu,matrix]"

F[`~1`, mu] = Vector[row](%id = 18446744078247885990)

(1.13)

Alternatively, one can define a tensor, specifying that the symmetries should be taken into account to minimize the number of its components passing the keyword minimizetensorcomponents to Define .

 

Example:

 

Define a tensor with the symmetries of the Riemann  tensor, that is, a tensor of 4 indices that is symmetric with respect to interchanging the positions of the 1st and 2nd pair of indices and antisymmetric with respect to interchanging the position of its 1st and 2nd indices, or 3rd and 4th indices, and define it minimizing the number of tensor components

 

Define(R[alpha, beta, mu, nu], symmetric = {[[1, 2], [3, 4]]}, antisymmetric = {[1, 2], [3, 4]}, minimizetensorcomponents)

`Defined objects with tensor properties`

 

{Physics:-Dgamma[mu], F[mu, nu], Physics:-Psigma[mu], R[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu]}

(1.14)

We now have

R[1, 2, 3, 4]+R[2, 1, 3, 4]

0

(1.15)

R[alpha, beta, mu, nu]-R[mu, nu, alpha, beta]

0

(1.16)
• 

One can always retrieve the symmetry properties in the abstract notation used by the Define command using the new Library:-GetTensorSymmetryProperties, its output is ordered, first the symmetric then the antisymmetric properties

 

Library:-GetTensorSymmetryProperties(R)

{[[1, 2], [3, 4]]}, {[1, 2], [3, 4]}

(1.17)
• 

After making the symmetries explicit (and also before that), it is frequently useful to know the number of independent components of a given tensor. For this purpose you can use the new Library:-NumberOfIndependentTensorComponents

 

Library:-NumberOfIndependentTensorComponents(R)

21

(1.18)

and besides taking into account the symmetries, in the case of the Riemann  tensor, after taking into account the first Bianchi identity this number of components is further reduced to 20.

 

A third way of using the new minimizetensorcomponents functionality is using Setup , so that, automatically, every subsequent definition of tensors with symmetries is performed minimizing the number of its components using the indicated symmetries

 

Example:

Setup(minimizetensorcomponents = true)

[minimizetensorcomponents = true]

(1.19)

So from hereafter you can define tensors taking into account their symmetries explicitly and without having to include the keyword minimizetensorcomponents at each definition

 

Define(C[alpha, beta], antisymmetric)

`Defined objects with tensor properties`

 

{C[mu, nu], Physics:-Dgamma[mu], F[mu, nu], Physics:-Psigma[mu], R[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu]}

(1.20)

 

C[]

C[mu, nu] = Matrix(%id = 18446744078408747598)

(1.21)
• 

Two new related functionalities are provided via Library:-RelabelTensorComponents and Library:-RedefineTensorComponent, the first one to have the number of tensor components directly reflected in the names of the components, the second one to redefine only one of these components

Library:-RelabelTensorComponents(C)

Matrix(%id = 18446744078408729774)

(1.22)

 

Suppose now we want to make one of these components equal to 1, say C__2

Library:-RedefineTensorComponent(C[1, 2] = 1)

C[mu, nu] = Matrix(%id = 18446744078270104390)

(1.23)

This nice development is work in collaboration with Pascal Szriftgiser from Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France.

``

 

Download MinimizeTensorComponents.mw

 

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

 

Hello, 

I study mainly subjects that fall under umbrella of number theory, but i have specified a little further in the worksheet. This is really a request for assistance, because in as much as i have met so many brilliant people online via social media etc,  I would always love to meet more, and especially ones who are more experienced in this field. 

 

Basically i am too cheap and old to think about going to a good university, so I am trying to get free advice from the people who have probably completed doctorates in the relevant field. Got to be honest I say.

 

Anyway my contact email is at the top of the attached worksheet.

 

First thing that stood out to me about the distributions produced in this worksheet is how sparse the number of points is for N=17 relative to all the other values of N.

EXAMPLE_FOR_MAPLE3.mw

 

Edit: Another example worksheet added.

MAPLE_EXAMPLE_13.mw

I wanted to use MAPLE to preform symbolic quantum computations. The role
of quantum operators and their tensor product is very important in simplying
the understating of such new calculus at least for the beginners. For instance,
(using "o" for the tensor product and "." for the scalar product, H being the Hadamard
operator on a qubit, I the identity operator, and CNOT the 2 qubit controled not
operator)
1) generating the Bells states |Bxy> two stages of operators are needed
     (CNOT) .  (H o  I)  . |x> o |y>

2) performing quantum teleportation of |psi>
     (H o I o I) . (CNOT o I ) . |psi>o |B00>
    followed by a measurements on the first two qubits for driving the application of
    quantum gates to the third qubit.

All these tensor products of operators can be easily written in MAPLE.

Here is a first version of the ExpandQop procedure that will be usefull the purpose of
expanding correctly the tensor product of two quantum operator expressed in Dirac notation.

I hope this is usefull.

LL 

 

######################################################################
# Author: Louis Lamarche                                             #
#         Institute of Research of Hydro-Quebec (IREQ)               #
#         Science des données et haute performance                   #
#         2018/02/20                                                 #
#                                                                    #
#         Function name: ExpandQop (x)                               #
#               Purpose: Compute the tensor product of two quantum   #
#                        operators in Dirac notations                #
#              Argument: x - a simple quantum operator               #
#   Future improvements: Manage all +, -, *, /, ^, mod  operations   #
#                        in the argument                             #
#               Version: 1.0                                         #
######################################################################
restart;

with(Physics):
interface(imaginaryunit=i):
Setup(mathematicalnotation=true);

[mathematicalnotation = true]

(1)

Setup(quantumoperators={A,B,C,Cn});
Setup(noncommutativeprefix={a,b});

[quantumoperators = {A, B, C, Cn}]

 

[noncommutativeprefix = {a, b}]

(2)

ExpandQop:=proc(x)
    local ret,j,lkb,cbk,rkb,no;
    ret:=1; lkb:=0; cbk:=0; rkb:=0; no:=nops(x);
    if (no > 3 ) then
        for j from 1 to no do
            if (j>1) then
                 if(lkb=0) then
                     if( type(op(j-1,x),Ket) and
                         type(op(j,x),Bra) ) then lkb:=j-1; fi;
                 else
                     if( type(op(j-1,x),Ket) and
                         type(op(j,x),Bra) ) then rkb:=j;   fi;
                 fi;
                 if( type(op(j-1,x),Bra) and type(op(j,x),Ket) )
                                             then cbk:=j;   fi;
            fi;
        end do;
        if ( (lkb < cbk) and (cbk<rkb) ) then
            for j from 1     to lkb   do ret := ret*op(j,x); end do;
            for j from cbk   to no    do ret := ret*op(j,x); end do;
            for j from lkb+1 to cbk-1 do ret := ret*op(j,x); end do;
        else
            ret:=x;
        fi;
    else
        ret:=x;
    fi;
    return ret;
end proc:

# Let A be an operator in a first Hilbert space of dimension n
#  using the associated orthonormal ket and bra vectors
#
#
kets1:=Ket(a1)*Ket(a2)*Ket(a3)*Ket(a4)*Ket(a5):
A:=kets1*Dagger(kets1);


# Let B be an operator in a second Hilbert (Ketspace of dimension m
#  using the associated orthonormal ket and bra vectors
#
#
kets2:=Ket(b1)*Ket(b2)*Ket(b3):
B:=kets2*Dagger(kets2);


# The tensor product of the two operators acts on a n+m third
# Hilbert space   unsing the appropriately ordered ket
# and bra  vectors of the two preceding spaces. The rule for
# building this operator in Dirac notation is as follows,
#
#


print("Maple do not compute the tensor product of operators,");
print("C=A*B gives:");
C:=A*B;

print("ExpandQop(C) gives the expected result:");
Cn:=ExpandQop(C);

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

"Maple do not compute the tensor product of operators,"

 

"C=A*B gives:"

 

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

 

"ExpandQop(C) gives the expected result:"

 

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

(3)

kets3:=kets1*kets2;
bras3:=Dagger(kets3);
print("Matrix element computed with C appears curious");
'bras3.C. kets3'="...";
bras3.C.kets3;
print("Matrix element computed with Cn as expected");
'bras3.Cn.kets3'=bras3.Cn.kets3;

Physics:-`*`(Physics:-Ket(a1), Physics:-Ket(a2), Physics:-Ket(a3), Physics:-Ket(a4), Physics:-Ket(a5), Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3))

 

Physics:-`*`(Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), Physics:-Bra(a5), Physics:-Bra(a4), Physics:-Bra(a3), Physics:-Bra(a2), Physics:-Bra(a1))

 

"Matrix element computed with C"

 

Physics:-`.`(bras3, C, kets3) = "..."

 

Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(b3))*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))

 

"Matrix element computed with Cn as expected"

 

Physics:-`.`(bras3, Cn, kets3) = Physics:-Bracket(Physics:-Bra(a1), Physics:-Ket(a1))^2*Physics:-Bracket(Physics:-Bra(a2), Physics:-Ket(a2))^2*Physics:-Bracket(Physics:-Bra(a3), Physics:-Ket(a3))^2*Physics:-Bracket(Physics:-Bra(a4), Physics:-Ket(a4))^2*Physics:-Bracket(Physics:-Bra(a5), Physics:-Ket(a5))^2*Physics:-Bracket(Physics:-Bra(b1), Physics:-Ket(b1))^2*Physics:-Bracket(Physics:-Bra(b2), Physics:-Ket(b2))^2*Physics:-Bracket(Physics:-Bra(b3), Physics:-Ket(b3))^2

(4)

 


 

Download ExpandQop.mw

 

 

Lesson_on_functions.mws

As the title says, a lesson on functions:  eg the -> operator, f(2), eval, evalf etc 

Using the learning sequence as an alternative to learn problems related to "balance of a body" is shown in this video; thanks to the kindness that Maple offers us in its fundamental programming syntax.

Balance_of_a_body_with_learning_sequence.mw

Lenin Araujo

Ambassador Of Maple

 

 

One case where an "expansion beyond all orders" may be needed is investigating the asymptotic behavior of the difference of two functions with coinciding dominant series.

We are interested in the asymptotic behavior of F(z) for large positive z:

h1 := proc (z) options operator, arrow; hypergeom([1/2], [5/4, 3/2, 7/4], z) end proc; h2 := proc (z) options operator, arrow; (3/4)*sqrt(2*Pi)*hypergeom([1/4], [3/4, 5/4, 3/2], z)/(GAMMA(3/4)*z^(1/4)) end proc; F := proc (z) options operator, arrow; h1(z)-h2(z)+(3/8)*sqrt(Pi)/sqrt(z) end proc

series does not succeed:

series(F(z), z = infinity, 20)

O((1/z)^(23/3))*exp(3/(1/z)^(1/3))

(1)

The reason is that the dominant terms containing the factor exp(3*z^(1/3)) in the two hypergeometric functions cancel exactly, and we have to look for the subdominant terms.

The order of the leading terms can be found from DETools:-formal_sol:

deq1 := FunctionAdvisor(DE, h1(z), f(z))[2, 1]

diff(diff(diff(diff(f(z), z), z), z), z) = -(15/2)*(diff(diff(diff(f(z), z), z), z))/z-(195/16)*(diff(diff(f(z), z), z))/z^2+(1/32)*(32*z-105)*(diff(f(z), z))/z^3+(1/2)*f(z)/z^3

(2)

DETools:-formal_sol(deq1, f(z), z = infinity, order = 0)

[(1/z)^(1/2), -exp(-3/(-1/z)^(1/3))/z, -exp(3*(-1)^(1/3)/(-1/z)^(1/3))/z, -exp(-3*(-1)^(2/3)/(-1/z)^(1/3))/z]

(3)

As expected, one of the solutions (the third one for positive z) contains the exp(3*z^(1/3)) factor, the leading term being of the order exp(3*z^(1/3))/z.

Another, subdominant, solution is algebraic and, in fact, is a series containing only one term, as 1/z^(1/2) is an exact solution. It will turn out that the algebraic part in F(z) also cancels out.

Thus we have to look for the subsubdominant terms, which contain decaying exponentials. We will accomplish this by applying the steepest descent method to the integral representations of h1(z) and h2(z).

ms := convert([h1(z), h2(z)], MeijerG)

[(3/32)*Pi*2^(1/2)*MeijerG([[1/2], []], [[0], [-1/4, -1/2, -3/4]], -z), (3/32)*2^(1/2)*Pi*MeijerG([[3/4], []], [[0], [1/4, -1/4, -1/2]], -z)/z^(1/4)]

(4)

m2g := proc (m, y) local a, b, c, d; a, b := op(op(1, m)); c, d := op(op(2, m)); -((1/2)*I)*mul(`~`[GAMMA](`~`[`-`](1+y, a)))*mul(`~`[GAMMA](`~`[`-`](c, y)))*op(3, m)^y/(Pi*mul(`~`[GAMMA](`~`[`-`](b, y)))*mul(`~`[GAMMA](`~`[`-`](1+y, d)))) end proc

gs := applyrule(conditional(e::anything, _op(0, e) = MeijerG) = 'm2g(e, y)', ms)

[-((3/64)*I)*2^(1/2)*GAMMA(1/2+y)*GAMMA(-y)*(-z)^y/(GAMMA(5/4+y)*GAMMA(3/2+y)*GAMMA(7/4+y)), -((3/64)*I)*2^(1/2)*GAMMA(1/4+y)*GAMMA(-y)*(-z)^y/(z^(1/4)*GAMMA(3/4+y)*GAMMA(5/4+y)*GAMMA(3/2+y))]

(5)

gs[2] := combine(eval(gs[2], [1/z^(1/4) = exp(I*Pi*(1/4))/(-z)^(1/4), y = y+1/4]), power)

-((3/64)*I)*GAMMA(1/2+y)*((1/2)*2^(1/2)+((1/2)*I)*2^(1/2))*(-z)^y*GAMMA(-1/4-y)*2^(1/2)/(GAMMA(1+y)*GAMMA(3/2+y)*GAMMA(7/4+y))

(6)

h1(z) and h2(z)are the integrals of gs[1] and of gs[2] over the same path, which is a loop encircling the poles ofGAMMA(-y) and of GAMMA(-1/4-y). Now a standard technique is to extend the integration contour far to the left, while still keeping both endpoints at "+infinity". Then the arguments of the gamma functions can be made large everywhere on the integration path, and the gamma functions can be replaced by their asymptotic approximations.

When moving the contour, we have to take into account the pole of the integrand at y = -1/2. The other poles of GAMMA(1/2+y) will be cancelled by the zeros of 1/GAMMA(3/2+y), which is why the algebraic part of the expansion will contain the single term of the order 1/z^(1/2).

This is the negative of the third term in F(z):

`assuming`([simplify((2*Pi*I)*residue(gs[1]-gs[2], y = -1/2))], [z > 0])

-(3/8)*Pi^(1/2)/z^(1/2)

(7)

Expanding the gamma functions produces terms containing exp(-I*Pi*y) and exp(I*Pi*y)

`assuming`([simplify(convert(MultiSeries:-series((gs[1]-gs[2])/z^y, y = -infinity, 1), polynom))], [z > 0]); collect(convert(%, exp), exp)

(-3/64-(3/64)*I)*(-1/y)^(1/2)*exp(-3*(ln(-y)-1)*y)*exp(-I*Pi*y)/(y^3*Pi^(1/2))+(3/64-(3/64)*I)*(-1/y)^(1/2)*exp(-3*(ln(-y)-1)*y)*exp(I*Pi*y)/(y^3*Pi^(1/2))

(8)

As we shall see, those terms have saddle points y0(z) = exp(`&+-`((1/3)*(2*Pi*I)))*z^(1/3) located in the left half-plane and contribute exponentially small factors exp(3*y0(z)). The terms for which the saddle point would be located at y = z^(1/3) have cancelled out, thus cancelling the exponentially large contributions. Another possible way to achieve the same result was to write h1(z)-h2(z) as a single Meijer G-function -(3/32)*MeijerG([[1/2], []], [[-1/4, 0], [-3/4, -1/2]], z).

We write the first term above in the form g(y)*exp(f(y)):

f := proc (z, y) options operator, arrow; -3*y*(ln(-y)-1)-I*Pi*y+y*ln(z) end proc

g := proc (y) options operator, arrow; (-3/64-(3/64)*I)*sqrt(-1/y)/(sqrt(Pi)*y^3) end proc

diff(f(z, y), y)

-3*ln(-y)-I*Pi+ln(z)

(9)

For this to become zero, we need argument(-y) = -(1/3)*Pi, and thus y = exp((1/3)*(2*I)*Pi)*z^(1/3). We can visualize the paths where the imaginary part of f(z, y) stays constant. The path of the steepest descent is the one that goes through the saddle point in the direction exp(I*Pi*(1/3)); the blue color indicates smaller values of the real part of f(z, y):

y0 := proc (z) options operator, arrow; exp(((2/3)*I)*Pi)*z^(1/3) end proc

(proc () local z; z := 2; plots:-display(plots:-contourplot(Re(f(z, u+I*v)), u = -5 .. 5, v = -5 .. 5, contours = ([seq])(Re(f(z, y0(z)))+i, i = -30 .. 6, 6), filledregions, coloring = [blue, red], grid = [100, 100]), plots:-implicitplot(Im(f(z, u+I*v)-f(z, y0(z))), u = -5 .. 5, v = -5 .. 5, gridrefine = 5, color = green), plot([cos((1/3)*Pi)*xi+Re(y0(z)), sin((1/3)*Pi)*xi+Im(y0(z)), xi = -3 .. 3], linestyle = dot, color = white), axes = boxed) end proc)()

 

The real part of f(z, y) has a maximum along this path at y0(z).

`assuming`([(`@`(`@`(simplify, evalc), series))(f(z, y0(z)+exp(I*Pi*(1/3))*xi), xi = 0, 3)], [z > 0]); quad := convert(%, polynom)

series((3/2)*(-1+I*3^(1/2))*z^(1/3)-((3/2)/z^(1/3))*xi^2+O(xi^3),xi,3)

(10)

Now we can compute the lead asymptotic term contributed by the saddle point y0(z):

lt1 := `assuming`([(`@`(simplify, evalc))(g(y0(z))*exp(I*Pi*(1/3))*(int(exp(quad), xi = -infinity .. infinity)))], [z > 0])

-(1/64)*exp(-(3/2)*z^(1/3))*3^(1/2)*((-1+I)*cos((3/2)*z^(1/3)*3^(1/2))+(-1-I)*sin((3/2)*z^(1/3)*3^(1/2)))*2^(1/2)/z

(11)

We repeat the same procedure for the second term of the integrand.

f := proc (z, y) options operator, arrow; -3*y*(ln(-y)-1)+I*Pi*y+y*ln(z) end proc

g := proc (y) options operator, arrow; (3/64-(3/64)*I)*sqrt(-1/y)/(sqrt(Pi)*y^3) end proc

diff(f(z, y), y)

-3*ln(-y)+I*Pi+ln(z)

(12)

y0 := proc (z) options operator, arrow; exp(-((2/3)*I)*Pi)*z^(1/3) end proc

The direction should be chosen as exp((1/3)*(2*I)*Pi) to be consistent with the direction of the integration contour, which goes from the lower to the upper half-plane.

lterm := proc (gy, fy, eq, dir) options operator, arrow; (eval(gy*exp(fy), eq))*dir*sqrt(-2*Pi/((eval(diff(fy, `$`(y, 2)), eq))*dir^2)) end proc

lt2 := `assuming`([(`@`(simplify, evalc))(lterm(g(y), f(z, y), y = y0(z), exp((1/3)*(2*I)*Pi)))], [z > 0])

(1/64)*exp(-(3/2)*z^(1/3))*((1+I)*cos((3/2)*z^(1/3)*3^(1/2))+(1-I)*sin((3/2)*z^(1/3)*3^(1/2)))*3^(1/2)*2^(1/2)/z

(13)

Combining the two results yields the leading term of F(z). The next terms can be obtained by expanding gs[1] and gs[2] to higher orders.

Fasympt := unapply(simplify(lt1+lt2), z)

proc (z) options operator, arrow; (1/32)*exp(-(3/2)*z^(1/3))*3^(1/2)*2^(1/2)*(cos((3/2)*z^(1/3)*3^(1/2))+sin((3/2)*z^(1/3)*3^(1/2)))/z end proc

(14)

(proc () Digits := 50; plot(`~`[`*`](exp((3/2)*z^(1/3)), [F(z), Fasympt(z)]), z = 1000 .. 10000, linestyle = [solid, dot], thickness = [1, 5], axes = frame) end proc)()

 

Download steep.mw

Just a simple little worksheet to see if I have enough propane to heat my house for the rest of the winter.


 

Do I have enough propane for the winter?

NULL

I've taken some measurements from my propane tank throughout the winter.  Now we can use Maple to see if we have enough to last the rest of the winter.

``

a := [["nov 27, 2017", 73.5], ["dec 9, 2017", 72], ["dec 16, 2017", 69], ["dec 31, 2017", 62], ["jan 12, 2018", 60], ["jan 19, 2018", 56], ["jan 26, 2018", 54], ["feb 4,2018", 51]]

[["nov 27, 2017", 73.5], ["dec 9, 2017", 72], ["dec 16, 2017", 69], ["dec 31, 2017", 62], ["jan 12, 2018", 60], ["jan 19, 2018", 56], ["jan 26, 2018", 54], ["feb 4,2018", 51]]

(1)

with(Finance)  ``

pts := [seq([DayCount(a[1, 1], a[i, 1]), a[i, 2]], i = 1 .. nops(a))]

[[0, 73.5], [12, 72], [19, 69], [34, 62], [46, 60], [53, 56], [60, 54], [69, 51]]

(2)

with(plots)

listplot(pts)

 

Adding a 30% and 20% level to the graph.  We probably shouldn't be too worried about the cold in June so DayCount("Nov 27, 2017", "Jun 1, 2018") = 186 we'll extend these reference lines out to 186.

plot({pts, [[0, 20], [186, 20]], [[0, 30], [186, 30]]}, view = [default, 0 .. 80])

 

 

30% is the recommended level your propane company wants you to fill up at.  The technician who installed the tank said 20% is all right.  It's up to you if you want to go to 10% but if you run out of propane the company has to come in and do a leak test on your system which is an added cost you don't want.  So let's predict at what point we need to start worrying about filling up our propane tank.  To do that, of course, all we need is a forecast line.  For that we'll just calculate a best fit.

 

a1 := [seq(DayCount(a[1, 1], a[i, 1]), i = 1 .. nops(a))]

[0, 12, 19, 34, 46, 53, 60, 69]

(3)

a2 := a[() .. (), 2]

[73.5, 72, 69, 62, 60, 56, 54, 51]

(4)

X := convert(a1, Vector)

Y := convert(a2, Vector)

with(Statistics)

L1 := LinearFit([1, x], X, Y, x)

HFloat(74.79237702730747)-HFloat(0.34416046490941915)*x

(5)

Plotting it all together

plot({L1, pts, [[0, 20], [186, 20]], [[0, 30], [186, 30]]}, x = 0 .. 200, y = 0 .. 80, labels = ["Days", ""], tickmarks = [default, [seq(10*i = cat(10*i, "%"), i = 1 .. 8)]])

 

Projecting the line to 30% we get

solve(L1 = 30)

130.1496877

(6)

AdvanceDate(a[1, 1], trunc(solve(L1 = 30)))

Record(monthDay = 6, month = 4, year = 2018, format = "%B %e, %Y", ModulePrint = proc (m) Finance:-FormatDate(m) end proc)

(7)

April is still a bit chilly so maybe if we wait until 20%, of course it's getting warmer all this time so our usage should go down.  

AdvanceDate(a[1, 1], trunc(solve(L1 = 20)))

Record(monthDay = 5, month = (), year = 2018, format = "%B %e, %Y", ModulePrint = proc (m) Finance:-FormatDate(m) end proc)

(8)

It isn't warm enough to turn off the furnace yet but it looks like we'll have enough to get us into the warm months

AdvanceDate(a[1, 1], trunc(solve(L1 = 10)))

Record(monthDay = 3, month = 6, year = 2018, format = "%B %e, %Y", ModulePrint = proc (m) Finance:-FormatDate(m) end proc)

(9)

We'll hit 10% well into late spring and almost right into summer of course it's a rough estimate however it looks like we won't have to fill up during the high price winter season.  I can tell my wife to relax, we should have enough propane for the winter.

 

 

NULL


 

Download Propane_usage-.mw

On the example of a manipulator with three degrees of freedom.
A mathematical model is created that takes into account degrees of freedom of the manipulator and the trajectory of the movement from the initial point to the final one (in the figure, the ends of the red curve). In the text of the program, these are the equations fi, i = 1..5.
Obviously, the straight line could be the simplest trajectory, but we will consider a slightly different variant. The solution of the system of equations is the coordinates of the points of the manipulator (x1, x2, x3) and (x4, x5, x6) in all trajectory. After that, knowing the lengths of the links and the coordinates of the points at each moment of time, any angles of the manipulator are calculated. The same selected trajectory is reproduced from these angles. The possible angles are displayed by black color.
All the work on creating a mathematical model and calculating the angles can be done without the manipulator itself, is sufficient to have only the instruction with technical characteristics.
To display some angles, the procedure created by vv is used.
MAN_2.mw

This post is devoted to the rigorous proof of Miquel's five circles theorem, which I learned about from this question. The proof is essentially very simple and takes only 15 lines of code. The figure below, in which all the labels coincide with the corresponding names in the code, illustrates the basic ideas of the code. First, we symbolically define common points of intersection of blue circles with a red unit circle  (these parameters  s1 .. s5  are the polar coordinates of these points). All other parameters of this configuration can be expressed through them. Then we find the centers  M  and  N  of two circles. Then we find the coordinates of the point  K  from the condition that  CK  is perpendicular to  MN . Then we find the point  and using the result obtained, we easily find the coordinates  of all the points  A1 .. A5. Then we find the coordinates of the point   P  as the point of intersection of the lines  A1A2  and  A3A4 . Finally, we verify that the point  P  lies on a circle with center at the point  N , which completes the proof.

                      

 

Below - the code of the proof. Note that the code does not use any special (in particular geometric) packages, only commands from the Maple kernel. I usually try any geometric problems to solve in this style, it is more reliable,  and often shorter.

restart;
t1:=s1/2+s2/2: t2:=s2/2+s3/2:
M:=[cos(t1),sin(t1)]: N:=[cos(t2),sin(t2)]:
C:=[cos(s2),sin(s2)]: K:=(1-t)*~M+t*~N:
CK:=K-C: MN:=M-N:
t0:=simplify(solve(CK[1]*MN[1]+CK[2]*MN[2]=0, t)):
E:=combine(simplify(C+2*eval(CK,t=t0))):
s0:=s5-2*Pi: s6:=s1+2*Pi:
assign(seq(A||i=eval(E,[s2=s||i,s1=s||(i-1),s3=s||(i+1)]), i=1..5)):
Dist:=(p,q)->sqrt((p[1]-q[1])^2+(p[2]-q[2])^2):
LineEq:=(P,Q)->(y-P[2])*(Q[1]-P[1])=(x-P[1])*(Q[2]-P[2]):
Line1:=LineEq(A1,A2):
Line2:=LineEq(A3,A4):
P:=combine(simplify(solve({Line1,Line2},[x,y])))[]:
Circle:=(x-N[1])^2+(y-N[2])^2-Dist(N,C)^2:
is(eval(Circle, P)=0);  
# The final result

                                                                    true


It may seem that this proof is easy to repeat manually. But this is not so. Maple brilliantly coped with very cumbersome trigonometric transformations. Look at the coordinates of point  , expressed through the initial parameters  s1 .. s5 :

simplify(eval([x,y], P));  # The coordinates of the point  P

  

  

 

ProofMiquel.mw

My September 9, 2016, blog post ("Next Number" Puzzles) pointed out the meaninglessness of the typical "next-number" puzzle. It did this by showing that two such puzzles in the STICKELERS column by Terry Stickels had more than one solution. In addition to the solution proposed in the column, another was found in a polynomial that interpolated the given members of the sequence. Of course, the very nature of the question "What is the next number?" is absurd because the next number could be anything. At best, such puzzles should require finding a pattern for the given sequence, admitting that there need not be a unique pattern.

The STICKELERS column continued to publish additional "next-number" puzzles, now no longer of interest. However, the remarkable puzzle of December 30, 2017, caused me to pull from the debris on my retirement desk the puzzle of July 15, 2017, a puzzle I had relegated to the accumulating dust thereon.

The members of the given sequence appear across the top of the following table that reproduces the graphic used to provide the solution.

It turns out that the pattern in the graphic can be expressed as 100 – (-1)k k(k+1)/2, k=0,…, a pattern Maple helped find. By the techniques in my earlier blog, an alternate pattern is expressed by the polynomial

which interpolates the nodes (1, 100), (2, 101) ... so that f(8) = -992.

The most recent puzzle consists of the sequence members 0, 1, 8, 11, 69, 88; the next number is given as 96 because these are strobogrammatic numbers, numbers that read the same upside down. Wow! A sequence with apparently no mathematical structure! Is the pattern unique? Well, it yields to the polynomial

which can also be expressed as

Hence, g(x) is an integer for any nonnegative integer x, and g(6) = -401, definitely not a strobogrammatic number. However, I do have a faint recollection that one of Terry's "next-number" puzzles had a pattern that did not yield to interpolation. Unfortunately, the dust on my desk has not yielded it up.
 

A few days ago, I drew attention to the question in which OP talked about the generation of triangles in a plane, for which the lengths of all sides, the area and radius of the inscribed circle are integers. In addition, all vertices must have different integer coordinates (6 different integers), the lengths of all sides are different and the triangles should not be rectangular. I prepared the answer to this question, but the question disappeared somewhere, so I designed my answer as a separate post.

The triangles in the plane, for which the lengths of all sides and the area  are integers, are called as Heronian triangles. See this very interesting article in the wiki about such triangles
https://en.wikipedia.org/wiki/Integer_triangle#Heronian_triangles

The procedure finds all triangles (with the fulfillment of all conditions above), for which the lengths of the two sides are in the range  N1 .. N2 . The left side of the range is an optional parameter (by default  N1=5). It is not recommended to take the length of the range more than 100, otherwise the operating time of the procedure will greatly increase. The procedure returns the list in which each triangle is represented by a list of  [list of coordinates of the vertices, area, radius of the inscribed circle, list of lengths of the sides]. Without loss of generality, one vertex coincides with the origin (obviously, by a shift it is easy to place it at any point). 

The procedure works as follows: one vertex at the origin, then the other two must lie on circles with integer and different radii  x^2+y^2=r^2. Using  isolve  command, we find all integer points on these circles, and then in the for loops we select the necessary triangles.


 

 

restart;
HeronianTriangles:=proc(N2::posint,N1::posint:=5)
local k, r, S, L, Ch, Dist, IsOnline, c, P, p, A, B, C, a, b, s, ABC, cc, s1, T ;
uses combinat, geometry;
if N2<N1 then error "Should be N2>=N1" fi;
if N2<34 then return [] fi;
k:=0:
for r from max(N1,5) to N2 do
S:=[isolve(x^2+y^2=r^2)];
if nops(S)>4 then k:=k+1; L[k]:=select(s->s[1]<>0 and s[2]<>0,map(t->rhs~(convert(t,list)), S)); fi;
od:
L:=convert(L, list):
if type(L[1],symbol) then return [] fi;

Ch:=combinat:-choose([$1..nops(L)], 2):
Dist:=(A::list,B::list)->simplify(sqrt((A[1]-B[1])^2+(A[2]-B[2])^2));
IsOnline:=(A::list,B::list)->`if`(A[1]*B[2]-A[2]*B[1]=0, true, false);
k:=0:
for c in Ch do
for A in L[c[1]] do
for B in L[c[2]] do
if not IsOnline(A,B) and nops({A[],B[]})=4 then if type(Dist(A,B),posint) then
 k:=k+1; P[k]:=[A,B] fi; fi;
od: od: od:
P:=convert(P, list):
if type(P[1],symbol) then return [] fi;

k:=0:
for p in P do
point('A',0,0), point('B',p[1]), point('C',p[2]);
a:=simplify(distance('A','B')); b:=simplify(distance('A','C')); c:=simplify(distance('B','C'));
s:=sort([a,b,c]); s1:={a,b,c};
triangle(ABC,['A','B','C']);
incircle(cc,ABC);
r:=radius(cc);
if type(r,integer) and s[3]^2<>s[1]^2+s[2]^2 and nops(s1)=3 then k:=k+1; T[k]:=[[[0,0],p[]],area(ABC),r, [a,b,c]] fi;
od:
T:=convert(T,list);
if type(T[1],symbol) then return [] fi;
T;
end proc:

Examples of use of the procedure  HeronianTriangles

T:=HeronianTriangles(100): # All the Geronian triangles, whose lengths of two sides do not exceed 100
nops(T);

256

(1)

Tp:=select(p->p[1,2,1]>0 and p[1,2,2]>0 and p[1,3,1]>0 and p[1,3,2]>0, T);

[[[[0, 0], [16, 30], [28, 21]], 252, 6, [34, 35, 15]], [[[0, 0], [30, 16], [21, 28]], 252, 6, [34, 35, 15]], [[[0, 0], [21, 28], [15, 36]], 168, 4, [35, 39, 10]], [[[0, 0], [28, 21], [36, 15]], 168, 4, [35, 39, 10]], [[[0, 0], [27, 36], [13, 84]], 900, 10, [45, 85, 50]], [[[0, 0], [36, 27], [84, 13]], 900, 10, [45, 85, 50]], [[[0, 0], [33, 44], [48, 36]], 462, 7, [55, 60, 17]], [[[0, 0], [44, 33], [36, 48]], 462, 7, [55, 60, 17]], [[[0, 0], [33, 44], [96, 28]], 1650, 15, [55, 100, 65]], [[[0, 0], [44, 33], [28, 96]], 1650, 15, [55, 100, 65]], [[[0, 0], [16, 63], [72, 21]], 2100, 20, [65, 75, 70]], [[[0, 0], [63, 16], [21, 72]], 2100, 20, [65, 75, 70]], [[[0, 0], [39, 52], [18, 80]], 1092, 12, [65, 82, 35]], [[[0, 0], [52, 39], [80, 18]], 1092, 12, [65, 82, 35]], [[[0, 0], [32, 60], [56, 42]], 1008, 12, [68, 70, 30]], [[[0, 0], [60, 32], [42, 56]], 1008, 12, [68, 70, 30]], [[[0, 0], [42, 56], [30, 72]], 672, 8, [70, 78, 20]], [[[0, 0], [56, 42], [72, 30]], 672, 8, [70, 78, 20]]]

(2)

Tr:=map(p->p+[2,1],Tp[1,1]);
with(geometry):
point(A,Tr[1]), point(B,Tr[2]), point(C,Tr[3]):
triangle(ABC,[A,B,C]):
simplify(distance(A,B)), simplify(distance(A,C)), simplify(distance(B,C));
local O:
incircle(c,ABC, centername=O):
draw([A,B,C, ABC, c(color=blue)], color=red, thickness=2, symbol=solidcircle, tickmarks = [spacing(1)$2], gridlines, scaling=constrained, view=[0..31,0..33], size=[800,550], printtext=true, font=[times, 18], axesfont=[times, 10]);

[[2, 1], [18, 31], [30, 22]]

 

34, 35, 15

 

 



Examples of triangles with longer sides

T:=HeronianTriangles(1000,980):  # All the Geronian triangles, whose lengths of two sides lie in the range  980..1000
nops(T);

56

(3)

Tp:=select(p->p[1,2,1]>0 and p[1,2,2]>0 and p[1,3,1]>0 and p[1,3,2]>0, T);  # Triangles lying in the first quarter x>0, y>0
nops(%);

[[[[0, 0], [540, 819], [680, 714]], 85680, 80, [981, 986, 175]], [[[0, 0], [819, 540], [714, 680]], 85680, 80, [981, 986, 175]], [[[0, 0], [216, 960], [600, 800]], 201600, 168, [984, 1000, 416]], [[[0, 0], [960, 216], [800, 600]], 201600, 168, [984, 1000, 416]], [[[0, 0], [380, 912], [324, 945]], 31806, 31, [988, 999, 65]], [[[0, 0], [912, 380], [945, 324]], 31806, 31, [988, 999, 65]], [[[0, 0], [594, 792], [945, 324]], 277992, 216, [990, 999, 585]], [[[0, 0], [792, 594], [324, 945]], 277992, 216, [990, 999, 585]]]

 

8

(4)

 


 

Download Integer_Triangle1.mw

Edit.

Implementation of Maple apps for the creation of mathematical exercises in
engineering

In this research work has allowed to show the implementation of applications developed in the Maple software for the creation of mathematical exercises given the different levels of education whether basic or higher.
For the majority of teachers in this area, it seems very difficult to implement apps in Maple; that is why we show the creation of exercises easily and permanently. The purpose is to get teachers from our institutions to use applications ready to be evaluated in the classroom. The results of these apps (applications with components made in Maple) are supported on mobile devices such as tablets and / or laptops and taken to the cloud to be executed online from any computer. The generation of patterns is a very important alternative leaving aside random numbers, which would allow us to lose results
onscreen. With this; Our teachers in schools or universities would evaluate their students in parallel on the blackboard without losing the results of any student and thus achieve the competencies proposed in the learning sessions.
In these apps would be the algorithms for future research updates and integrated with systems in content management. Therefore what we show here is extremely important for the evaluation on the blackboard in bulk to students without losing any scientific criteria.

FAST_UNT_2018.mw

FAST_UNT_2018.pdf

Lenin Araujo Castillo

Ambasador of Maple

 

In the beginning, Maple had indexed names, entered as x[abc]; as early as Maple V Release 4 (mid 1990s), this would display as xabc. So, x[1] could be used as x1, a subscripted variable; but assigning a value to x1 created a table whose name was x. This had, and still has, undesirable side effects. See Table 0 for an illustration in which an indexed variable is assigned a value, and then the name of the concomitant table is also assigned a value. The original indexed name is destroyed by these steps.
 

 

At one time this quirk could break commands such as dsolve. I don't know if it still does, but it's a usage that those "in the know" avoid. For other users, this was a problem that cried out for a solution. And Maplesoft did provide such a solution by going nuclear - it invented the Atomic Variable, which subsumed the subscript issue by solving a larger problem.

The larger problem is this: Arbitrary collections of symbols are not necessarily valid Maple names. For example, the expression  is not a valid name, and cannot appear on the left of an assignment operator. Values cannot be assigned to it. The Atomic solution locks such symbols together into a valid name originally called an Atomic Identifier but now called an Atomic Variable. Ah, so then xcan be either an indexed name (table entry) or a non-indexed literal name (Atomic Variable). By solving the bigger problem of creating assignable names, Maplesoft solved the smaller problem of subscripts by allowing literal subscripts to be Atomic Variables.

It is only in Maple 2017 that all vestiges of "Identifier" have disappeared, replaced by "Variable" throughout. The earliest appearance I can trace for the Atomic Identifier is in Maple 11, but it might have existed in Maple 10. Since Maple 11, help for the Atomic Identifier is found on the page 

 

In Maple 17 this help could be obtained by executing help("AtomicIdentifier"). In Maple 2017, a help page for AtomicVariables exists.

In Maple 17, construction of these Atomic things changed, and a setting was introduced to make writing literal subscripts "simpler." With two settings and two outcomes for a "subscripted variable" (either indexed or non-indexed), it might be useful to see the meaning of "simpler," as detailed in the worksheet AfterMath.mw.

Hi,

2 examples to explore components to illustrated Fractals
 

NULL

Flocon de VonKoch

 

 

 

NULL

NULL

NULL

NULL

 

 

``


 

Download AppliFloconVonKoch.mw

It passed through my mind it would be interesting to collect the links to the most relevant Mapleprimes posts about Quantum Mechanics using the Physics package of the last couple of years, to have them all accessible from one place. These posts give an idea of what kind of computation is already doable in quantum mechanics, how close is the worksheet input to what we write with paper and pencil, and how close is the typesetting of the output to what we see in textbooks.

At the end of each page linked below, you will see another link to download the corresponding worksheet, that you can open using Maple (say the current version or the version 1 or 2 years ago).

This other set of three consecutive posts develops one problem split into three parts:

This other link is interesting as a quick and compact entry point to the use of the Physics package:

There is an equivalent set of Mapleprimes posts illustrating the Physics package tackling problems in General Relativity, collecting them is for one other time.

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

First 27 28 29 30 31 32 33 Last Page 29 of 77