ecterrab

14837 Reputation

24 Badges

21 years, 88 days

MaplePrimes Activity


These are Posts that have been published by ecterrab


 

Quantum Commutation Rules Basics

 

Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

(1) Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France

(2) Maplesoft

NULL

NULL

In Quantum Mechanics, in the coordinates representation, the component of the momentum operator along the x axis is given by the differential operator


 "`p__x`=-i `ℏ`(∂)/(∂x)  "

 

The purpose of the exercises below is thus to derive the commutation rules, in the coordinates representation, between an arbitrary function of the coordinates and the related momentum, departing from the differential representation

 

p[n] = -i*`ℏ`*`∂`[n]

These two exercises illustrate how to have full control of the computational process by using different elements of the Maple language, including inert representations of abstract vectorial differential operators, Hermitian operators, algebra rules, etc.

 

These exercises also illustrate a new feature of the Physics package, introduced in Maple 2017, that is getting refined (the computation below requires the Maplesoft updates of the Physics package) which is the ability to perform computations algebraically, using the product operator, but with differential operators, and transform the products into the application of the operators only when we want that, as we do with paper and pencil.

 

%Commutator(g(x, y, z), p_) = I*`ℏ`*Nabla(F(X))

 
> 

restart; with(Physics); with(Physics[Vectors]); interface(imaginaryunit = i)

 

Start setting the problem:

– 

 all ofx, y, z, p__x, p__y, p__z are Hermitian operators

– 

 all of x, y, z commute between each other

– 

 tell the system only that the operators x, y, z are the differentiation variables of the corresponding (differential) operators p__x, p__y, p__z but do not tell what is the form of the operators

 

> 

Setup(mathematicalnotation = true, differentialoperators = {[p_, [x, y, z]]}, hermitianoperators = {p, x, y, z}, algebrarules = {%Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, z) = 0}, quiet)

[algebrarules = {%Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, z) = 0}, differentialoperators = {[p_, [x, y, z]]}, hermitianoperators = {p, x, y, z}, mathematicalnotation = true]

(1.1)

Assuming F(X) is a smooth function, the idea is to apply the commutator %Commutator(F(X), p_) to an arbitrary ket of the Hilbert space Ket(psi, x, y, z), perform the operation explicitly after setting a differential operator representation for `#mover(mi("p",mathcolor = "olive"),mo("→"))`, and from there get the commutation rule between F(X) and `#mover(mi("p",mathcolor = "olive"),mo("→"))`.

 

Start introducing the commutator, to proceed with full control of the operations we use the inert form %Commutator

> 

alias(X = (x, y, z))

> 

CompactDisplay(F(X))

` F`(X)*`will now be displayed as`*F

(1.2)
> 

%Commutator(F(X), p_)*Ket(psi, X)

Physics:-`*`(%Commutator(F(X), p_), Physics:-Ket(psi, x, y, z))

(1.3)

For illustration purposes only (not necessary), expand this commutator

> 

Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)) = expand(Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)))

Physics:-`*`(%Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = Physics:-`*`(F(X), p_, Physics:-Ket(psi, x, y, z))-Physics:-`*`(p_, F(X), Physics:-Ket(psi, x, y, z))

(1.4)

Note that  `#mover(mi("p",mathcolor = "olive"),mo("→"))`, F(X) and the ket Ket(psi, x, y, z) are operands in the products above and that they do not commute: we indicated that the coordinates x, y, z are the differentiation variables of `#mover(mi("p",mathcolor = "olive"),mo("→"))`. This emulates what we do when computing with these operators with paper and pencil, where we represent the application of a differential operator as a product operation.

 

This representation can be transformed into the (traditional in computer algebra) application of the differential operator when desired, as follows:

> 

Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)) = Library:-ApplyProductsOfDifferentialOperators(Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)))

Physics:-`*`(%Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = Physics:-`*`(F(X), p_(Physics:-Ket(psi, x, y, z)))-p_(Physics:-`*`(F(X), Physics:-Ket(psi, x, y, z)))

(1.5)

Note that, in `#mover(mi("p",mathcolor = "olive"),mo("→"))`(F(X)*Ket(psi, x, y, z)), the application of `#mover(mi("p",mathcolor = "olive"),mo("→"))` is not expanded: at this point nothing is known about  `#mover(mi("p",mathcolor = "olive"),mo("→"))` , it is not necessarily a linear operator. In the Quantum Mechanics problem at hands, however, it is. So give now the operator  `#mover(mi("p",mathcolor = "olive"),mo("→"))` an explicit representation as a linear vectorial differential operator (we use the inert form %Nabla, %Nabla, to be able to proceed with full control one step at a time)

> 

p_ := proc (f) options operator, arrow; -I*`ℏ`*%Nabla(f) end proc

proc (f) options operator, arrow; -Physics:-`*`(Physics:-`*`(I, `ℏ`), %Nabla(f)) end proc

(1.6)

The expression (1.5) becomes

> 

Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)) = Physics[`*`](F(X), p_(Physics[Ket](psi, x, y, z)))-p_(Physics[`*`](F(X), Physics[Ket](psi, x, y, z)))

Physics:-`*`(%Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = -I*`ℏ`*Physics:-`*`(F(X), %Nabla(Physics:-Ket(psi, x, y, z)))+I*`ℏ`*%Nabla(Physics:-`*`(F(X), Physics:-Ket(psi, x, y, z)))

(1.7)

Activate now the inert operator VectorCalculus[Nabla] and simplify taking into account the algebra rules for the coordinate operators {%Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, z) = 0}

