Maple Questions and Posts

These are Posts and Questions associated with the product, Maple

Hello there!

I am facing a strange problem. I am defining a function as a finite series which involves the constants c[k] and the variable x. when I try to calculate the value of the function at zero it should give me c[1] but its giving me 0. it seems a bug in maple. i have tried it on two different versions.i.e. 2016 and 17


Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/bug.mw .
 

Download bug.mw

 

Automatic handling of collision of tensor indices in products

 

 

The design of products of tensorial expressions that have contracted indices got enhanced. The idea: repeated indices in certain subexpressions are actually dummies. So suppose T[a, b] and B[b] are tensors, then in T[trace] = T[a, `~a`], a is just a dummy, therefore T[a, `~a`]*B[a] = T[b, `~b`]*B[a] is a well defined object. The new design automatically maps input like T[a, `~a`]*B[a] into T[b, `~b`]*B[a]. This significantly simplifies the manipulation of indices when working with products of tensorial expressions: each tensorial expression being multiplied has its repeated indices automatically transformed into different ones when they would collide with the free or repeated indices of the other expressions being multiplied.

 

This functionality is available within the Physics update distributed at the Maplesoft R&D Physics webpage (but for what you see under Preview that makes use of a new kernel feature of the Maple version under development).

 

restart

with(Physics); Setup(spacetimeindices = lowercaselatin, quiet)

[spacetimeindices = lowercaselatin]

(1)

Define(T[a, b], B[b])

`Defined objects with tensor properties`

 

{B[b], Physics:-Dgamma[a], Physics:-Psigma[a], T[a, b], Physics:-d_[a], Physics:-g_[a, b], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c, d]}

(2)

This shows the automatic handling of collision of indices

T[a, a]*B[a]

T[b, `~b`]*B[a]

(3)

T[a, a]^2

T[a, `~a`]*T[b, `~b`]

(4)

``

Preview only in the upcomming version under development

 

Consider now the case of three tensors

Define(A[a], C[a])

`Defined objects with tensor properties`

 

{A[a], B[b], C[a], Physics:-Dgamma[a], Physics:-Psigma[a], T[a, b], Physics:-d_[a], Physics:-g_[a, b], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c, d]}

(5)

A[a]*B[a]*C[a]

A[a]*B[a]*C[a]

(6)

The product above has indeed the index a repeated more than once, therefore none of its occurrences got automatically transformed into contravariant in the output, and Check  detects the problem interrupting with an error  message

Check(A[a]*B[a]*C[a])

Error, (in Physics:-Check) wrong use of the summation rule for repeated indices: `a repeated 3 times`, in A[a]*B[a]*C[a]

 

 

However, it is now also possible to indicate, using parenthesis, that the product of two of these tensors actually form a subexpression, so that the following two tensorial expressions are well defined, and the colliding dummy is automatically replaced making that explicit

A[a]*B[a]*C[a]

A[b]*B[`~b`]*C[a]

(7)

A[a]*B[a]*C[a]

A[a]*B[b]*C[`~b`]

(8)

 

 

This change in design makes concretely simpler the use of indices in that it eliminates the need for manually replacing dummies. For example, consider the tensorial expression for the angular momentum in terms of the coordinates and momentum vectors, in 3 dimensions

 

Setup(coordinates = cartesian, dimension = 3, metric = euclidean, quiet)

[coordinatesystems = {X}, dimension = 3, metric = {(1, 1) = 1, (2, 2) = 1, (3, 3) = 1}]

(9)

Define L[j], p[k] respectively representing angular and linear momentum

Define(L[j], p[k])

`Defined objects with tensor properties`

 

{Physics:-Dgamma[a], L[j], Physics:-Psigma[a], Physics:-d_[a], Physics:-g_[a, b], p[k], Physics:-KroneckerDelta[a, b], Physics:-LeviCivita[a, b, c], Physics:-SpaceTimeVector[a](X)}

(10)

Introduce the tensorial expression for L[a]

L[a] = LeviCivita[a, b, c]*X[b]*p[c]

L[a] = Physics:-LeviCivita[a, b, c]*Physics:-SpaceTimeVector[b](X)*p[c]

(11)

The left-hand side has one free index, a, while the right-hand side has two dummy indices b and c

Check(L[a] = Physics[LeviCivita][a, b, c]*Physics[SpaceTimeVector][b](X)*p[c], all)