> 

Simplify(value(Physics[`*`](%Commutator(F(X), p_), Physics[Ket](psi, x, y, z)) = -I*`ℏ`*Physics[`*`](F(X), %Nabla(Physics[Ket](psi, x, y, z)))+I*`ℏ`*%Nabla(Physics[`*`](F(X), Physics[Ket](psi, x, y, z)))))

Physics:-`*`(Physics:-Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = I*`ℏ`*_i*Physics:-`*`(diff(F(X), x), Physics:-Ket(psi, x, y, z))+I*`ℏ`*_j*Physics:-`*`(diff(F(X), y), Physics:-Ket(psi, x, y, z))+I*`ℏ`*_k*Physics:-`*`(diff(F(X), z), Physics:-Ket(psi, x, y, z))

(1.8)

To make explicit the gradient in disguise on the right-hand side, factor out the arbitrary ket Ket(psi, x, y, z)

> 

Factor(Physics[`*`](Physics[Commutator](F(X), p_), Physics[Ket](psi, x, y, z)) = I*`ℏ`*_i*Physics[`*`](diff(F(X), x), Physics[Ket](psi, x, y, z))+I*`ℏ`*_j*Physics[`*`](diff(F(X), y), Physics[Ket](psi, x, y, z))+I*`ℏ`*_k*Physics[`*`](diff(F(X), z), Physics[Ket](psi, x, y, z)))

Physics:-`*`(Physics:-Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = I*`ℏ`*Physics:-`*`((diff(F(X), y))*_j+(diff(F(X), z))*_k+(diff(F(X), x))*_i, Physics:-Ket(psi, x, y, z))

(1.9)

Combine now the expanded gradient into its inert (not-expanded) form

> 

subs((Gradient = %Gradient)(F(X)), Physics[`*`](Physics[Commutator](F(X), p_), Physics[Ket](psi, x, y, z)) = I*`ℏ`*Physics[`*`]((diff(F(X), y))*_j+(diff(F(X), z))*_k+(diff(F(X), x))*_i, Physics[Ket](psi, x, y, z)))

Physics:-`*`(Physics:-Commutator(F(X), p_), Physics:-Ket(psi, x, y, z)) = I*`ℏ`*Physics:-`*`(%Gradient(F(X)), Physics:-Ket(psi, x, y, z))

(1.10)

Since (1.10) is true for allKet(psi, x, y, z), this ket can be removed from both sides of the equation. One can do that either taking coefficients (see Coefficients ) or multiplying by the "formal inverse" of this ket, arriving at the (expected) form of the commutation rule between F(X) and `#mover(mi("p",mathcolor = "olive"),mo("→"))`

> 

(Physics[`*`](Physics[Commutator](F(X), p_), Ket(psi, x, y, z)) = I*`ℏ`*Physics[`*`](%Gradient(F(X)), Ket(psi, x, y, z)))*Inverse(Ket(psi, x, y, z))

Physics:-Commutator(F(X), p_) = I*`ℏ`*%Gradient(F(X))

(1.11)

Tensor notation, "[`X__m`,P[n]][-]=i `ℏ` g[m,n]"

 

The computation rule for position and momentum, this time in tensor notation, is performed in the same way, just that, additionally, specify that the space indices to be used are lowercase latin letters, and set the relationship between the differential operators and the coordinates directly using tensor notation.

You can also specify that the metric is Euclidean, but that is not necessary: the default metric of the Physics package, a Minkowski spacetime, includes a 3D subspace that is Euclidean, and the default signature, (- - - +), is not a problem regarding this computation.

 

> 

restart; with(Physics); interface(imaginaryunit = i)

> 

Setup(mathematicalnotation = true, coordinates = cartesian, spaceindices = lowercaselatin, algebrarules = {%Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, z) = 0}, hermitianoperators = {P, X, p}, differentialoperators = {[P[m], [x, y, z]]}, quiet)

[algebrarules = {%Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, z) = 0}, coordinatesystems = {X}, differentialoperators = {[P[m], [x, y, z]]}, hermitianoperators = {P, p, t, x, y, z}, mathematicalnotation = true, spaceindices = lowercaselatin]

(2.1)

Define now the tensor P[m]

> 

Define(P[m], quiet)

{Physics:-Dgamma[mu], P[m], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[a, b], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(2.2)

Introduce now the Commutator, this time in active form, to show how to reobtain the non-expanded form at the end by resorting the operands in products

> 

Commutator(X[m], P[n])*Ket(psi, x, y, z)

Physics:-`*`(Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]), Physics:-Ket(psi, x, y, z))

(2.3)

Expand first (not necessary) to see how the operator P[n] is going to be applied

> 

Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)) = expand(Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)))

Physics:-`*`(Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]), Physics:-Ket(psi, x, y, z)) = Physics:-`*`(Physics:-SpaceTimeVector[m](X), P[n], Physics:-Ket(psi, x, y, z))-Physics:-`*`(P[n], Physics:-SpaceTimeVector[m](X), Physics:-Ket(psi, x, y, z))

(2.4)

Now expand and directly apply in one ago the differential operator P[n]

> 

Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)) = Library:-ApplyProductsOfDifferentialOperators(Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)))

Physics:-`*`(Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]), Physics:-Ket(psi, x, y, z)) = Physics:-`*`(Physics:-SpaceTimeVector[m](X), P[n](Physics:-Ket(psi, x, y, z)))-P[n](Physics:-`*`(Physics:-SpaceTimeVector[m](X), Physics:-Ket(psi, x, y, z)))

(2.5)

Introducing the explicit differential operator representation for P[n], here again using the inert %d_[n] to keep control of the computations step by step

> 

P[n] := proc (f) options operator, arrow; -I*`ℏ`*%d_[n](f) end proc

proc (f) options operator, arrow; -Physics:-`*`(Physics:-`*`(I, `ℏ`), %d_[n](f)) end proc

(2.6)

The expanded and applied commutator (2.5) becomes

> 

Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)) = Physics[`*`](Physics[SpaceTimeVector][m](X), P[n](Ket(psi, x, y, z)))-P[n](Physics[`*`](Physics[SpaceTimeVector][m](X), Ket(psi, x, y, z)))

Physics:-`*`(Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]), Physics:-Ket(psi, x, y, z)) = -I*`ℏ`*Physics:-`*`(Physics:-SpaceTimeVector[m](X), %d_[n](Physics:-Ket(psi, x, y, z)))+I*`ℏ`*%d_[n](Physics:-`*`(Physics:-SpaceTimeVector[m](X), Physics:-Ket(psi, x, y, z)))

(2.7)

Activate now the inert operators %d_[n] and simplify taking into account Einstein's rule for repeated indices

> 

Simplify(value(Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)) = -I*`ℏ`*Physics[`*`](Physics[SpaceTimeVector][m](X), %d_[n](Ket(psi, x, y, z)))+I*`ℏ`*%d_[n](Physics[`*`](Physics[SpaceTimeVector][m](X), Ket(psi, x, y, z)))))

Physics:-`*`(Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]), Physics:-Ket(psi, x, y, z)) = I*`ℏ`*Physics:-g_[m, n]*Physics:-Ket(psi, x, y, z)

(2.8)

Since the ket Ket(psi, x, y, z) is arbitrary, we can take coefficients (or multiply by the formal Inverse  of this ket as done in the previous section). For illustration purposes, we use   Coefficients  and note hwo it automatically expands the commutator

> 

Coefficients(Physics[`*`](Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]), Ket(psi, x, y, z)) = I*`ℏ`*Physics[g_][m, n]*Ket(psi, x, y, z), Ket(psi, x, y, z))

Physics:-`*`(Physics:-SpaceTimeVector[m](X), P[n])-Physics:-`*`(P[n], Physics:-SpaceTimeVector[m](X)) = I*`ℏ`*Physics:-g_[m, n]

(2.9)

One can undo this (frequently undesired) expansion of the commutator by sorting the products on the left-hand side using the commutator between X[m] and P[n]

> 

Library:-SortProducts(Physics[`*`](Physics[SpaceTimeVector][m](X), P[n])-Physics[`*`](P[n], Physics[SpaceTimeVector][m](X)) = I*`ℏ`*Physics[g_][m, n], [P[n], X[m]], usecommutator)

Physics:-Commutator(Physics:-SpaceTimeVector[m](X), P[n]) = I*`ℏ`*Physics:-g_[m, n]

(2.10)

And that is the result we wanted to compute.

 

Additionally, to see this rule in matrix form,

> 

TensorArray(-(Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]) = I*`ℏ`*Physics[g_][m, n]))

Matrix(%id = 18446744078261558678)

(2.11)

In the above, we use equation (2.10) multiplied by -1 to avoid a minus sign in all the elements of (2.11), due to having worked with the default signature (- - - +); this minus sign is not necessary if in the Setup at the beginning one also sets  signature = `+ + + -`

 

For display purposes, to see this matrix expressed in terms of the geometrical components of the momentum `#mover(mi("p",mathcolor = "olive"),mo("→"))` , redefine the tensor P[n] explicitly indicating its Cartesian components

> 

Define(P[m] = [p__x, p__y, p__z], quiet)