`The repeated indices per term are: `[{`...`}, {`...`}, `...`]*`; the free indices are: `*{`...`}

 

([{}], {a}) = ([{b, c}], {a})

(12)

If we want to compute`#mrow(msup(mfenced(mover(mi("L"),mo("→")),open = "‖",close = "‖"),mn("2")),mo("="),msubsup(mi("L"),mi("a"),mn("2")))`we can now take the square of (11) directly, and the dummy indices on the right-hand side are automatically handled, there is now no need to manually substitute the repeated indices to avoid their collision

(L[a] = Physics[LeviCivita][a, b, c]*Physics[SpaceTimeVector][b](X)*p[c])^2

L[a]^2 = Physics:-LeviCivita[a, b, c]*Physics:-SpaceTimeVector[b](X)*p[c]*Physics:-LeviCivita[a, d, e]*Physics:-SpaceTimeVector[d](X)*p[e]

(13)

The repeated indices on the right-hand side are now a, b, c, d, e

Check(L[a]^2 = Physics[LeviCivita][a, b, c]*Physics[SpaceTimeVector][b](X)*p[c]*Physics[LeviCivita][a, d, e]*Physics[SpaceTimeVector][d](X)*p[e], all)

`The repeated indices per term are: `[{`...`}, {`...`}, `...`]*`; the free indices are: `*{`...`}

 

([{a}], {}) = ([{a, b, c, d, e}], {})

(14)

NULL


 

Download AutomaticHandlingCollisionOfTensorIndices.mw

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

The problem, find the general solution of y '' + 4y = t cos (2t).

Maple input:

de:=diff(y(t),t,t)+4*y(t)=t*cos(2*t);
sol:=dsolve(de,y(t));

Maple output:

sin(2*t)*_C2+cos(2*t)*_C1+(1/8)*t^2*sin(2*t)-(1/64)*sin(2*t)+(1/16)*t*cos(2*t)

The odd thing is the inclusion of the term -(1/64)*sin(2*t). It is not incorrect since you can collect this term with sin(2*t)*_C2. Is there a reason why it's there, and how can i remove it without inspecting it? Note that Wolfram doesn't have this extra term.

https://www.wolframalpha.com/input/?i=solve+y%27%27+%2B+4y+%3D+t*cos(2*t)

I attached the worksheet and added a more detailed calculation.

diffeq.mw

 

Hi,

I was wondering if anyone has a clever way to code the Cayley Omega process?
For those who are wondering, the Omega process is a differential operator. Given an n-dimensional space (x[1],x[2],x[3],...,x[n]), and n forms Q[1](x[1][1],x[1][2],x[1][3],etc) Q[2](x[2][1],x[2][2],x[2][3],etc) ... Q[n](x[n][1],x[n][1],x[n][1],etc), the operator is the determinant of the matrix who entries are the partial differential operators del/delx[i][j].

Thoughts? Suggestions?

 

Thanks.

Here is a major upgrade of the procedure I submitted in february.

https://www.mapleprimes.com/posts/209030-Procedure-For-Computing-The-Tensor-Product

There is a line after the procedure to save it in the file "ExpandQop.m"
In future post I will use it in order to minimize the size of the examples.

Louis Lamarche
 

######################################################################
# NOTICE                                                             #
# Author: Louis Lamarche                                             #
#         Institute of Research of Hydro-Quebec (IREQ)               #
#         Science des données et haute performance                   #
#         2018, March 7                                              #
#                                                                    #
# Function name: ExpandQop (x)                                       #
#       Purpose: Compute the tensor product of two quantum           #
#                operators in Dirac notations                        #
#      Argument: x: a quantum operator                               #
#  Improvements: Manage all +, -, *, /, ^, mod  operations           #
#                in the argument                                     #
#       Version: 2                                                   #
#                                                                    #
#  Copyrigth(c) Hydro-Quebec.                                        #
#        Note 1: Permission to use this softwate is granted if you   #
#                acknowledge its author and copyright                #
#        Note 2: Permission to copy this softwate is granted if you  #
#                leave this 21 lines notice intact. Thank you.       #
######################################################################
restart;

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

[mathematicalnotation = true]

(1)

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

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

 

[noncommutativeprefix = {a, b, q}]

(2)