{Physics:-Dgamma[mu], P[m], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[a, b], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(2.12)
> 

TensorArray(-(Physics[Commutator](Physics[SpaceTimeVector][m](X), P[n]) = I*`ℏ`*Physics[g_][m, n]))

Matrix(%id = 18446744078575996430)

(2.13)

Finally, in a typical situation, these commutation rules are to be taken into account in further computations, and for that purpose they can be added to the setup via

> 

"Setup(?)"

[algebrarules = {%Commutator(x, p__x) = I*`ℏ`, %Commutator(x, p__y) = 0, %Commutator(x, p__z) = 0, %Commutator(x, y) = 0, %Commutator(x, z) = 0, %Commutator(y, p__x) = 0, %Commutator(y, p__y) = I*`ℏ`, %Commutator(y, p__z) = 0, %Commutator(y, z) = 0, %Commutator(z, p__x) = 0, %Commutator(z, p__y) = 0, %Commutator(z, p__z) = I*`ℏ`}]

(2.14)

For example, from herein computations are performed taking into account that

> 

(%Commutator = Commutator)(x, p__x)

%Commutator(x, p__x) = I*`ℏ`

(2.15)
> 

NULL

NULL


 

Download DifferentialOperatorCommutatorRules.mw

 

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

This presentation is about magnetic traps for neutral particles, first achieved for cold neutrons and nowadays widely used in cold-atom physics. The level is that of undergraduate electrodynamics and tensor calculus courses. Tackling this topic within a computer algebra worksheet as shown below illustrates well the kind of advanced computations that can be done today with the Physics package. A new feature minimizetensorcomponents and related functionality is used along the presentation, that requires the updated Physics library distributed at the Maplesoft R&D Physics webpage.
 

 

Magnetic traps in cold-atom physics

 

Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

(1) Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France

(2) Maplesoft

 

We consider a device constructed with a set of electrical wires fed with constant electrical currents. Those wires can have an arbitrary complex shape. The device is operated in a regime such that, in some region of interest, the moving particles experience a magnetic field that varies slowly compared to the Larmor spin precession frequency. In this region, the effective potential is proportional to the modulus of the field: LinearAlgebra[Norm](`#mover(mi("B"),mo("→"))`(x, y, z)), this potential has a minimum and, close to this minimum, the device behaves as a magnetic trap.

 

 

 

Figure 1: Schematic representation of a Ioffe-Pritchard magnetic trap. It is made of four infinite rods and two coils.

_________________________________________

 

Following [1], we show that:

 

  

a) For a time-independent magnetic field  `#mover(mi("B"),mo("→"))`(x, y, z) in vacuum, up to order two in the relative coordinates X__i = [x, y, z] around some point of interest, the coefficients of orders 1 and 2 in this expansion, `v__i,j` and `c__i,j,k` , respectively the gradient and curvature, contain only 5 and 7 independent components.

  

b) All stationary points of LinearAlgebra[Norm](`#mover(mi("B"),mo("→"))`(x, y, z))^2 (nonzero minima and saddle points) are confined to a curved surface defined by det(`∂`[j](B[i])) = 0.

  

c) The effective potential, proportional to LinearAlgebra[Norm](`#mover(mi("B"),mo("→"))`(x, y, z)), has no maximum, only a minimum.

 

Finally, we draw the stationary condition surface for the case of the widely used Ioffe-Pritchard magnetic trap.

  

 

  

Reference

  

[1] R. Gerritsma and R. J. C. Spreeuw, Topological constraints on magnetostatic traps,  Phys. Rev. A 74, 043405 (2006)

  

 

The independent components of `v__i,j` and `c__i,j,k` entering B[i] = u[i]+v[i, j]*X[j]+(1/2)*c[i, j, k]*X[j]*X[k]

   

The stationary points are within the surface det(`∂`[j](B[i])) = 0

   

U = LinearAlgebra[Norm](`#mover(mi("B",fontweight = "bold"),mo("→",fontweight = "bold"))`)^2 has only minima, no maxima

   

Drawing the Ioffe-Pritchard Magnetic Trap

   


 

MagneticTraps.mw or in pdf format with the sections open: MagneticTraps.pdf

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


 

Maple 2017.1 and 2017.2 introduced several improvements in the solution of PDE & Boundary conditions problems (exact solutions).  Maple 2017.3 includes more improvements in this same area.

 

The following is a set of 25 examples of different PDE & Boundary Conditions problems that are solvable in Maple 2017.3 but not in Maple 2017.2 or previous releases. In the examples that follow, in some cases the PDE is different, in other cases the boundary conditions are of a different kind or solving the problem involves different computational strategies.

As usual, at the end there is a link pointing to the corresponding worksheet.

 

> 

pde[1] := diff(u(x, t), t) = k*(diff(u(x, t), x, x)); bc[1] := u(x, 0) = 6+4*cos(3*Pi*x/L), (D[1](u))(0, t) = 0, (D[1](u))(L, t) = 0

diff(u(x, t), t) = k*(diff(diff(u(x, t), x), x))

 

u(x, 0) = 6+4*cos(3*Pi*x/L), (D[1](u))(0, t) = 0, (D[1](u))(L, t) = 0

(1)
> 

`assuming`([pdsolve([pde[1], bc[1]], u(x, t))], [L > 0, k > 0])

u(x, t) = 6+4*cos(3*Pi*x/L)*exp(-9*k*Pi^2*t/L^2)

(2)
> 

pde[2] := diff(g(t, x), t) = diff(g(t, x), x, x)+a*g(t, x); bc[2] := g(0, x) = 1

diff(g(t, x), t) = diff(diff(g(t, x), x), x)+a*g(t, x)

 

g(0, x) = 1

(3)
> 

pdsolve([pde[2], bc[2]])

g(t, x) = exp(a*t)

(4)
> 

pde[3] := diff(u(x, t), t) = diff(u(x, t), x, x)

> 

bc[3] := u(x, 0) = f(x), u(-1, t) = 0, u(1, t) = 0

diff(u(x, t), t) = diff(diff(u(x, t), x), x)

 

u(x, 0) = f(x), u(-1, t) = 0, u(1, t) = 0

(5)
> 

pdsolve([pde[3], bc[3]], u(x, t))

u(x, t) = Sum((Int(f(x)*sin(n*Pi*x), x = -1 .. 1))*sin(n*Pi*x)*exp(-Pi^2*n^2*t), n = 1 .. infinity)

(6)
> 

pde[4] := diff(u(x, t), t) = k*(diff(u(x, t), x, x)); bc[4] := u(0, t) = 0, u(L, t) = 0, u(x, 0) = piecewise(0 < x and x <= (1/2)*L, 1, (1/2)*L < x and x < L, 2)