opexp:= op(0,10^x):            # exponentiation id
opnp := op(0,10*x):            # normal product id
optp := op(0,Ket(q0)*Ket(q1)): # tensor product id
opdiv:= `Fraction`:            # fraction       id          
opsym:= op(0,x):               # symbol         id
opint:= op(0,10):              # integer        id
opflt:= op(0,10.0):            # float          id
opcpx:= op(0,i):               # complex        id
opbra:= op(0,Bra(q)):          # bra            id
opket:= op(0,Ket(q)):          # ket            id
opmod:= op(0, a mod b):        # mod            id
ExpandQop:=proc(x)
    local ret,j,lkb,cbk,rkb,no,lop;
    lkb:=0; cbk:=0; rkb:=0;
    lop:=op(0,x);
    no:=nops(x);
    if lop = opsym or lop = opint or lop = opflt or
       lop = opbra or lop = opket or lop = opcpx then
         ret:=x;
    else
    if lop = opexp then
        ret:=x;
    else       
    if lop = opnp then
        ret:=1;
        for j from 1 to no do
            ret:=ret*ExpandQop(op(j,x));
        end do;        
    else
    if lop = `+` then
        ret:=0;
        for j from 1 to no do
            ret:=ret+ExpandQop(op(j,x));
        end do;
    else
    if lop = `-` then
        ret:=0;
        for j from 1 to no do
            ret:=ret-ExpandQop(op(j,x));
        end do;
    else
    if lop = opdiv then
       ret:=1;
       for j from 1 to no do
           ret:=ret/ExpandQop(op(j,x));
       end do;
    else
    if lop = opmod then
       ret:=x;
    else
    if lop = optp then
        ret:=1;
       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;
    else ret:=x;
    fi; # optp
    fi; # opmod
    fi; # opdiv
    fi; # `-`
    fi; # `+`
    fi; # `opnp
    fi; # `opexp`
    fi; # opsym, opint, opflt, opbra, opket, opcpx

    return ret;
end proc:

# For saving
# save opexp,opnp,optp,opdiv,opint,opflt,opcpx,opbra,opket,opmod, ExpandQop,"ExpandQop.m"

# 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 elements computed with C appears curious");
'bras3.C. kets3'="...";
bras3.C.kets3;
print("Matrix elements 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 elements computed with C appears curious"

 

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 elements 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)

print("Example");
En:=ExpandQop(10*(1-x+y+z)*i*(1/sqrt(2))*A*B);

"Example"

 

-(5*I)*2^(1/2)*(-1+x-y-z)*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))

(5)

print("Another example");
'F'='A*B/sqrt(2)+B*A/sqrt(2)';
F:=A*B/sqrt(2)+B*A/sqrt(2):
'op(1,F)'=op(1,F);
'op(2,F)'=op(2,F);

'Fn'='ExpandQop(F)';
Fn:=ExpandQop(F):
'op(1,Fn)'=op(1,Fn);
'op(2,Fn)'=op(2,Fn);

"Another example"

 

F = Physics:-`*`(Physics:-`*`(A, B), Physics:-`^`(sqrt(2), -1))+Physics:-`*`(Physics:-`*`(B, A), Physics:-`^`(sqrt(2), -1))

 

op(1, F) = (1/2)*2^(1/2)*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))

 

op(2, F) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), Physics:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1), 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))

 

Fn = ExpandQop(F)

 

op(1, Fn) = (1/2)*2^(1/2)*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))

 

op(2, Fn) = (1/2)*2^(1/2)*Physics:-`*`(Physics:-Ket(b1), Physics:-Ket(b2), Physics:-Ket(b3), 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:-Bra(b3), Physics:-Bra(b2), Physics:-Bra(b1))

(6)

 


 

Download ExpandQopV2.mw

The set and list produced by map (see below) contain duplicates.  How to remove duplicates?
 

p := (1+5^(1/2))*(1/2)

1/2+(1/2)*5^(1/2)

(1)

with(Bits)

[And, FirstNonzeroBit, GetBits, Iff, Implies, Join, Nand, Nor, Not, Or, Settings, Split, String, Xor]

(2)

with(LinearAlgebra)