pde[4] := diff(u(x, t), t) = k*(diff(u(x, t), x, x))

 

u(0, t) = 0, u(L, t) = 0, u(x, 0) = piecewise(0 < x and x <= (1/2)*L, 1, (1/2)*L < x and x < L, 2)

(7)
> 

`assuming`([pdsolve([pde[4], bc[4]], u(x, t))], [L > 0])

u(x, t) = Sum((2*cos((1/2)*Pi*n)+2+4*(-1)^(1+n))*sin(n*Pi*x/L)*exp(-k*Pi^2*n^2*t/L^2)/(Pi*n), n = 1 .. infinity)

(8)
> 

pde[5] := diff(u(x, t), t) = k*(diff(u(x, t), x, x))

> 

bc[5] := (D[1](u))(0, t) = 0, (D[1](u))(L, t) = 0, u(x, 0) = -3*cos(8*Pi*x/L)

diff(u(x, t), t) = k*(diff(diff(u(x, t), x), x))

 

(D[1](u))(0, t) = 0, (D[1](u))(L, t) = 0, u(x, 0) = -3*cos(8*Pi*x/L)

(9)
> 

`assuming`([pdsolve([pde[5], bc[5]], u(x, t))], [0 < L, 0 < k])

u(x, t) = -3*cos(8*Pi*x/L)*exp(-64*k*Pi^2*t/L^2)

(10)
> 

pde[6] := diff(u(x, t), t) = k*(diff(u(x, t), x, x))+f(x, t); bc[6] := u(0, t) = 0, u(l, t) = 0, u(x, 0) = g(x)

diff(u(x, t), t) = k*(diff(diff(u(x, t), x), x))+f(x, t)

 

u(0, t) = 0, u(l, t) = 0, u(x, 0) = g(x)

(11)
> 

pdsolve([pde[6], bc[6]], u(x, t))

u(x, t) = Sum(2*(Int(g(tau1)*sin(Pi*n1*tau1/l), tau1 = 0 .. l))*sin(Pi*n1*x/l)*exp(-k*Pi^2*n1^2*t/l^2)/l, n1 = 1 .. infinity)+Int(Sum(2*(Int(f(x, tau1)*sin(Pi*n*x/l), x = 0 .. l))*sin(Pi*n*x/l)*exp(-k*Pi^2*n^2*(t-tau1)/l^2)/l, n = 1 .. infinity), tau1 = 0 .. t)

(12)
> 

pde[7] := diff(u(x, t), t) = diff(u(x, t), x, x); bc[7] := u(x, 0) = f(x), u(-1, t) = 0, u(1, t) = 0

diff(u(x, t), t) = diff(diff(u(x, t), x), x)

 

u(x, 0) = f(x), u(-1, t) = 0, u(1, t) = 0

(13)
> 

pdsolve([pde[7], bc[7]], u(x, t))

u(x, t) = Sum((Int(f(x)*sin(n*Pi*x), x = -1 .. 1))*sin(n*Pi*x)*exp(-Pi^2*n^2*t), n = 1 .. infinity)

(14)
> 

pde[8] := diff(u(x, t), t) = k*(diff(u(x, t), x, x))-h*u(x, t); bc[8] := u(x, 0) = sin(x), u(-Pi, t) = u(Pi, t), (D[1](u))(-Pi, t) = (D[1](u))(Pi, t)

diff(u(x, t), t) = k*(diff(diff(u(x, t), x), x))-h*u(x, t)

 

u(x, 0) = sin(x), u(-Pi, t) = u(Pi, t), (D[1](u))(-Pi, t) = (D[1](u))(Pi, t)

(15)
> 

pdsolve([pde[8], bc[8]], u(x, t))

u(x, t) = sin(x)*exp(-t*(k+h))

(16)
> 

pde[9] := diff(u(x, t), t) = diff(u(x, t), x, x)

> 

bc[9] := u(0, t) = 20, u(1, t) = 50, u(x, 0) = 0

diff(u(x, t), t) = diff(diff(u(x, t), x), x)

 

u(0, t) = 20, u(1, t) = 50, u(x, 0) = 0

(17)
> 

pdsolve([pde[9], bc[9]], u(x, t))

u(x, t) = 20+Sum((-40+100*(-1)^n)*sin(n*Pi*x)*exp(-Pi^2*n^2*t)/(Pi*n), n = 1 .. infinity)+30*x

(18)
> 

pde[10] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0; bc[10] := u(x, 0) = 0, u(x, 1) = 0, u(0, y) = y^2-y, u(1, y) = 0

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

u(x, 0) = 0, u(x, 1) = 0, u(0, y) = y^2-y, u(1, y) = 0

(19)
> 

pdsolve([pde[10], bc[10]], u(x, y))

u(x, y) = Sum(-4*((-1)^n-1)*sin(Pi*y*n)*(exp(Pi*n*(2*x-1))-exp(Pi*n))*exp(-Pi*n*(x-1))/((exp(2*Pi*n)-1)*Pi^3*n^3), n = 1 .. infinity)

(20)
> 

pde[11] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0

> 

bc[11] := (D[1](u))(0, y) = 0, (D[1](u))(L, y) = 0, u(x, H) = f(x), u(x, 0) = 0

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

(D[1](u))(0, y) = 0, (D[1](u))(L, y) = 0, u(x, H) = f(x), u(x, 0) = 0

(21)
> 

`assuming`([pdsolve([pde[11], bc[11]], u(x, y))], [0 < L, 0 < H])

u(x, y) = Sum(2*(Int(cos(Pi*x*n/L)*f(x), x = 0 .. L))*cos(Pi*x*n/L)*exp(Pi*n*(H-y)/L)*(exp(2*Pi*y*n/L)-1)/(L*(exp(2*Pi*H*n/L)-1)), n = 1 .. infinity)

(22)
> 

pde[12] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0

> 

bc[12] := (D[1](u))(L, y) = 0, u(x, H) = 0, u(x, 0) = 0, (D[1](u))(0, y) = g(y)

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

(D[1](u))(L, y) = 0, u(x, H) = 0, u(x, 0) = 0, (D[1](u))(0, y) = g(y)

(23)
> 

`assuming`([pdsolve([pde[12], bc[12]], u(x, y))], [0 < x, x <= L, 0 < y, y <= H])

u(x, y) = Sum(-2*(Int(sin(Pi*y*n/H)*g(y), y = 0 .. H))*sin(Pi*y*n/H)*(exp(-Pi*n*(L-2*x)/H)+exp(Pi*L*n/H))*exp(Pi*n*(L-x)/H)/(Pi*n*(exp(2*Pi*L*n/H)-1)), n = 1 .. infinity)

(24)
> 

pde[13] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0

> 

bc[13] := (D[1](u))(0, y) = 0, u(x, 0) = 0, u(x, H) = 0, u(L, y) = g(y)

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

(D[1](u))(0, y) = 0, u(x, 0) = 0, u(x, H) = 0, u(L, y) = g(y)

(25)
> 

`assuming`([pdsolve([pde[13], bc[13]], u(x, y))], [0 < L, 0 < H])

u(x, y) = Sum(2*(Int(sin(Pi*y*n/H)*g(y), y = 0 .. H))*sin(Pi*y*n/H)*exp(Pi*n*(L-x)/H)*(exp(2*Pi*x*n/H)+1)/(H*(exp(2*Pi*L*n/H)+1)), n = 1 .. infinity)

(26)
> 

pde[14] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0

> 

bc[14] := u(0, y) = g(y), u(L, y) = 0, (D[2](u))(x, 0) = 0, u(x, H) = 0

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

u(0, y) = g(y), u(L, y) = 0, (D[2](u))(x, 0) = 0, u(x, H) = 0

(27)
> 

`assuming`([pdsolve([pde[14], bc[14]], u(x, y))], [0 < x, x <= L, 0 < y, y <= H])

u(x, y) = Sum(-2*(exp(-(L-2*x)*(1/2+n)*Pi/H)-exp((1/2)*Pi*(1+2*n)*L/H))*cos((1/2)*Pi*(1+2*n)*y/H)*(Int(cos((1/2)*Pi*(1+2*n)*y/H)*g(y), y = 0 .. H))*exp((1/2)*Pi*(1+2*n)*(L-x)/H)/(H*(exp(Pi*(1+2*n)*L/H)-1)), n = 0 .. infinity)

(28)
> 

pde[15] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0

> 

bc[15] := u(0, y) = 0, u(L, y) = 0, u(x, 0) = (D[2](u))(x, 0), u(x, H) = f(x)

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

u(0, y) = 0, u(L, y) = 0, u(x, 0) = (D[2](u))(x, 0), u(x, H) = f(x)

(29)
> 

`assuming`([pdsolve([pde[15], bc[15]], u(x, y))], [0 < x, x <= L, 0 < y, y <= H])

u(x, y) = Sum(2*exp(Pi*n*(H-y)/L)*sin(Pi*x*n/L)*((Pi*n+L)*exp(2*Pi*y*n/L)+Pi*n-L)*(Int(sin(Pi*x*n/L)*f(x), x = 0 .. L))/(L*((Pi*n+L)*exp(2*Pi*H*n/L)+Pi*n-L)), n = 1 .. infinity)

(30)
> 

pde[16] := diff(u(x, t), t, t) = c^2*(diff(u(x, t), x, x))

> 

bc[16] := u(0, t) = 0, (D[1](u))(L, t) = 0, (D[2](u))(x, 0) = 0, u(x, 0) = f(x)

diff(diff(u(x, t), t), t) = c^2*(diff(diff(u(x, t), x), x))

 

u(0, t) = 0, (D[1](u))(L, t) = 0, (D[2](u))(x, 0) = 0, u(x, 0) = f(x)

(31)
> 

`assuming`([pdsolve([pde[16], bc[16]], u(x, t))], [0 < x, x <= L])

u(x, t) = Sum(2*cos((1/2)*c*Pi*(1+2*n)*t/L)*(Int(sin((1/2)*Pi*(1+2*n)*x/L)*f(x), x = 0 .. L))*sin((1/2)*Pi*(1+2*n)*x/L)/L, n = 0 .. infinity)

(32)
> 

pde[17] := diff(w(x1, x2, x3, t), t) = diff(w(x1, x2, x3, t), x1, x1)+diff(w(x1, x2, x3, t), x2, x2)+diff(w(x1, x2, x3, t), x3, x3)

> 

bc[17] := w(x1, x2, x3, 0) = exp(x1)+x2*x3^5

diff(w(x1, x2, x3, t), t) = diff(diff(w(x1, x2, x3, t), x1), x1)+diff(diff(w(x1, x2, x3, t), x2), x2)+diff(diff(w(x1, x2, x3, t), x3), x3)

 