[`&x`, Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CARE, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, CompanionMatrix, CompressedSparseForm, ConditionNumber, ConstantMatrix, ConstantVector, Copy, CreatePermutation, CrossProduct, DARE, DeleteColumn, DeleteRow, Determinant, Diagonal, DiagonalMatrix, Dimension, Dimensions, DotProduct, EigenConditionNumbers, Eigenvalues, Eigenvectors, Equal, ForwardSubstitute, FrobeniusForm, FromCompressedSparseForm, FromSplitForm, GaussianElimination, GenerateEquations, GenerateMatrix, Generic, GetResultDataType, GetResultShape, GivensRotationMatrix, GramSchmidt, HankelMatrix, HermiteForm, HermitianTranspose, HessenbergForm, HilbertMatrix, HouseholderMatrix, IdentityMatrix, IntersectionBasis, IsDefinite, IsOrthogonal, IsSimilar, IsUnitary, JordanBlockMatrix, JordanForm, KroneckerProduct, LA_Main, LUDecomposition, LeastSquares, LinearSolve, LyapunovSolve, Map, Map2, MatrixAdd, MatrixExponential, MatrixFunction, MatrixInverse, MatrixMatrixMultiply, MatrixNorm, MatrixPower, MatrixScalarMultiply, MatrixVectorMultiply, MinimalPolynomial, Minor, Modular, Multiply, NoUserValue, Norm, Normalize, NullSpace, OuterProductMatrix, Permanent, Pivot, PopovForm, ProjectionMatrix, QRDecomposition, RandomMatrix, RandomVector, Rank, RationalCanonicalForm, ReducedRowEchelonForm, Row, RowDimension, RowOperation, RowSpace, ScalarMatrix, ScalarMultiply, ScalarVector, SchurForm, SingularValues, SmithForm, SplitForm, StronglyConnectedBlocks, SubMatrix, SubVector, SumBasis, SylvesterMatrix, SylvesterSolve, ToeplitzMatrix, Trace, Transpose, TridiagonalForm, UnitVector, VandermondeMatrix, VectorAdd, VectorAngle, VectorMatrixMultiply, VectorNorm, VectorScalarMultiply, ZeroMatrix, ZeroVector, Zip]

(3)

with(VectorCalculus)

[`&x`, `*`, `+`, `-`, `.`, `<,>`, `<|>`, About, AddCoordinates, ArcLength, BasisFormat, Binormal, Compatibility, ConvertVector, CrossProduct, Curl, Curvature, D, Del, DirectionalDiff, Divergence, DotProduct, Flux, GetCoordinateParameters, GetCoordinates, GetNames, GetPVDescription, GetRootPoint, GetSpace, Gradient, Hessian, IsPositionVector, IsRootedVector, IsVectorField, Jacobian, Laplacian, LineInt, MapToBasis, Nabla, Norm, Normalize, PathInt, PlotPositionVector, PlotVector, PositionVector, PrincipalNormal, RadiusOfCurvature, RootedVector, ScalarPotential, SetCoordinateParameters, SetCoordinates, SpaceCurve, SurfaceInt, TNBFrame, Tangent, TangentLine, TangentPlane, TangentVector, Torsion, Vector, VectorField, VectorPotential, VectorSpace, Wronskian, diff, eval, evalVF, int, limit, series]

(4)

b := proc (x, w) options operator, arrow; (-1)^GetBits(x, w, output = number) end proc

proc (x, w) options operator, arrow; (-1)^Bits:-GetBits(x, w, output = number) end proc

(5)

l := proc (x, t, u, v) options operator, arrow; frac(x)*Vector([b(floor(x), 0)*t, b(floor(x), 1)*u, b(floor(x), 2)*v])+(1-frac(x))*Vector([b(floor(x), 0)*v, b(floor(x), 1)*t, b(floor(x), 2)*u]) end proc

proc (x, t, u, v) options operator, arrow; VectorCalculus:-`+`(VectorCalculus:-`*`(frac(x), VectorCalculus:-Vector([VectorCalculus:-`*`(b(floor(x), 0), t), VectorCalculus:-`*`(b(floor(x), 1), u), VectorCalculus:-`*`(b(floor(x), 2), v)])), VectorCalculus:-`*`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(frac(x))), VectorCalculus:-Vector([VectorCalculus:-`*`(b(floor(x), 0), v), VectorCalculus:-`*`(b(floor(x), 1), t), VectorCalculus:-`*`(b(floor(x), 2), u)]))) end proc

(6)

map(l, {0, 1, 2, 3, 4, 5, 6, 7, 8}, 0, 1, p)

{Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian])}

(7)

map(l, [0, 1, 2, 3, 4, 5, 6, 7, 8], 0, 1, p)

[Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian])]

(8)

q := ListTools:-MakeUnique(%)

q := [Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian])]

(9)

Equal(q[1], q[9])

true

(10)

qq := [op({q[]})]

qq := [Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = -1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = -1/2-(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian]), Vector(3, {(1) = 1/2+(1/2)*5^(1/2), (2) = 0, (3) = 1}, attributes = [coords = cartesian])]

(11)

LinearAlgebra:-Equal(qq[1], qq[5])

true

(12)

NULL


 

Download cp.mw

 we know that  maple evaluate an expression as a tree, for example if i have the expression:

f:=x^2*exp(3/2)*sin(Pi/3-1/x);

i want to get all the tree, is there a command to obtain this !!

I get the following message:

"Error, (in Bits:-GetBits) argument 1, the number, must be a nonnegative integer"

Presumably "argument 1" does not mean that the offending argument is 1 (the problem is that the offending argument "must be a nonnegative integer", but 1 is a nonnegative integer).  So how do I get maple to report what the offending argument is? 

Could you help me to write maple code for solving following matrix differential equations?

where

In here, the matrices M,C,K,P are as follows ( M,C,K are nxn matrices and  V,P are nx1 matrices) 

In here, l,P,A,rho, alpha,v,N,E,I are constants.

Thank you very much.

I writed the matrices in maple. You can find it in the below.

 The Code.mw

--- Prolog.ue ---

The best things in life come free of charge.

Happiness, love, and wisdom of expertise are first few that hit my mind.

As a business economist, I keep my eyes keenly open to opportunities for growth; such as Maple 2017 training session.

It was a Saturday afternoon in Waterloo, ON, this chilly Feburary which was blessed by snowstorm warning.

 

--- Encountering with Maple ---

I was aware of Maple for many years back when my academic career began.

In fact, Maple was available in the lab computers at university. 

But I did not know what to do with it.

Nor did I use any mathematics softwares until recently, but I had this thought : one day I could learn.

The motivation for this `learn how to use it' did not occur to me for a long time (14 years!!).

Things changed this year when I enrolled to an Electrical Engineering program at Lassonde.

Mind you, I have already been using various types of languages and tools such as: Python, C, Java, OpenOfficeSuites, Stata, SAS, Latex just to mention a few.

These stuffs also run on multiple platforms which I am sure you have heard of if you're reading this post; Windows, OSX and Linux. And Maple supports all these major operating systems.

 

--- Why do I like Maple ---

During the first week of school, Dr. Smith would ask us to purchase and practice using MATLAB because it had a relatively easy learning curve for beginners like python and we were going to use it for labs.

Furthermore, students get a huge discount (i.e. $1500 to $50) for these softwares.

Then, the professor also added; "Maple is also a great tool to use, but we won't use it for this class".

ME: ' Why not ? '

The curiosity inside me gave in and I decided to try both!

After all, my laziness in taking derivatives by hand or the possibility of making mistake would disappear if I can verify results using software.

That's it...!

Being able to check correct answer was already worth more than $50.

I can not emphasize this point enough; 

For people in the industry being paid for their time, or students like me who got a busy schedule can not afford to waste any time. (i.e. need to minimize homework effort & frustration, while maximizing the educational attainment & final grades)

Right? Time is money.

Don't we all just want more spare time for things we care?

Googling through many ambiguous Yahoo Answers or online forums like Stackoverflow replies are often misleading and time consuming. 

I have spent years (estimated 3000+ hours) going through those wildly inaccurate webpages hoping for some clearly written information with sub-optimal outcome.

Diverting many hours of study time is not something a first year S.T.E.M. students can afford.

 

--- Maple Training ---

Now you know about my relationships with Maple; Let me describe how the training session went.

I will begin with the sad news first, =(

First of all, there was no coffee available when I arrived. It arrived only after lunch.

Although it was a free event aside other best things in life, this was only a material factor I didn't enjoy at the site. 

Still a large portion of Canadians start their work with a zolt of caffeine in my defence.

Secondly, there was a kind of assumption which expected attendee were familiar with software behavior.

A handful of people were having trouble opening example file, perhaps because of their browser setting or link to preferred software by OS.

Not being able to follow the tutorials as the presenter demonstrated various facets of software substantially diminished the  efficacy of training session for those who could not be on the same page.

These minor annoyances were the only drawbacks I experinced from the event.

 

Here comes the happy side, =)

1. The staffs were considerate enough to provide vegetarion options for inclusive lunch as well as answering all my curious, at times orthogonal questions regarding Maplesoft company.

2. Highly respectable professionals were presenting themselves; 

That is, Prof. Illias Kotsireas, Dr. Erik Postma and Dr. Jürgen Gerhard.

I can not appreciate enough of their contribution for the training in an eloquent and humble manners.

To put it other way, leading of the presentation was well structured and planned out.

In the beginning, Prof. Kotsireas presented `Introduction to Maple' which included terminology and basic behaviors of Maple (i.e. commands and features) with simple examples you can quickly digest. Furthermore, Maple has internal function to interface with Latex! No more typing hours of $$s and many frac{}{}, \delta_{} to publish. In order for me to study all this would have been two-weeks kind of commitment in which he summarized in a couple of hours time. Short-cut keys that are often used by his project was pretty interesting, which will improve work efficiency.