w(x1, x2, x3, 0) = exp(x1)+x2*x3^5

(33)
> 

pdsolve([pde[17], bc[17]])

w(x1, x2, x3, t) = Sum(t^n*((proc (U) options operator, arrow; diff(diff(U, x1), x1)+diff(diff(U, x2), x2)+diff(diff(U, x3), x3) end proc)@@n)(exp(x1)+x2*x3^5)/factorial(n), n = 0 .. infinity)

(34)
> 

pde[18] := diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = -x

> 

bc[18] := u(x, 0) = x

diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = -x

 

u(x, 0) = x

(35)
> 

pdsolve([pde[18], bc[18]], u(x, t))

u(x, t) = x/tan(t+(1/4)*Pi)

(36)
> 

pde[19] := diff(u(x, t), t)-u(x, t)^2*(diff(u(x, t), x)) = 3*u(x, t)

> 

bc[19] := u(x, 0) = x

diff(u(x, t), t)-u(x, t)^2*(diff(u(x, t), x)) = 3*u(x, t)

 

u(x, 0) = x

(37)
> 

pdsolve([pde[19], bc[19]], u(x, t))

u(x, t) = -exp(3*t)*((-6*exp(6*t)*x+6*x+9)^(1/2)-3)/(exp(6*t)-1)

(38)
> 

pde[20] := diff(u(x, t), t)-x*u(x, t)*(diff(u(x, t), x)) = 0

> 

bc[20] := u(x, 0) = x

diff(u(x, t), t)-x*u(x, t)*(diff(u(x, t), x)) = 0

 

u(x, 0) = x

(39)
> 

pdsolve([pde[20], bc[20]], u(x, t))

u(x, t) = -LambertW(-t*x)/t

(40)
> 

pde[21] := diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = 0

> 

bc[21] := u(x, 0) = x

diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = 0

 

u(x, 0) = x

(41)
> 

pdsolve([pde[21], bc[21]], u(x, t))

u(x, t) = x/(t+1)

(42)
> 

pde[22] := (t+u(x, t))*(diff(u(x, t), x))+t*(diff(u(x, t), t)) = 0; bc[22] := u(x, 1) = x

(t+u(x, t))*(diff(u(x, t), x))+t*(diff(u(x, t), t)) = 0

 

u(x, 1) = x

(43)
> 

pdsolve([pde[22], bc[22]], u(x, t))

u(x, t) = (t-x-1)/(ln(1/t)-1)

(44)
> 

pde[23] := (diff(r*(diff(u(r, theta), r)), r))/r+(diff(u(r, theta), theta, theta))/r^2 = 0; bc[23] := u(a, theta) = f(theta), u(r, -Pi) = u(r, Pi), (D[2](u))(r, -Pi) = (D[2](u))(r, Pi)

(diff(u(r, theta), r)+r*(diff(diff(u(r, theta), r), r)))/r+(diff(diff(u(r, theta), theta), theta))/r^2 = 0

 

u(a, theta) = f(theta), u(r, -Pi) = u(r, Pi), (D[2](u))(r, -Pi) = (D[2](u))(r, Pi)

(45)
> 

`assuming`([pdsolve([pde[23], bc[23]], u(r, theta), HINT = boundedseries)], [a > 0])

u(r, theta) = (1/2)*(2*(Sum(r^n*(sin(n*theta)*(Int(f(theta)*sin(n*theta), theta = -Pi .. Pi))+cos(n*theta)*(Int(f(theta)*cos(n*theta), theta = -Pi .. Pi)))*a^(-n)/Pi, n = 1 .. infinity))*Pi+Int(f(theta), theta = -Pi .. Pi))/Pi

(46)
> 

pde[24] := diff(g(t, x), t) = diff(g(t, x), x, x)+a*g(t, x); bc[24] := g(0, x) = f(x)

diff(g(t, x), t) = diff(diff(g(t, x), x), x)+a*g(t, x)

 

g(0, x) = f(x)

(47)
> 

pdsolve([pde[24], bc[24]])

g(t, x) = exp(a*t)*invfourier(fourier(f(x), x, s1)*exp(-t*s1^2), s1, x)

(48)
> 

pde[25] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0; bc[25] := u(0, y) = y*(-y+1), u(1, y) = 0, (D[2](u))(x, 0) = 0, (D[2](u))(x, 1) = 0

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

u(0, y) = y*(-y+1), u(1, y) = 0, (D[2](u))(x, 0) = 0, (D[2](u))(x, 1) = 0

(49)
> 

pdsolve([pde[25], bc[25]], u(x, y))

u(x, y) = Sum(2*((-1)^n+1)*cos(Pi*y*n)*(exp(Pi*n*(2*x-1))-exp(Pi*n))*exp(-Pi*n*(x-1))/((exp(2*Pi*n)-1)*Pi^2*n^2), n = 1 .. infinity)

(50)
> 

``


 

Download ImprovementsInPdsolve.mw

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

And the Nobel prize in physics 2017 went for work in General Relativity ! Actually, experimental work involving sophisticated detectors and Numerical Relativity, one of the branches of GR. The prize was awarded to Rainer Weiss (85 years old, 1/2 of the prize), Barry Barish (81 years old, 1/4 of the prize) and Kip Thorne (77 years old, 1/4 of the prize) who have "shaken the world again" with their work on Ligo experiment, which was able to detect ripples in the fabric of spacetime.