After a brief lunch, which was supplied more than enough for all, Dr. Erik Postma delivered a critical component of simluation. That is, `Random Number Generation'. Again, he showed us some software-related tricks such as `Text mode' vs. `Math mode'.  The default RNG embedded in the software allows reproducible results unless we set seed and randomize further. Main part of the presentation was regarding `Optimization of solution through simulation'. He iteratively improved efficiency of test model, which I will not go in depth here. However, visually and quantitatively showing the output was engaging the attendees and Maple has modularized this process (method available for all the users!!).

Finally, we got some coffee break that allowed to me to push through all the way to the end. I believe if we had some coffee earlier less attendees would have left.

The last part of the training was presented by Dr. Jürgen Gerhard. In this part, we were using various applications of Maple in solving different types of problems. We tackled combinatorics of Fibonacci sequence by formula manipulation. In this particular example he showed us how to optimize logic of a function that made a huge impact in processing time and memory usage. Followed by graph theory example, damped harmonic oscillator, 2 DOF chaotic system, optimization and lastly proof of orthocentre by coding. I will save the examples for you to enjoy in future sessions. 

The way they went through examples were super easy to follow. This can only be done with profound understanding of the subject and a lot of prior effort in preparing the presentation.
 

I appreciate much efforts put together by whom organized this event, allocating their own precious weekend time and allowing many to gain opportunity to learn directly from the person in the house.

 

--- Epilogue ---

My hope for Maple usage lies in enhancing education outcome for first year students, especially in the field of Science and Economics. This is a free opportunity for economic empowerment which is uncaptured.

Engineering students are already pretty good at problem solving, and will figure things out as I witnessed my colleagues have.

However, students of natural sciences and B.A. programs tend to skimp on utilizing tools due to lack of exposure.

Furthermore, I am supporting their development of SaaS, software as service, which delivers modules like gRPC does.

Also, I hope the optimization package from prior version written by Dr. Postma will become available to public sometime.

Here's a BIG thank you to staffs once again, and forgive me for any grammatical errors from rushed writing. I tried to incorporate as much observation as possible gathered from the event.

To contact me, my email is hyonwoo.kee (at) gmail.com;

I am looking to develop a histogram from a random data set that contains integers that are approximately normally distributed.  The intent is to develop a histogram from the dataset that looks like a fairly decent-approximation of a normalized curve (but not necessarily perfect)....  the next stage to this is a second histogram that shows a bimodal distribution (but not necessarily perfect)

 

To this point I have...

restart;

with(Statistics);

randomize();

N := RandomVariable(Normal(300, 10));

A := [seq(MapleTA:-Builtin:-decimal(0, Sample(N, 500)[i]), i = 1 .. 500)];

Histogram(A, frequencyscale = absolute, bincount = 25, binwidth = 1, tickmarks = [default, default]);

 

I am certain you will find what I have thus far is far from ideal, as such any help is appreciated.

 

Thanks,

Mark

 

To restrict the domain of a vector field, I have multiplied a coordinate with a non-real complex number (namely a sqrt(of negative expression)). This does work, as shown in this Maple 2017 worksheet program (below). My question is whether this is the best technique of accomplishing this result, or else how to do it better? Would be interested in suggestions for improvements. Here is my program so far:

restart;
#
with(plots):
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(` Vector Field inside Torus`);
print(` ------- ------- ------- ------- ------- ------- -------`);
print(` Assignment:  `);
print(` In a circular pipe of radius (my2r), water is flowing in the direction `);
print(` of the pipe, with speed (my2r)^2-(mya)^2, where (mya) is the distance  `);
print(` to the axis of the pipe.  `);
print(` Depict the vector field describing the flow if the pipe goes around in `);
print(` the shape of a torus with major radius (my1r).  `);
print(`   `);
print(`   `);
print(`   `);
print(` ------- ------- ------- ------- ------- ------- -------`);
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 1) major radius of torus:`);
#
my1r  := 5;     
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 2) minor radius of torus (pipe radius):`);
#
my2r := 4; 
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 3) definition of torus (polar coordinates):`);
#
c00 := [(my1r+my2r*cos(s))*cos(t),(my1r+my2r*cos(s))*sin(t),my2r*sin(s)];
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 4) 3D plot of solid torus (polar coordinates):`);
#
plot3d({c00},scaling=constrained,color=red);
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 5) 3D plot of wireframe torus (polar coordinates):`);
#
P1 := plot3d({c00},scaling=constrained,style=wireframe);
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 6) implicit definition of torus (cartesian coordinates):`);
#
c01 := (sqrt(x^2+y^2)-my1r)^2+z^2-my2r^2;
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 7) implicit 3D plot of solid torus (cartesian coordinates):`);
#
gx := my1r+my2r; # min and max of each coordinate
#
implicitplot3d(c01,x=-gx..gx,y=-gx..gx,z=-gx..gx,numpoints=9000);
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 8) vector field definition (cartesian coordinates):`);
#
my1vfx := -y;
my1vfy := x;
my1vfz := 0;
#
my1fld := [my1vfx,my1vfy,my1vfz];
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 9) 3D plot of vector field (cartesian coordinates):`);

#
fieldplot3d(my1fld,x=-gx..gx,y=-gx..gx,z=-gx..gx);
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 10) definition of vector field with unit length vectors (cartesian coordinates):`);
#
my1vl := sqrt(my1vfx^2+my1vfy^2+my1vfz^2); # vector length
#
my2fld := [my1vfx/my1vl,my1vfy/my1vl,my1vfz/my1vl];
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 11) 3D plot of vector field with unit length vectors (cartesian coordinates):`);
#
fieldplot3d(my2fld,x=-gx..gx,y=-gx..gx,z=-gx..gx);
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 12) definition of vector field with asked for length vectors (cartesian coordinates):`);
#
mya := sqrt((sqrt(x^2+y^2)-my1r)^2+z^2);
c01r := sqrt(my2r^2-mya^2); # also used for domain restricting vector field below
#
my1tsz := solve([c01],[z]);
#
assign(my1tsz[1][1]);
my1tz := z;
unassign('z');
#
assign(my1tsz[2][1]);
my2tz := z;
unassign('z');
#
my1vp := c01r/my2r; # vector length (maximum one unit)
#
my3fld := [my1vp*my1vfx/my1vl,my1vp*my1vfy/my1vl,my1vp*my1vfz/my1vl];
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 13) 3D plot of vector field with asked for length vectors (cartesian coordinates):`);
#
fieldplot3d(my3fld,x=-gx..gx,y=-gx..gx,z=-gx..gx);
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n(Section 14) same asked for vector field with 3-D arrow vectors:`);
print(`   `);
print(` (to get this to display properly it was necessary to do:  `);
print(` -> Maple 2017 -> Preferences... -> Precision ->   `);
print(`   [unselect] Limit expression length to   `);
print(`   Apply to Session`);
print(`   `);
#
gr := 15;
#
P3 := fieldplot3d(my3fld,x=-gx..gx,y=-gx..gx,z=-gx..gx,arrows=`3-D`,grid=[gr,gr,gr]);
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);
print(`\n display asked for vector field within wireframe torus:`);
#
display([P1,P3]);
#
#
print(` ------- ------- ------- ------- ------- ------- -------`);

 

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

Hi

However I run similar codes in Maple and Matlab, but the different results are observed.

By increasing number of basis (i.e. m) The convergence is observed, but results are different from each other.

The convergence for Maple: (m,P) =(10,10.8154),(20,10.8081),(30,10.8067),(50,10.8061)

The convergence for Matlab: (m,P) =(10,10.0474),(20,10.0194),(30,10.0143),(40,10.0125)

I have two questions:

1-which answer is correct?

2-why the different results for the same  number of basis  are obtained?

 

Maple code (Maple 2016)

restart;

tm := time():

with(LinearAlgebra):

Digits := 500:

beta := 1:

nu := 0.3:

lambda := 2:

G := 5:

ko := .5*0.1e7:

Ec := 0.380e12:

Em := 0.70e11:

ri := 0:

ro := 0.5:

ti := 0.1e-1:

K := 0.4e-1*0.1e10:

n := 1:

m := 20:

alpha := 0.1:

t := ti*(1+alpha*r/ro)^beta:

Er := (Ec-Em)*(r/ro)^n+Em:

W := simplify(add(a[n]*ChebyshevT(n, r), n = 0 .. m)):

sys := {eval(W, r = ro), eval(diff(W, r), r = ri)}:

W := subs(solve(sys, {a[0], a[1]}), W):

d1 := diff(W, r):

d2 := diff(d1, r):

W := subs(solve({eval(ko*d1-Ec*ti^3*(1+alpha)^(3*beta)*(d1*nu+d2*r)/(12*r*(-nu^2+1)), r = ro) = 0}, {a[2]}), W):

d1 := diff(W, r):

d2 := diff(d1, r):

Uf := int(K*(1+lambda*r/ro+G*(r/ro)^2)*W^2*r, r = ri .. ro):

NUM := int(((d2+d1/r)^2-(2*(1-nu))*d2*d1/r)*Er*t^3*r/(12*(-nu^2+1)), r = ri .. ro)+ko*ro*(eval(d1, r = ro))^2+Uf:

DEN := int(d1^2*r, r = ri .. ro):

PT := NUM-Em*ti^3*P*DEN/ro^2:

for c from 3 to m do

eq[c] := diff(PT, a[c]) end do:

MT := GenerateMatrix([seq(eq[j], j = 3 .. m)], [a[j]$j = 3 .. m])[1]:

solve(Determinant(MT)):

evalf[10](%[1]);

____________________________________________________________________________

Matlab Code (R2014a)

clear all;
digits(10);
tic
m = 20;
n = 1;
alpha = 0.1;
beta = 1;
nu = 0.3;
lambda = 2;
G = 5;
ko = 0.5*1e6;
Ec = 0.380e12;
Em = 0.70e11;
ro = 0.5;
ti = 0.1e-1;
K = 0.04*0.1e10;
ri=0;
syms r P;
c=sym('c%d',[1,m+1]);
c(1)=1;
c(2)=r;
for j=3:m+1
    c(j)=simplify(2*r*c(j-1)-c(j-2));
end
Er = (Ec-Em)*(r/ro)^n+Em;
a=sym('a%d',[1,m+1]);
t = ti*(1+alpha*r/ro)^beta;
W = simplify(a*transpose(c));
[a(1),a(2)] = solve(subs(W,r,ro)==0,subs(diff(W,r),r,ri)==0,a(1),a(2));
W = simplify(a*transpose(c));
d1 = diff(W, r);
d2 = diff(d1, r);
a(3)=solve(subs(ko*d1-Ec*ti^3*(1+alpha)^(3*beta)*(d1*nu+d2*r)/(12*r*(-nu^2+1)),r,ro)==0,a(3));
W = simplify(a*transpose(c));
d1 = diff(W, r);
d2 = diff(d1, r);
Uf = int(K*(1+lambda*r/ro+G*(r/ro)^2)*W^2*r,r,ri,ro);
NUM = int(((d2+d1/r)^2-(2*(1-nu))*d2*d1/r)*Er*t^3*r/(12*(-nu^2+1)),r,ri,ro)+ko*ro*(subs(d1,r,ro))^2+Uf;
DEN = int(d1^2*r, r , ri , ro);
PT = NUM-Em*ti^3*P*DEN/ro^2;
for q=1:m-2
    eq(q)= diff(PT, a(q+3));
end;
for q = 1:m-2
    for u = 1:m-2
        T = coeffs(eq(q),a(u+3));
        tt(q,u) = T(2);
    end;
end;
F=sort(vpa(solve(det(tt))));
F(1)
toc

___________________________________________________________

First 840 841 842 843 844 845 846 Last Page 842 of 2216