General Relativity continues to be at the center of work in theoretical and experimental physics. I take this opportunity to note that, in Maple 2017, among the several improvements in the Physics package regarding General Relativity, there is a new package, Physics:-ThreePlusOne, all dedicated to the symbolic manipulations necessary to formulate problems in Numerical Relativity.

The GR functionality implemented in Physics, Physics:-Tetrads and Physics:-ThreePlusOne is unique in computer algebra systems and reflects the Maplesoft intention, for several years now, to provide the very best possible computer algebra environment for Physics, regarding current research activity as well as related education in advanced mathematical-physics methods.

For what is going on in theoretical physics nowadays and its connection with General Relativity, in very short: the unification of gravity with the other forces, check for instance this map from Aug/2015 (by the way a very nice summary for whoever is interested):

It is also interesting the article behind this map of topics as well as this brilliant and accessible presentation by Nima Arkani-Hamed (Princeton):  Quantum Mechanics and Spacetime in the 21st Century, given in the Perimeter Institute for Theoretical Physics (Waterloo), November 2014. 

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

I'm back from presenting work in the "23rd Conference on Applications of Computer Algebra -2017" . It was a very interesting event. This fifth presentation, about "The Appell doubly hypergeometric functions", describes a very recent project I've been working at Maple, i.e. the very first complete computational implementation of the Appell doubly hypergeometric functions. This work appeared in Maple 2017. These functions have a tremendous potential in that, at the same time, they have a myriad of properties, and include as particular cases most of the existing mathematical language, and so they have obvious applications in integration, differential equations, and applied mathematics all around. I think these will be the functions of this XXI century, analogously to what happened with hypergeometric functions in the previous century.

At the end, there is a link to the presentation worksheet, with which one could open the sections and reproduce the presentation examples.
 

The four double-hypergeometric Appell functions,

a complete implementation in a computer algebra system

 

Edgardo S. Cheb-Terrab

Physics, Differential Equations and Mathematical Functions, Maplesoft

 

Abstract:
The four multi-parameter Appell functions, AppellF1 , AppellF2 , AppellF3  and AppellF4  are doubly hypergeometric functions that include as particular cases the 2F1 hypergeometric  and some cases of the MeijerG  function, and with them most of the known functions of mathematical physics. Appell functions have been popping up with increasing frequency in applications in quantum mechanics, molecular physics, and general relativity. In this talk, a full implementation of these functions in the Maple computer algebra system, including, for the first time, their numerical evaluation over the whole complex plane, is presented, with details about the symbolic and numerical strategies used.

Appell Functions (symbolic)

 

 

The main references:

• 

P. Appel, J.Kamke de Feriet, "Fonctions hypergeometriques et Hyperspheriques", 1926

• 

H. Srivastava, P.W. Karlsson, "Multiple Gaussian Hypergeometric Series", 1985

• 

24 papers in the literature, ranging from 1882 to 2015

 

Definition and Symmetries

   

Polynomial and Singular Cases

   

Single Power Series with Hypergeometric Coefficients

   

Analytic Extension from the Appell Series to the Appell Functions

   

Euler-Type and Contiguity Identities

   

Appell Differential Equations

   

Putting all together

   

Problem: some formulas in the literature are wrong or miss the conditions indicating when are they valid (exchange with the Mathematics director of the DLMF - NIST)

   

Appell Functions (numeric)

 

 

Goals

 

• 

Compute these Appell functions over the whole complex plane

• 

Considering that this is a research problem, implement different methods and flexible optional arguments to allow for:

a) comparison between methods (both performance and correctness),

b) investigation of a single method in different circumstances.

• 

Develop a computational structure that can be reused with other special functions (abstract code and provide the main options), and that could also be translated to C (so: only one numerical implementation, not 100 special function numerical implementations)

Limitation: the Maple original evalf command does not accept optional arguments

 

The cost of numerically evaluating an Appell function

 

• 

If it is a special hypergeometric case, then between 1 to 2 hypergeometric functions

• 

Next simplest case (series/recurrence below) 3 to 4 hypergeometric functions plus adding somewhat large formulas that involve only arithmetic operations up to 20,000 times (frequently less than 100 times)

• 

Next simplest case: the formulas themselves are power series with hypergeometric function coefficients; these cases frequently converge rapidly but may involve the numerical evaluation of up to hundreds of hypergeometric functions to get the value of a single Appell function.

 

Strategy for the numerical evaluation of Appell functions (or other functions ...)

 

 

The numerical evaluation flows orderly according to:

1) check whether it is a singular case

2) check whether it is a special value

3) compute the value using a series derived from a recurrence related to the underlying ODE

4) perform an sum using an infinite sum formula, checking for convergence

5) perform the numerical integration of the ODE underlying the given Appell function

6) perform a sequence of concatenated Taylor series expansions

Examples

   

Series/recurrence

   

Numerical integration of an underlying differential equation (ODEs and dsolve/numeric)

   

Concatenated Taylor series expansions covering the whole complex plane

   

Subproducts

 

Improvements in the numerical evaluation of hypergeometric functions

   

Evalf: an organized structure to implement the numerical evaluation of special functions in general

   

To be done

   


 

Download Appell_Functions.mw   
Download Appell_Functions.pdf

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

First 8 9 10 11 12 13 14 Last Page 10 of 18