Education

What is new in Physics in Maple 2025

by: ecterrab 14320 Product: Maple

17 hours ago

6 0

What is new in Physics in Maple 2025

This post is, basically, the page of what is new in Physics distributed with Maple 2025. Why post it? Because this time we achieved a result that is a breakthrough/milestone in Computer Algebra: for the first time we can systematically compute Einstein's equations from first principles, using a computer, starting from a Lagrangian for gravity. This is something to celebrate.

Although being able to perform this computation on a computer algebra sheet is relevant mostly for physicists working in general relativity, the developments performed in the tensor manipulation and functional differentiation routines of the Maple system to cover this computation were tremendous, in number of lines of code, efforts and time, resulting in improvements not just fore general relativity but for physics all around, and in an area out of reach of neural network AIs . Other results, like the new routines for linearizing gravity, requested other times here in Mapleprimes, are also part of the relevant novelties.

Lagrange Equations and simplification of tensorial expressions in curved spacetimes

Linearized Gravity

Relative Tensors

New Physics:-Library commands

See Also

Lagrange Equations and simplification of tensorial expressions in curved spacetimes

LagrangeEquations is a Physics command introduced in 2023 taking advantage of the functional differentiation capabilities of the Physics package . This command can handle tensors and vectors of the Physics package as well as derivatives using vectorial differential operators (see d_ and Nabla ), works by performing functional differentiation (see Fundiff ), and handles 1^st, and higher order derivatives of the coordinates in the Lagrangian automatically. LagrangeEquations receives an expression representing a Lagrangian and returns a sequence of Lagrange equations with as many equations as coordinates are indicated. The number of parameters can also be many. For example, in electrodynamics, the "coordinate" is a tensor field , there are then four coordinates, one for each of the values of the index , and there are four parameters .

New in Maple 2025, the "coordinates" can now also be the components of the metric tensor in a curved spacetime, in which case the equations returned are Einstein's equations. Also new, instead of a coordinate or set of them, you can pass the keyword EnergyMomentum, in which case the output is the conserved energy-momentum tensor of the physical model represented by the given Lagrangian L.

Examples

(1)

The model in classical field theory and corresponding field equations, as in previous releases

>

(2)

>

(1/2)*Physics:-d_[mu](Phi(X), [X])*Physics:-d_[`~mu`](Phi(X), [X])-(1/2)*m^2*Phi(X)^2+(1/4)*lambda*Phi(X)^4

(3)

Lagrange's equations

>

(4)

New: The energy-momentum tensor can be computed as the Lagrange equations taking the metric as the coordinate, not equating to 0 the result, but multiplying the variation of the action "(delta S)/(delta g[]^(mu,nu))" by (in flat spacetimes ). For that purpose, you can use the EnergyMomentum keyword. You can optionally indicate the indices to be used in the output as well as their covariant or contravariant character

>

Physics:-EnergyMomentum[mu, nu] = (-(1/4)*lambda*Phi(X)^4+(1/2)*m^2*Phi(X)^2-(1/2)*Physics:-d_[`~beta`](Phi(X), [X])*Physics:-d_[beta](Phi(X), [X]))*Physics:-g_[mu, nu]+Physics:-d_[mu](Phi(X), [X])*Physics:-d_[nu](Phi(X), [X])

(5)

To further compute using the above as the definition for , you can use the Define command

>

${Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-EnergyMomentum[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}$

(6)

After which the system knows about the symmetry properties and the components of

>

Physics:-EnergyMomentum[mu, nu] = (-(1/4)*lambda*Phi(X)^4+(1/2)*m^2*Phi(X)^2-(1/2)*Physics:-d_[`~beta`](Phi(X), [X])*Physics:-d_[beta](Phi(X), [X]))*Physics:-g_[mu, nu]+Physics:-d_[mu](Phi(X), [X])*Physics:-d_[nu](Phi(X), [X])

(7)

>

(8)

>

Physics:-EnergyMomentum[mu, nu] = Matrix(%id = 36893488151952536988)

(9)

New: LagrangeEquations takes advantage of the extension of Fundiff to compute functional derivatives in curved spacetimes introduced for Maple 2025, and so it also handles the case of a scalar field in a curved spacetime. Set for instance an arbitrary metric

>

Physics:-g_[mu, nu] = Matrix(%id = 36893488152257720060)

(10)

For the action to be a true scalar in spacetime, the Lagrangian density now needs to be multiplied by the square root of the determinant of the metric

>

(-%g_[determinant])^(1/2)*((1/2)*Physics:-d_[mu](Phi(X), [X])*Physics:-d_[`~mu`](Phi(X), [X])-(1/2)*m^2*Phi(X)^2+(1/4)*lambda*Phi(X)^4)

(11)

New: With the extension of the tensorial simplification algorithms for curved spacetimes, the Lagrange equations can be computed arriving directly to the compact form

>

(12)

Comparing with the result (4) for the same Lagrangian in a flat spacetime, we see the only difference is that the dAlembertian is now expressed in terms of covariant derivatives D_ .

The EnergyMomentum tensor is computed in the same way as when the spacetime is flat

>

Physics:-EnergyMomentum[mu, nu] = (-(1/4)*lambda*Phi(X)^4+(1/2)*m^2*Phi(X)^2-(1/2)*Physics:-d_[`~beta`](Phi(X), [X])*Physics:-d_[beta](Phi(X), [X]))*Physics:-g_[mu, nu]+Physics:-d_[mu](Phi(X), [X])*Physics:-d_[nu](Phi(X), [X])

(13)

General Relativity

New: the most significant development in LagrangeEquations is regarding General Relativity. It can now compute Einstein's equations directly from the Lagrangian, not using tabulated cases, and properly handling several (traditional or not) alternative ways of presenting the Lagrangian.

Einstein's equations concern the case of a curved spacetime with metric as, for instance, the general case of an arbitrary metric set lines above. In the Lagrangian formulation, the coordinates of the problem are the components of the metric , and as in the case of electrodynamics the parameters are the spacetime coordinates . The simplest case is that of Einstein's equation in vacuum, for which the Lagrangian density is expressed in terms of the trace of the Ricci tensor by

>

(14)

Einstein's equations in vacuum:

>

-(1/2)*Physics:-g_[mu, nu]*Physics:-Ricci[alpha, `~alpha`]+Physics:-Ricci[mu, nu] = 0

(15)

where in the above instead of passing as second argument, we passed to get the equations using those free indices. The tensorial equation computed is also the definition of the Einstein tensor

>

Physics:-Einstein[mu, nu] = -(1/2)*Physics:-g_[mu, nu]*Physics:-Ricci[alpha, `~alpha`]+Physics:-Ricci[mu, nu]

(16)

The Lagrangian used to compute Einstein's equations (15) contains first and second derivatives of the metric. To see that, rewrite L in terms of Christoffel symbols

>

(17)

Recalling the definition

>

Physics:-Christoffel[alpha, mu, nu] = (1/2)*Physics:-d_[nu](Physics:-g_[alpha, mu], [X])+(1/2)*Physics:-d_[mu](Physics:-g_[alpha, nu], [X])-(1/2)*Physics:-d_[alpha](Physics:-g_[mu, nu], [X])

(18)

in the two terms containing derivatives of Christoffel symbols contain second order derivatives of . Now, it is always possible to add a total spacetime derivative to without changing Einstein's equations (assuming the variation of the metric in the corresponding boundary integrals vanishes), and in that way, in this particular case of , obtain a Lagrangian involving only 1st order derivatives. The total derivative, expressed using the inert command to see it before the differentiation operation is performed, is

>

(19)

Adding this term to , performing the differentiation operation and simplifying we get

>

(20)

>

(21)

>

(22)

which is a Lagrangian depending only on 1st order derivatives of the metric through Christoffel symbols. As expected, the equations of motion resulting from this Lagrangian are the same Einstein equations computed in (15)

>

-(1/2)*Physics:-Ricci[iota, `~iota`]*Physics:-g_[mu, nu]+Physics:-Ricci[mu, nu] = 0

(23)

To illustrate the new Maple 2025 tensorial simplification capabilities note that is no just ≡ (17) after discarding its two terms involving derivatives of Christoffel symbols. To verify this, split into the terms containing or not derivatives of Christoffel

>

(24)

Comparing, the total derivative TD≡ (19) is not just , but

>

(25)

Things like these, , can now be verified directly with the new tensorial simplification capabilities: take the left-hand side minus the right-hand side, evaluate the inert derivative and simplify to see the equality is true

>

(26)

>

(27)

>

(28)

That said, it is also true that results in the Lagrangian , and since the equations of movement don't depend on the sign of the Lagrangian, for this Lagrangian adding the term happens to be equivalent to just discarding the terms of involving derivatives of Christoffel symbols.

Also new in Maple 2025, due to the extension of Fundiff to compute in curved spacetimes, it is now also possible to compute Einstein's equations from first principles by constructing the action,

>

Int(Int(Int(Int((-%g_[determinant])^(1/2)*Physics:-Ricci[alpha, `~alpha`], x = -infinity .. infinity), y = -infinity .. infinity), z = -infinity .. infinity), t = -infinity .. infinity)

(29)

and equating to zero the functional derivative with respect to the metric. To avoid displaying the resulting large expression, end the input line with ":"

>

Simplifying this result, we get an expression in terms of Christoffel symbols and its derivatives

>

(30)

In this result, we see derivatives of Christoffel symbols, expressed using the D_ command for covariant differentiation. Although, such objects have not the geometrical meaning of a covariant derivative, computationally, they here represent what would be a covariant derivative if the Christoffel symbols were a tensor. For example,

>

(31)

With this computational meaning for the derivatives of Christoffel symbols appearing in (30), rewrite EEC≡(30) in terms of the Ricci and Riemann tensors. For that, consider the definition

>

(32)

Rewrite the noncovariant derivatives in terms of derivatives using the computational representation (31), simplify and isolate one of them

>

(33)

>

(34)

>

(35)

Analogously, derive an expression to rewrite derivatives of Christoffel symbols using the Riemann tensor

>

(36)

>

(37)

>

(38)

>

(39)

Substitute these two equations, in sequence, into Einstein's equations EEC≡(30)

>

(40)

Simplify to arrive at the traditional compact form of Einstein's equations

>

(1/2)*Physics:-Ricci[chi, `~chi`]*Physics:-g_[`~alpha`, `~beta`]-Physics:-Ricci[`~alpha`, `~beta`] = 0

(41)

Linearized Gravity

Generally speaking, linearizing gravity is about discarding in Einstein's field equations the terms that are quadratic in the metric and its derivatives, an approximation valid when the gravitational field is weak (the deviation from a flat Minkowski spacetime is small). Linearizing gravity is used, e.g. in the study of gravitational waves. In the context of Maple's Physics , the formulation of linearized gravity can be done using the general relativity tensors that come predefined in Physics plus a new in Maple 2025 Physics:-Library:-Linearize command.

In what follows it is shown how to linearize the Ricci tensor and through it Einstein 's equations. To compare results, see for instance the Wikipedia page for Linearized gravity. Start setting coordinates, you could use Cartesian, spherical, cylindrical, or define your own.

>

restart; with(Physics); Setup(coordinates = cartesian)

(42)

The default metric when Physics is loaded is the Minkowski metric, representing a flat (no curvature) spacetime

>

Physics:-g_[mu, nu] = Matrix(%id = 36893488151987392132)

(43)

The weakly perturbed metric

Suppose you want to define a small perturbation around this metric. For that purpose, define a perturbation tensor , that in the general case depends on the coordinates and is not diagonal, the only requirement is that it is symmetric (to have it diagonal, change symmetric by diagonal; to have it constant, change by )

>

h[mu, nu] = Matrix(%id = 36893488152568729716)

(44)

In the above it is understood that , so that quadratic or higher powers of it or its derivatives can be approximated to 0 (discarded). Define the components of accordingly

>

${Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-d_[mu], Physics:-g_[mu, nu], h[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}$

(45)

Define also a tensor representing the unperturbed Minkowski metric

>

eta[mu, nu] = Matrix(%id = 36893488151987392132)

(46)

>

${Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-d_[mu], eta[mu, nu], Physics:-g_[mu, nu], h[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}$

(47)

The weakly perturbed metric is given by

>

(48)

Make this be the definition of the metric

>

Physics:-g_[mu, nu] = Matrix(%id = 36893488152066920564)

${Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], eta[mu, nu], Physics:-g_[mu, nu], Physics:-gamma_[i, j], h[mu, nu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}$

(49)

Linearizing the Ricci tensor

The linearized form of the Ricci tensor is computed by introducing this weakly perturbed metric (48) in the expression of the Ricci tensor as a function of the metric. This can be accomplished in different ways, the simpler being to use the conversion network between tensors , but for illustration purposes, showing steps one at time, a substitution of definitions one into the other one is used

>

(50)

>

Physics:-Christoffel[`~alpha`, mu, nu] = (1/2)*Physics:-g_[`~alpha`, `~beta`]*(Physics:-d_[nu](Physics:-g_[beta, mu], [X])+Physics:-d_[mu](Physics:-g_[beta, nu], [X])-Physics:-d_[beta](Physics:-g_[mu, nu], [X]))

(51)

>

(52)

Introducing (48)≡, and also the inert form of the Ricci tensor to facilitate simplification some steps below,

>

(53)

This expression contains several terms quadratic in the small perturbation and its derivatives. The new in Maple 2025 routine to filter out those terms is Physics:-Library:-Linearize , which requires specifying the symbol representing the small quantities

>

(54)

Important: in this result, is the flat Minkowski metric, not the perturbed metric . However, in the context of a linearized formulation, raises and lowers tensor indices the same way as . Hence, to further simplify contracted products of in (54) , it is practical to reintroduce representing that Minkowski metric and simplify using the internal algorithms for a flat metric

>

Physics:-g_[mu, nu] = Matrix(%id = 36893488152066895868)

(55)

To proceed simplifying, replace in the expression (54) for the Ricci tensor the intermediate Minkowski by

>

(56)

Simplifying, results in the linearized form of the Ricci tensor shown in the Wikipedia page for Linearized gravity.

>

%Ricci[mu, nu] = -(1/2)*Physics:-d_[mu](Physics:-d_[nu](h[tau, `~tau`], [X]), [X])-(1/2)*Physics:-dAlembertian(h[mu, nu], [X])+(1/2)*Physics:-d_[nu](Physics:-d_[tau](h[mu, `~tau`], [X]), [X])+(1/2)*Physics:-d_[mu](Physics:-d_[tau](h[nu, `~tau`], [X]), [X])

(57)

Linearizing Einstein's equations

Einstein's equations are the components of Einstein's tensor , whose definition in terms of the Ricci tensor is

>

Physics:-Einstein[mu, nu] = Physics:-Ricci[mu, nu]-(1/2)*Physics:-g_[mu, nu]*Physics:-Ricci[alpha, `~alpha`]

(58)

Compute the trace directly from the linearized form (57) of the Ricci tensor,

>

%Ricci[mu, nu]*Physics:-g_[`~mu`, `~nu`] = (-(1/2)*Physics:-d_[mu](Physics:-d_[nu](h[tau, `~tau`], [X]), [X])-(1/2)*Physics:-dAlembertian(h[mu, nu], [X])+(1/2)*Physics:-d_[nu](Physics:-d_[tau](h[mu, `~tau`], [X]), [X])+(1/2)*Physics:-d_[mu](Physics:-d_[tau](h[nu, `~tau`], [X]), [X]))*Physics:-g_[`~mu`, `~nu`]

(59)

>

(60)

The linearized Einstein equations are constructed reproducing the definition (58) using (57) and (60)

>

(61)

which is the same formula shown in the Wikipedia page for Linearized gravity.

You can now redefine the general introduced in (44) in different ways (see discussion in the Wikipedia page), or, depending on the case, just substitute your preferred gauge in this formula (61) for the general case. For example, the condition for the Harmonic gauge also known as Lorentz gauge reduces the linearized field equations to their simplest form

>

Physics:-d_[mu](h[`~mu`, nu], [X]) = (1/2)*Physics:-d_[nu](h[alpha, `~alpha`], [X])

(62)

>

(63)

>

%Ricci[mu, nu]-(1/2)*Physics:-g_[mu, nu]*%Ricci[alpha, `~alpha`] = -(1/2)*Physics:-dAlembertian(h[mu, nu], [X])+(1/4)*Physics:-dAlembertian(h[alpha, `~alpha`], [X])*Physics:-g_[mu, nu]

(64)

Relative Tensors

In General Relativity, the context of a curved spacetime, it is sometimes necessary to work with relative tensors, for which the transformation rule under a transformation of coordinates involves powers of the determinant of the transformation - see Chapter 4 of "Lovelock, D., and Rund, H. Tensors, Differential Forms and Variational Principles, Dover, 1989." Physics in Maple 2025 includes a complete, new implementation of relative tensors.

To indicate that a tensor being defined is relative pass its relative weight. For example, set a curved spacetime,

>

restart; with(Physics); g_[sc]

$`Default differentiation variables for d_, D_ and dAlembertian are:`*{X = (r, theta, phi, t)}$

Physics:-g_[mu, nu] = Matrix(%id = 36893488152573138332)

(65)

Define now two tensors of one index, one of them being relative

>

${Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], T[mu], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[i, j], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}$

(66)

>

${Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], R[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], T[mu], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[i, j], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}$

(67)

Transformation of Coordinates

Consider a transformation of coordinates, from spherical to where

>

TR := r = (1+m/(2*rho))^2*rho

r = (1+(1/2)*m/rho)^2*rho

(68)

The transformed components of and are, respectively,

>

Vector[column](%id = 36893488151820263660)

(69)

>

Vector[column](%id = 36893488151823155676)

(70)

where, when comparing both results, we see that the transformed components for are all multiplied by with and is the determinant of the transformation:

>

Matrix(%id = 36893488151964220700)

(71)

>

J = (1/4)*(-m^2+4*rho^2)/rho^2

(72)

Relative weight

The relative weight of a scalar, tensor or tensorial expression can be computed using the Physics:-Library:-GetRelativeWeight command. For the two tensors and used above,

>

(73)

>

(74)

The relative weight of a tensor does not depend on the covariant or contravariant character of its indices

>

(75)

The LeviCivita tensor is a special case, has its relative weight defined when Physics is loaded, and because in a curved spacetime it is not a tensor its relative weight depends on the covariant or contravariant character of its indices

>

(76)

>

(77)

The relative weight w of a product is equal to the sum of relative weights of each factor

>

(78)

>

(79)

The relative weight w of a power is equal to the relative weight of the base multiplied by the power

>

1/(R[mu]*R[`~mu`])

(80)

(81)

The relative weight w of a sum is equal to the relative weight of one of its terms and exists if all the terms have the same w.

>

(82)

>

(83)

The relative weight of any determinant is always equal to 2

>

(84)

>

(85)

Relative Term in covariant derivatives

When computing the covariant derivative of a relative scalar, tensor or tensorial expression that has non-zero relative weight w, a relative term is added, that can be computed using the Physics:-Library:-GetRelativeWeight command.

>

(86)

>

(87)

Consequently,

>

(88)

To understand this zero value on the right-hand side, express the left-hand side in terms of d_

>

(89)

evaluate the inert %d_

>

(90)

The factor in parentheses is equal to , where the covariant derivative of the metric is equal to zero, so

>

(91)

Consider the covariant derivative of and defined in (66) and (67)

>

(92)

>

(93)

The corresponding covariant derivatives

>

(94)

>

%D_[mu](T[`~mu`](X)) = 2*T[`~mu`](X)*Physics:-d_[mu](r, [X])/r+Physics:-d_[mu](theta, [X])*cos(theta)*T[`~mu`](X)/sin(theta)+Physics:-d_[mu](T[`~mu`](X), [X])

(95)

>

(96)

>

%D_[mu](R[`~mu`](X)) = 2*R[`~mu`](X)*Physics:-d_[mu](r, [X])/r+Physics:-d_[mu](theta, [X])*cos(theta)*R[`~mu`](X)/sin(theta)+Physics:-d_[mu](R[`~mu`](X), [X])-Physics:-Christoffel[`~nu`, mu, nu]*R[`~mu`](X)

(97)

where in the above we see the additional (relative) term

>

(98)

New Physics:-Library commands

ConvertToF , Linearize , GetRelativeTerm , GetRelativeWeight.

Examples

•

ConvertToF receives an algebraic expression involving tensors and/or tensor functions and rewrites them in terms of the tensor of name F when that is possible. This routine is similar, however more general than the standard convert which only handles the existing conversion network for the tensors of General Relativity in that ConvertToF also uses any tensor definition you introduce using Define , expressing a tensor in terms of others.

Load any curved spacetime metric automatically setting the coordinates

restart; with(Physics); g_[sc]

Physics:-g_[mu, nu] = Matrix(%id = 36893488152573203868)

(99)

For example, rewrite the Christoffel symbols in terms of the metric g_ ; this works as in previous releases

>

Physics:-Christoffel[mu, alpha, beta] = (1/2)*Physics:-d_[beta](Physics:-g_[alpha, mu], [X])+(1/2)*Physics:-d_[alpha](Physics:-g_[beta, mu], [X])-(1/2)*Physics:-d_[mu](Physics:-g_[alpha, beta], [X])

(100)

Define a representing the 4D electromagnetic potential as a function of the coordinates and representing the electromagnetic field tensors

>

${A[mu], Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[i, j], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}$

(101)

>

${A[mu], Physics:-D_[mu], Physics:-Dgamma[mu], F[mu, nu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-gamma_[i, j], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}$

(102)

Rewrite the following expression in terms of the electromagnetic potential

>

(103)

In the example above, the output is similar to this other one

>

(104)

The rewriting, however, works also with tensorial expressions

>

(105)

>

(106)

•	Linearize receives a tensorial expression T and an indication of the small quantities h in T , and discards terms quadratic or of higher order in h. For an example of this new routine in action, see the section Linearized Gravity above.

•	GetRelativeTerm and GetRelativeWeight are illustrated in the section Relative Tensors above.

	See Also
	Index of New Maple 2025 Features , Physics , Computer Algebra for Theoretical Physics, The Physics project, The Physics Updates

Download New_in_Physics_2025.mw

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

An Introduction to Limits with Maple

by: Paras31 205 Product: Maple 2024

March 21 2025

2 0

LIMITS

Limits in maths are defined as the values that a function approaches the output for the given input values. Limits play a vital role in calculus and mathematical analysis and are used to define integrals, derivatives, and continuity. It is used in the analysis process, and it always concerns the behavior of the function at a particular point. The limit of a sequence is further generalized in the concept of the limit of a topological net and related to the limit and direct limit in the theory category. Generally, the integrals are classified into two types namely, definite and indefinite integrals. For definite integrals, the upper limit and lower limits are defined properly. Whereas indefinite integrals are expressed without limits, and it will have an arbitrary constant while integrating the function.

Sometimes we can't work something out directly ... but we can see what it should be as we get closer and closer!

Example 1

>

"restart; f(x):=(|x|-3)/(x-3);"

(1)

>

Error, (in f) numeric exception: division by zero

Now 0/0 is a difficulty! We don't really know the value of 0/0 (it is "indeterminate"), so we need another way of answering this.

So instead of trying to work it out for x=3 let's try approaching it closer and closer:

>

(2)

>

(3)

>

(4)

>

(5)

>

(6)

>

(7)

>

(8)

Example 2

Sometimes some functions are not continuous. That is, they appear to be approaching two different values when they are approached from two sides.

>

proc (x) options operator, arrow, function_assign; piecewise(0 < x and x < 2, 1/(2*x-x^2), 2 < x and x <= 3, 2-x, 3 < x and x < 4, x-4, 4 <= x, Pi, undefined) end proc

(9)

>

Suppose we want to approach 2 and see the function’s limit. This naturally leads to directions from which we can approach. Left-hand side and the right-hand side limits.

The right-hand side limit is the value of the function that it takes while approaching it from the right-hand side of the desired point. Similarly, the left-hand side limit is the value of function while approaching it from the left-hand side.

>

(10)

>

(11)

>

(12)

>

(13)

And the ordinary limit "does not exist".

>

(14)

>

(15)

>

(16)

>

(17)

And the ordinary limit "does not exist".

>

Example 3

Estimate the value of the following limit limit(h(x)*where, x = 2), h(x) = piecewise(x <> 2, x+12, x = 2, 4) .

>

"h(x):={[[x+12,x<>2],[4,x=2]];"

proc (x) options operator, arrow, function_assign; piecewise(x <> 2, x+12, x = 2, 4) end proc

(18)

>

(19)

The limit is NOT 2025!Remember from the first example that limits do not care what the function is actually doing at the point in question. Limits are only concerned with what is going on around the point. Since the only thing about the function that we actually changed was its behavior at this will not change the limit.

Example 4

>

proc (x) options operator, arrow, function_assign; piecewise(x < 0, -x+5, 0 <= x, 2*x) end proc

(20)

>

(21)

>

(22)

>

(23)

Example 5

>

proc (x) options operator, arrow, function_assign; piecewise(x < 5, x+4, 5 <= x, x^2-2) end proc

(24)

>

(25)

>

(26)

>

(27)

>

(28)

>

(29)

Example 6

>

proc (x) options operator, arrow, function_assign; piecewise(x <= 1, (x-8)/(x-3), 3 <= x, sqrt(x^2+x+2), undefined) end proc

(30)

>

(31)

>

(32)

>

(33)

>

(34)

Example 7

Estimate the value of the following limit. where, H(t) = piecewise(t < 0, 0, t >= 0, 1)

>

proc (t) options operator, arrow, function_assign; piecewise(t < 0, 0, 0 <= t, 1) end proc

(35)

This function is often called either the Heaviside or step function. We could use a table of values to estimate the limit, but it’s probably just as quick in this case to use the graph so let’s do that. Below is the graph of this function.

>

(36)

>

(37)

We can see from the graph that if we approach from the right side the function is moving in towards a value of 1. Well actually it’s just staying at 1, but in the terminology that we’ve been using in this section it’s moving in towards 1.

Also, if we move in towards from the left the function is moving in towards a value of 0.

According to our definition of the limit the function needs to move in towards a single value as we move in towards (from both sides). This isn’t happening in this case and so in this example we will also say that the limit doesn’t exist.

Download limits.mw

ΕΜΒΑΔΟΝ ΕΠΙΠΕΔΟΥ ΧΩΡΙΟΥ – Maple Worksheet for...

by: Paras31 205 Product: Maple 2024

March 10 2025

1 2

Hello everyone,

I have created a Maple worksheet titled "ΕΜΒΑΔΟΝ ΕΠΙΠΕΔΟΥ ΧΩΡΙΟΥ", designed to help my students prepare for their final exams as they qualify for university. This worksheet focuses on area calculations in plane geometry, using Maple to visualize and solve problems efficiently.

This worksheet is aimed at high school students preparing for university entrance exams, as well as teachers who want to integrate Maple into their teaching.

I would love to hear your thoughts and feedback!

Have you used Maple for similar exam preparation?
εμβαδόν_χωρίου.mw

Matt's simple prime checker for a^2+2

by: Matt C Anderson 185 Product: Maple 13

March 02 2025

1 0

Hi Maple lovers, and others,

I wrote up some simple Maple code that
checks if an expression (a^3+2) is a prime number.

The output fits on one page of paper.
n_cubed_plus_2.mw

n_cubed_plus_2.pdf

Also, I have an account at mersenneforum.org

Kindest Regards,

Matt

PS Does this help? The code works.

Spring 2025 issue of Maple Transactionss

by: rcorless 710 Product: Maple

February 25 2025

7 0

I'm happy to announce the publication of Volume 5, Issue #1 of Maple Transactions. You can find it at

mapletransactions.org

We have a survey paper by Veselin Jungic and Naomi Borwein on teaching Experimental Mathematics courses as our Featured Contribution. Many of you will find it interesting and useful.

In the refereed paper section we have a paper on Metaprogramming with Maple and C by Ilias Kotsireas; a paper on fast transposed Vandermonde solving by Hyukho Kwon & Michael Monagan; a paper by David Ulgenes (an honours student in Oslo) on Gamma, Pseudogamma, and Inverse Gamma functions; and a paper by John Campbell on applications of Gosper's nonlocal derangement identity (which, if you don't know that the word "derangement" has a technical mathematical meaning, may give you the wrong impression!).

As usual I've also written something, and I hope you like it: it's about Chladni figures and standing modes in an elliptical drum, and visualizing such in Maple. It uses Mathieu functions in Maple and noodles a bit about zerofinding (but winds up using fsolve because that's so convenient).

Keep the papers coming. This is the 12th issue of Maple Transactions, and I remind you that it has a "Diamond Class" designation, which means there are no page charges to authors, and the articles are free to read for everyone. This means that there's some volunteer labour needed, of course: you have to write the articles, and what we want is that you write articles that people in the Maple community actually want to read.

I'd also like to thank the copyeditor, Michelle Hatzel, for her very hard work on this issue. She's really made a difference, and I think you will be able to see it.

Another pancake cutting theorem

by: Kitonum 21308 Product: Maple 2018

February 24 2025

3 0

In the book by W.G. Chinn, N.E. Steenrod "First Concepts of Topology" the another remarkable theorem was proved: any two flat bounded regions can be cut by a single straight line so that each of these regions is divided into two regions of equal area (the second pancake problem). This is an existence theorem which does not provide any way to find this cut. In this post I made an attempt to find such cut for any 2 convex regions on the plane bounded by a piecewise smooth self-non-intersecting curves.
The Each_Into_2_Equal_Areas procedure returns a list of coordinates of 4 endpoints of the cutting segments. This procedure significantly uses my old procedures Area and Picture , which can be found in detail at the link https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure- . The formal arguments of the Each_Into_2_Equal_Areas procedure are the lists L1 and L2 specifying the boundaries of the regions to be cut. When specifying L1 and L2 , the boundary can be passed clockwise or counterclockwise, but it is necessary that the parameter t (when specifying each link) should go in ascending order. This can always be achieved by replacing t with -t if necessary. The Pic procedure draws a picture of the source regions and cutting segments. For ease of use, the code for the Area and Picture procedure is also provided. It is also worth noting that the procedure also works for "not too" non-convex regions (see examples below).

restart;
Area := proc(L) 
local i, var, e, e1, e2, P; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then 
P[i] := (1/2)*add(L[i, j, 1]*L[i, j+1, 2]-L[i, j, 2]*L[i, j+1, 1], j = 1 .. nops(L[i])-1) else 
var := lhs(L[i, 2]); 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := (1/2)*(int(e^2, L[i, 2])) else 
if var = y then P[i] := (1/2)*simplify(int(e-var*(diff(e, var)), L[i, 2])) else 
P[i] := (1/2)*simplify(int(var*(diff(e, var))-e, L[i, 2])) end if end if else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; 
P[i] := (1/2)*simplify(int(e1*(diff(e2, var))-e2*(diff(e1, var)), L[i, 2])) end if end if end do; 
abs(add(P[i], i = 1 .. nops(L))); 
end proc:

Picture := proc(L, C, N::posint := 100, Boundary::list := [linestyle = 1]) 
local i, var, var1, var2, e, e1, e2, P, Q, h; 
global Border;
uses plottools; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then P[i] := op(L[i]) else 
var := lhs(L[i, 2]); var1 := lhs(rhs(L[i, 2])); var2 := rhs(rhs(L[i, 2])); h := (var2-var1)/N; 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := seq(subs(var = var1+h*i, [e*cos(var), e*sin(var)]), i = 0 .. N) else 
P[i] := seq([var1+h*i, subs(var = var1+h*i, e)], i = 0 .. N) fi else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; P[i] := seq(subs(var = var1+h*i, [e1, e2]), i = 0 .. N) fi; fi; od; 
Q := [seq(P[i], i = 1 .. nops(L))]; Border := curve([op(Q), Q[1]], op(Boundary)); [polygon(Q, C), Border] 
end proc:

Each_Into_2_Equal_Areas:=proc(L1::list, L2::list)
local D, n, m, L10, L20, S1,S2, f, L11, L21, i, j, k, s, A, B, C , sol;

f:=(X,Y)->expand((y-X[2])*(Y[1]-X[1])-(x-X[1])*(Y[2]-X[2]));
L10:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L1);
L20:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L2);
S1:=Area(L1); S2:=Area(L2);  
n:=nops(L1); m:=nops(L2);

for i from 1 to n do
for j from i to n do

for k from 1 to m do
for s from k to m do

if not ((nops({i,j})=1 and type(L1[i],listlist)) or (nops({k,s})=1 and type(L2[k],listlist))) then

A:=eval(L10[i,1],t=t1): 
B:=eval(L10[j,1],t=t2):
C:=eval(L20[k,1],t=t3): 
D:=eval(L20[s,1],t=t4):

L11:=`if`(j=i,[subsop([2,2]=t1..t2,L10[i]),[B,A]],`if`(j=i+1,[subsop([2,2]=t1..op([2,2,2],L10[i]),L10[i]),subsop([2,2]=op([2,2,1],L10[j])..t2,L10[j]),[B,A]], [subsop([2,2]=t1..op([2,2,2],L10[i]),L10[i]),op(L10[i+1..j-1]),subsop([2,2]=op([2,2,1],L10[j])..t2,L10[j]),[B,A]])):

L21:=`if`(s=k,[subsop([2,2]=t3..t4,L20[k]),[D,C]],`if`(s=k+1,[subsop([2,2]=t3..op([2,2,2],L20[k]),L20[k]),subsop([2,2]=op([2,2,1],L20[s])..t4,L20[s]),[D,C]], [subsop([2,2]=t3..op([2,2,2],L20[k]),L20[k]),op(L20[k+1..s-1]),subsop([2,2]=op([2,2,1],L20[s])..t4,L20[s]),[D,C]])):

sol:=fsolve(simplify({Area(L11)-S1/2,Area(L21)-S2/2,eval(f(A,B),[x=C[1],y=C[2]]),eval(f(A,B),[x=D[1],y=D[2]])}),{t1=op([2,2,1],L10[i])..op([2,2,2],L10[i]),t2=op([2,2,1],L10[j])..op([2,2,2],L10[j]),t3=op([2,2,1],L20[k])..op([2,2,2],L20[k]),t4=op([2,2,1],L20[s])..op([2,2,2],L20[s])});

if type(sol,set(`=`)) then  return eval([A,B,C,D],sol) fi;

fi;
od: od: od: od:
end proc:

Pic := proc(L1,L2,col1,col2,Size:=[800,400])
local P1, P2, P3, T, P;
uses plots, plottools;
P1, P2 := Picture(L1, color=col1), Picture(L2, color=col2):
P3 := line(Sol[1..2][],color=red,thickness=3), line(Sol[3..4][],color=red,thickness=3), line(Sol[1],Sol[4],linestyle=2,thickness=3,color=red):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"]], font=[times,18], align=[left,above]);
P:=pointplot(Sol, symbol=solidcircle, color=red, symbolsize=10);
display(P1,P2,P3,T,P, scaling=constrained, size=Size, axes=none);
end proc:

Examples of use.

local D:
L1:=[[[8,0],[6,7]],[[6,7],[2,5]],[[2,5],[0,2]],[[0,2],[0,0]],[[0,0],[8,0]]]:
L2:=[[[5*cos(t)+16,5*sin(t)],t=Pi/2..Pi],[[5*cos(t)+16,5*sin(t)/2],t=Pi..2*Pi],[[21,0],[16,5]]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2): Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue",[900,400]);

The specified regions may overlap:

L1:=[[[8,0],[6,7]],[[6,7],[2,5]],[[2,5],[0,2]],[[0,2],[0,0]],[[0,0],[8,0]]]:
L2:=[[[5*cos(t)+9,5*sin(t)],t=Pi/2..Pi],[[5*cos(t)+9,5*sin(t)/2],t=Pi..2*Pi],[[14,0],[9,5]]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2):  Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue");

If there is a solution for which the cutting segments intersect the boundary of each of the regions at 2 points, then the procedure also works for such non-convex regions:

L1:=[[[cos(t),sin(t)],t=Pi/3..2*Pi-Pi/3],[[cos(-t)+1,sin(-t)],t=-Pi-Pi/3..-Pi+Pi/3]]:
L2:=[[[cos(t)+2,sin(t)],t=-Pi/6..Pi+Pi/6],[[cos(-t)+2,sin(-t)-1],t=-5*Pi/6..-Pi/6]]:
Sol:=Each_Into_2_Equal_Areas(L1,L2): Points:=[A,B,C,D]:
seq(Points[i]=Sol[i], i=1..4);
Pic(L1,L2,"Yellow","LightBlue");

A number of other interesting examples can be found in the attached file.

Each_Into_2_Equal_Areas1.mw

AI and writing Maple code

by: Harry Garst 259 Product: Maple

February 18 2025

3 0

After a long chat with ChatGPT I finally received a fully working code for a proc for a generalized Woodbury Identity for the inversion of the sum of two or more positive definite matrices.
I was inspired by a family member who is a trained professional programmer, who told me that in his professional work he uses ChatGTP for an initial draft of his program.
I did find out that ChatGTP makes errors: from simple ones like writing 'Simplify' instead of 'simplify' to serious conceptual errors, for example in recursive loops. However, ChatGTP seems to 'understand' the error after given specific feedback. Although, this does not mean that the next proposal does not contain the same logical error. But after a long chat I received a nice proc that seems to work.
My second surprise was that Gemini suggested a formula for the generalized Woodbury lemma that was unknown to me, and I was unable to find on Scholar Google or https://math.stackexchange.com. Based on a special case of that formula, I was able to write the second proc myself.
In conclusion, to start working it can be helpful to collaborate with AI friend, a little patience may help, AI may not be astute as someone on Mapleprimes wrote, but neither am I. I am now retired, and it is fun to play with Maple and AI.
By the way, the search term Woodbury did not give a single hit on Mapleprimes.With_a_little_help_from_my_friends.mw
kind regards,Harry

Teaching group theory with GroupTheory

by: epostma 1554 Product: Maple

February 11 2025

9 1

For Maple 17, in 2013, we introduced the package. It has seen a lot of improvements since its introduction, and I figured I would write something about how you can use it to teach a group theory course.

First of all, I think it would be a great idea to have your students just play with the package in Maple, and get some intuition for what makes group theory tick by getting their hands dirty. But that is not what I want to talk about today. Instead, I want to discuss how you might use it for giving your lectures - to show some visualizations while you show your beamer presentation or write on the black- or whiteboard.

When starting out with the definition of a group, you might want to compare the Cayley table of a group with that of a binary operation that doesn't define a group. A set with any binary operation is called a magma; to find good examples, we might want to use one group; one magma that has an identity and is associative but not invertible; and one magma that is has an identity and is invertible but not associative. Maybe we also want the group to not be cyclic (so the group order will have to not be prime), because those Cayley tables look a little boring. For this we can use the package.

We note that a magma that is invertible and has an identity element is called a loop. We can enumerate all magmas of a given order and with certain properties using the command , finding either the number of such objects (by default) or a list of the multiplication tables (with the option). Can we find interesting examples of order 4?

(1)

No: all loops of order 4 are groups. Let's try order 6.

(2)

Yes, we can find examples here. We find the Cayley tables of three suitable magmas; we call the non-cyclic group one , the non-associative loop one , and the non-invertible magma .

(3)

(4)

(5)

Most of these appear to have many more of one value than of another (e.g., many more threes). Let's find one where that is not so extreme. For example, we might want to minimize the standard deviation of the frequencies of the values occurring in the Cayley table.

(6)

(7)

Aha, the smallest standard deviation is a bit under two and a half, and the next possible value is greater than 3. Let's examine the Cayley tables that show this standard deviation.

(8)

(9)

The last of these looks interesting.

Now we can display these three Cayley tables:

It's clear that the multiplication shown in the third Cayley table is not invertible (because multiplying by any element other than 1, on either side, is not a permutation of the elements: rows and columns beyond the first don't have all colors). It's less visually obvious that the second Cayley table doesn't show an associative multiplication, but a quick loop showed that whereas . I'm not sure if there is something real underpinning this, but to my mind this is a nice parallel to the fact that when you're proving a set with a given operation is a group, it's easier to show that there are inverses than that the operation is associative.

A second nice visualization, also usable in the first week or two of a group theory course, is also named after Cayley: Cayley graphs. Here is a Cayley graph for the dihedral group of order 6:

(10)

Here is a well-known presentation of the alternating group on 4 letters. I find "grabbing and turning" the 3D plot shows the structure much better than simply an embedding of the graph in the plane.

(11)

Another great opportunity for a visualization comes a couple of weeks later, when you talk about subgroup lattices. Here are a few examples:

(12)

(13)

This is the plain subgroup lattice. Note how the subgroups in the first and second nontrivial "layer" from the bottom are grouped by conjugacy class. The default visualization of a subgroup lattice highlights the centre in light blue and the (in this case six) other normal subgroups in green.

We see that every conjugacy class of size one is a normal subgroup. We can also highlight other subgroups, if we want. Maybe we want to show which are the (cyclic) subgroups generated by one of the chosen generators:

(14)

Of course you can also show more advanced topics in this way. For example, here is the Frattini series of this group.

Download blogpost-algebra.mw

With two cuts into 4 equal areas

by: Kitonum 21308 Product: Maple 2018

January 02 2025

11 7

In the excellent book by W.G. Chinn, N.E. Steenrod "First Concepts of Topology" the theorem is proved which states that any bounded planar region can be cut into 4 regions of equal area by 2 perpendicular cuts (the pancake problem). This is an existence theorem which does not provide any way to find these cuts. In this post I made an attempt to find such cuts for any convex region on the plane bounded by a piecewise smooth self-non-intersecting curve.
The Into_4_Equal_Areas procedure returns a list of coordinates of 5 points: the first 4 points are the endpoints of the cutting segments, the fifth point is the intersection point of these segments. This procedure significantly uses my old procedure Area , which can be found in detail at the link https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure- . The formal argument of the Into_4_Equal_Areas procedure is a list L specifying the boundary of the region to be cut. When specifying L, the boundary can be passed clockwise or counterclockwise, but it is necessary that the parameter t (when specifying each link) should go in ascending order. This can always be achieved by replacing t with -t if necessary. The Pic procedure draws a picture of the source region and cutting segments. For ease of use, the code for the Area procedure is also provided. It is also worth noting that the procedure also works for "not too" non-convex regions (see examples below).

restart;
Area := proc(L) 
local i, var, e, e1, e2, P; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then 
P[i] := (1/2)*add(L[i, j, 1]*L[i, j+1, 2]-L[i, j, 2]*L[i, j+1, 1], j = 1 .. nops(L[i])-1) else 
var := lhs(L[i, 2]); 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := (1/2)*(int(e^2, L[i, 2])) else 
if var = y then P[i] := (1/2)*simplify(int(e-var*(diff(e, var)), L[i, 2])) else 
P[i] := (1/2)*simplify(int(var*(diff(e, var))-e, L[i, 2])) end if end if else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; 
P[i] := (1/2)*simplify(int(e1*(diff(e2, var))-e2*(diff(e1, var)), L[i, 2])) end if end if end do; 
abs(add(P[i], i = 1 .. nops(L))); 
end proc:

Into_4_Equal_Areas:=proc(L::list,N::symbol:='OneSolution', eps::numeric:=0)
local D, n, c, L1, L2, L3, f, L0, i, j, k, m, A, B, C, P, S, sol, Sol;
f:=(X,Y)->expand((y-X[2])*(Y[1]-X[1])-(x-X[1])*(Y[2]-X[2]));
L0:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L);
S:=Area(L); c:=0;
n:=nops(L);
for i from 1 to n do
for j from i to n do
for k from j to n do
for m from k to n do
if not ((nops({i,j,k})=1 and type(L[i],listlist)) or (nops({j,k,m})=1 and type(L[j],listlist)))then
A:=convert(subs(t=t1,L0[i,1]),Vector): 
B:=convert(subs(t=t2,L0[j,1]),Vector):
C:=convert(subs(t=t3,L0[k,1]),Vector): 
D:=convert(subs(t=t4,L0[m,1]),Vector):
P:=eval([x,y], solve({f(A,C),f(B,D)},{x,y})):
L1:=`if`(j=i,[subsop([2,2]=t1..t2,L0[i]),[convert(B,list),P],[P,convert(A,list)]],`if`(j=i+1,[subsop([2,2]=t1..op([2,2,2],L0[i]),L0[i]),subsop([2,2]=op([2,2,1],L0[j])..t2,L0[j]),[convert(B,list),P],[P,convert(A,list)]], [subsop([2,2]=t1..op([2,2,2],L0[i]),L0[i]),L0[i+1..j-1][],subsop([2,2]=op([2,2,1],L0[j])..t2,L0[j]),[convert(B,list),P],[P,convert(A,list)]])):
L2:=`if`(k=j,[subsop([2,2]=t2..t3,L0[j]),[convert(C,list),P],[P,convert(B,list)]],`if`(k=j+1,[subsop([2,2]=t2..op([2,2,2],L0[j]),L0[j]),subsop([2,2]=op([2,2,1],L0[k])..t3,L0[k]),[convert(C,list),P],[P,convert(B,list)]], [subsop([2,2]=t2..op([2,2,2],L0[j]),L0[j]),L0[j+1..k-1][],subsop([2,2]=op([2,2,1],L0[k])..t3,L0[k]),[convert(C,list),P],[P,convert(B,list)]])):
L3:=`if`(m=k,[subsop([2,2]=t3..t4,L0[k]),[convert(D,list),P],[P,convert(C,list)]],`if`(m=k+1,[subsop([2,2]=t3..op([2,2,2],L0[k]),L0[k]),subsop([2,2]=op([2,2,1],L0[m])..t4,L0[m]),[convert(D,list),P],[P,convert(C,list)]], [subsop([2,2]=t3..op([2,2,2],L0[k]),L0[k]),L0[k+1..m-1][],subsop([2,2]=op([2,2,1],L0[m])..t4,L0[m]),[convert(D,list),P],[P,convert(C,list)]])):
sol:=fsolve({Area(L1)-S/4,Area(L2)-S/4,Area(L3)-S/4,LinearAlgebra:-DotProduct(D-B,C-A, conjugate=false)},{t1=op([2,2,1],L0[i])-eps..op([2,2,2],L0[i])+eps,t2=op([2,2,1],L0[j])-eps..op([2,2,2],L0[j])+eps,t3=op([2,2,1],L0[k])-eps..op([2,2,2],L0[k])+eps,t4=op([2,2,1],L0[m])-eps..op([2,2,2],L0[m])+eps}) assuming real:
if type(sol,set(`=`)) then if N='OneSolution' then return convert~(eval([A,B,C,D,P],sol),list) else c:=c+1; Sol[c]:=convert~(eval([A,B,C,D,P],sol),list) fi;
 fi; fi;
od: od: od: od:
convert(Sol,list);
end proc:

Pic:=proc(L,Sol)
local P1, P2, T;
uses plots, plottools;
P1:=seq(`if`(type(s,listlist),line(s[],color=blue, thickness=2),plot([s[1][],s[2]],color=blue, thickness=2)),s=L):
P2:=line(Sol[1],Sol[3],color=red, thickness=2), line(Sol[2],Sol[4],color=red):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"],[Sol[5][],"P"]], font=[times,18], align=[left,above]);
display(P1,P2,T, scaling=constrained, size=[800,500], axes=none);
end proc:

Examples of use:

L:=[[[0,0],[1,4]],[[1,4],[6,7]],[[6,7],[12,0]],[[12,0],[0,0]]]:
Sol:=Into_4_Equal_Areas(L);
Pic(L, Sol);

# Check (areas of all 4 regions)
Area([[L[1,1],Sol[4],Sol[5],Sol[1],L[1,1]]]),
Area([[Sol[4],Sol[5],Sol[3],L[4,1],Sol[4]]]),
Area([[Sol[5],Sol[2],L[3,1],Sol[3],Sol[5]]]),
Area([[Sol[5],Sol[2],L[1,2],Sol[1],Sol[5]]]);

L:=[[[1+cos(-t),1+sin(-t)],t=-3*Pi/2..-Pi],[[0,1],[-1,0]],[[cos(t),sin(t)],t=Pi..2*Pi]]:
Sol:=Into_4_Equal_Areas(L);
Pic(L,Sol);

# The boundary is the Archimedes spiral and the arc of a circle

L:=[[[t*cos(t),t*sin(t)],t=0..2*Pi],[[Pi+5*cos(-t),sqrt(25-Pi^2)+5*sin(-t)],t=arccos(Pi/5)..Pi-arccos(Pi/5)]]:
Sol:=evalf(Into_4_Equal_Areas(L));
Pic(L,Sol);

L:=[[[0,0],[2,0]],[[2,0],[1,sqrt(3)]],[[1,sqrt(3)],[0,0]]]:
Sol:=evalf[5](Into_4_Equal_Areas(L, AllSolutions, 0.1)); # All 3 solutions
plots:-display(<Pic(L, Sol[1]) | Pic(L, Sol[2])  | Pic(L, Sol[3])>, size=[300,300]);

L:=[[[-t,-sin(-t)],t=-5*Pi/4..0],[[cos(t),sin(t)-1],t=Pi/2..3*Pi/2],[[t,cos(t)-3],t=0..3*Pi/2],[[3*Pi/2,-3],[5*Pi/4,sqrt(2)/2]]]:
Sol:=evalf(Into_4_Equal_Areas(L));
Pic(L,Sol);

More examples can be found in the attached file.

4_Equal_Area1.mw

[Edit]. The post has been edited. One inaccuracy in the code has been corrected, which could sometimes lead to errors. Two options have been added to the code of Into_4_Equal_Areas procedure. The first option is the formal argument N . If N=OneSolution (by default), the procedure returns one solution. If N=AllSolutions , the procedure returns all solutions that it can find. The eps option has also been added (by default, eps=0). It is advisable to use it when we are looking for all solutions, and the ends of the cutting segments fall on the boundaries of intervals (this option slightly expands the boundaries of intervals, otherwise the fsolve command sometimes misses solutions). Two new examples have also been added.

Winter issue of Maple Transactions now published

by: rcorless 710 Product: Maple

December 31 2024

3 2

Maple Transactions Volume 4 Issue 4 has now been published.

This issue has two Featured Contributions by people who have been plenary speakers at Maple Conferences in the past, namely Veselin Jungić and Juana Sendra. We hope you enjoy both articles. There is an accompanying video by Professor Sendra, which we will add a link to when it becomes ready.

As usual, there is an article in the Editor's Corner, but this one is a bit different. In this one, Michelle Hatzell (the new copyeditor for Maple Transactions, who is also a Masters' student working with me at Western) and I have written about a fun use of Maple's colour contour plots to make an image that might be used as the cover of an upcoming book, namely Perturbation Methods using backward error, which I'm just finishing now with Nic Fillion and which SIAM will publish next year. So, while there's some math in that paper, it's more about Maple's utilities for colour plotting; so you might find it useful. We also hope you like at least some of the images. Some are more attractive than others!

We have several Refereed Contributions, not all of which are ready at this time of publication but which will be added as they are revised and sent in. We have a nice paper on using continued fractions in a high school context, another on code generation, and another on using Digital Signal Processing in Engineering courses.

Finally we have a first publication in French, by Jalale Soussi. Actually we have the paper also in English: we chose to publish both, in our Communications section, each with links to the other. It is possible to publish in Maple Transactions solely in French, of course, but the author provided both, so why not?

Happy reading, and best wishes for 2025.

Sample Worksheet on ODEs with Maple – A Glimpse...

by: Paras31 205 Product: Maple 2024

November 19 2024

1 4

Hello Maple enthusiasts,

I am excited to share a sample worksheet on Ordinary Differential Equations (ODEs), created as part of my ongoing project—a book I am writing for undergraduate students. This book is designed to teach ODEs using Maple, offering an interactive and intuitive approach to solving differential equations.

As far as I know, there aren’t many books available in the Greek language that combine ODEs with Maple. In fact, I believe there’s only one other such resource, which highlights the lack of materials in this niche. My goal is to fill this gap by providing students and educators with a resource that is both practical and accessible, leveraging Maple's powerful capabilities to deepen understanding and simplify complex concepts.

The worksheet I’m sharing includes:

Step-by-step solutions to ODEs using Maple.
Graphical representations to visualize solutions, which I believe are invaluable for fostering comprehension.

I hope this preview sparks your interest and provides insight into the teaching style and structure of the upcoming book. I would love to hear your thoughts, feedback, or suggestions for topics you think should be included.

Solve the following differential equation

with initial condition

>

(1)

>

(2)

>

y(x) = c__1*exp((1/2)*x^2)

(3)

>	$particular_solution := dsolve({ic, ode}, y(x))$

(4)

>

(5)

>

Solve the following differential equation

diff(y(x), x) = x/y(x)
with initial condition

>

restart; with(plots); with(DEtools)

>

(6)

>

(7)

>

(8)

>

Solve the following differential equation

dy/dx = e^y/(x^2+1)
with initial condition

>

ode2 := diff(y(x), x) = exp(y(x))/(x^2+1)

diff(y(x), x) = exp(y(x))/(x^2+1)

(9)

>

(10)

>

(11)

>

Solve the following differential equation

dy/dx = y^2+y
with initial condition

>

(12)

>

(13)

>

y(x) = 1/(-1+exp(-x)*c__1)

(14)

>

(15)

>

Solve the following differential equation

y*dy/dx-x = 0
with initial condition and

>

(16)

>

(17)

>

(18)

>

Solve the following differential equation

dy/dx = 2*x+2*xy^2/y
with initial condition

>

(19)

>

(20)

>

(21)

>

(22)

>

Download separable_diff_equations.mw

Rolling an ellipse along spatial curves also on...

by: sija 165 Product: Maple 2024

November 17 2024

4 0

Hello,

I present the rolling of an ellipse along a spatial curve and along curves placed on a surface, also along a Mobius surface.

Attached are the Maple source files.

Best regards.

Source_-_mw.zip

Matt's teaching tool for triangles

by: Matt C Anderson 185 Product: Maple 13

November 13 2024

2 0

Hi all,

Hopefully this will be useful to educators who are introducing trigonometry to their students.

subjects - Maple procedures, complementary angles, and degree to radian conversion

Enjoy.

'big' Matt

cosine_teaching_tool.pdf

cosine_teaching_tool.mw

Have a nice day

Exploring Oscillators: From Simple Linear to Advanced...

by: Paras31 205 Product: Maple 2024

November 10 2024

0 0

Continuing from my previous post on "Exploring Oscillators: A Simple Pendulum Worksheet," we now dive deeper into the fascinating progression of oscillatory systems—from linear oscillations to complex, nonlinear dynamics with chaotic behaviors. This post examines how adjusting parameters in oscillatory models, such as the Duffing oscillator, reveal rich dynamics, including resonance and chaos.

LINEAR AND NON-LINEAR OSCILLATORS

Athanasios Paraskevopoulos

Imagine a block of mass m attached to a spring. The spring is fixed at one end, and the block is free to move horizontally on a frictionless surface.
When the block is displaced from its equilibrium position and released, it oscillates back and forth

According to Hooke's Law, the restoring force exerted by the spring is proportional to the displacement
from the equilibrium position and acts in the opposite direction:

Using Newton's second law,
, where is the acceleration, we get:

Since acceleration is the second derivative of displacement with respect to time, the equation becomes:

m*(diff(x, t, t)) = -k*x
Rearranging terms, we get:

diff(x, t, t)+k*x/m = 0
Letting where ω is the angular frequency, the equation simplifies to

diff(x, t, t)+omega^2*x = 0

•	Amplitude A: Maximum displacement from the equilibrium position.

•	Angular Frequency : Related to the spring constant and mass by .

•	Period T: The time for one complete oscillation, given by

•	Frequency f: The number of oscillations per unit time

In a linear simple harmonic oscillator (SHO), the Hamiltonian energy represents the total energy of the system, which is a combination of its kinetic and potential energy.

1.	Position and Momentum:

•	Let be the displacement of the oscillator from its equilibrium position.

•	Let be the momentum of the oscillator.

2.	Hamiltonian ():

•	The Hamiltonian function is expressed in terms of the kinetic energy () and the potential energy () of the system.

•	For a linear SHO, it is given by:

3.	Kinetic Energy (T):

Kinetic energy depends on the mass of the object and its velocity. In physics, it's expressed as:

where:

•	is the mass of the oscillator,

•	This equation shows that kinetic energy increases with the square of the momentum.

4.	Potential Energy ()

Potential energy in a simple harmonic oscillator arises from the restoring force of the spring. The potential energy is given by:

V = (1/2)*k*x^2

•	is the spring constant, and

•	is the displacement from the equilibrium position.

This means that the further the mass is displaced from equilibrium, the greater the potential energy stored in the system.

Suppose you have a spring with a spring constant and a mass . If the mass is displaced by from its equilibrium position and released, it will undergo SHM

>

k := 10; m := .2; omega := sqrt(k/m)

(1)

>

(2)

>

plot(T, x = 0 .. 1, title = "Period vs Initial Amplitude (x)", axes = boxed, labels = ["Initial Amplitude (x)", "Period (T)"])

>

f := omega/(2*Pi)

(3)

>

ODE__1 := diff(x(t), t, t)+omega^2*x(t) = 0; IC := x(0) = A, (D(x))(0) = 0

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

(4)

>	$sol := dsolve({IC, ODE__1}, x(t))$

(5)

>

plot_1 := subs(A = -.104, sol); plotsresult := plot([rhs(plot_1)], t = 0 .. 2, legend = ["-0.104cos( 7.07 t)"], color = [red])

>

x_t := rhs(subs(A = -.104, sol)); v_t := diff(x_t, t)

(6)

>

T := (1/2)*m*v_t^2; V := (1/2)*k*x_t^2; H := T+V

(7)

>

energy_plot := plot([eval(T), eval(V), eval(H)], t = 0 .. 2, color = [red, blue, green], legend = ["Kinetic Energy", "Potential Energy", "Total Energy"], title = "Energy Exchange in Simple Harmonic Oscillator", labels = ["Time (s)", "Energy (Joules)"])

>

directionfield := DEplot([diff(x(t), t) = v(t), diff(v(t), t) = -omega^2*x(t)], [x(t), v(t)], t = 0 .. 10, x = -2 .. 2, v = -10 .. 10, arrows = medium, title = "Direction Field for Simple Harmonic Oscillator", axes = boxed)

>	$sol1 := dsolve({ODE__1, x(0) = 1, (D(x))(0) = 0}, x(t)); sol2 := dsolve({ODE__1, x(0) = .5, (D(x))(0) = 0}, x(t)); sol3 := dsolve({ODE__1, x(0) = -1, (D(x))(0) = 0}, x(t)); sol4 := dsolve({ODE__1, x(0) = .7, (D(x))(0) = .5}, x(t))$

(8)

>

x1 := rhs(sol1); x2 := rhs(sol2); x3 := rhs(sol3); x4 := rhs(sol4); v1 := diff(x1, t); v2 := diff(x2, t); v3 := diff(x3, t); v4 := diff(x4, t)

(9)

>

Damped Oscillators:

Consider a mass-spring-damper system, where a block of mass is attached to a spring (with spring constant ) and a damper (with damping coefficient ).
The damper provides resistance to the motion, simulating friction or air resistance.

Damped Oscillator Equation:
m*(diff(x, t, t))+delta*(diff(x, t))+k*x = 0

This is in the form of a homogeneous second order differential equation and has a solution of the form

Substituting this form gives an auxiliary equation for

The roots of the quadratic auxiliary equation are

λ= (-delta+`&+-`(sqrt(delta^2-4*mk)))/(2*m)

The three resulting cases for the damped oscillator are

•	Overdamped

•	Critical damping

•	Underdamped

The period of a damped oscillator is primarily determined by the damped angular frequency `ω__d`; `ω__d` = sqrt(k/m-(delta/(2*m))^2)

Suppose a mass is attached to a spring with and a damper with a damping coefficient

>

k := 10; m := .2; `δ__c` := 2*sqrt(k*m); delta := .6; omega := sqrt(k/m); `ω__d` := sqrt(omega^2-(delta/(2*m))^2); T := 2*Pi/`ω__d`; f := `ω__d`/(2*Pi)

(10)

>

ODE__2 := m*(diff(x(t), t, t))+delta*(diff(x(t), t))+k*x(t) = 0; IC := x(0) = -.104, (D(x))(0) = 0

.2*(diff(diff(x(t), t), t))+.6*(diff(x(t), t))+10*x(t) = 0

(11)

>

x(t) = c__1*exp(-(3/2)*t)*sin((1/2)*191^(1/2)*t)+c__2*exp(-(3/2)*t)*cos((1/2)*191^(1/2)*t)

(12)

>	$sol := dsolve({IC, ODE__2}, x(t), numeric, output = listprocedure)$

>

H := proc (t) local pos, vel, kinetic, potential; pos := rhs(sol(t)[2]); vel := rhs(sol(t)[3]); kinetic := (1/2)*m*vel^2; potential := (1/2)*k*pos^2; return kinetic+potential end proc

>

timeRange := 0 .. 35; plot(H(t), t = timeRange, title = "Hamiltonian Energy vs Time", labels = ["Time (t)", "Energy (H)"])

>

for meth in SolutionMethods do meth, H(100) end do

(13)

Mass-Spring System with a Non-Linear Spring

•	Springs with large deflections are not linear and their forces can be approximated by a cubic term

•	If hardening spring (stiffer with larger deflection), if softening spring (more flexible with larger deflection).

•	So, the EOM fot a mass-spring system in free vibration will be:

m*(diff(x(t), t, t))+k*x(t)+`βx`^3 = 0

The period of oscillation will not be constant and depends on the applitude of oscillations and sign of

First case

>

restart; with(plots); m := .2; k := 10; beta := 1

(14)

>

.2*(diff(diff(x(t), t), t))+10*x(t)+x(t)^3 = 0

(15)

>

sol := dsolve({EOM, x(0) = -.104, (D(x))(0) = 0}, x(t), numeric, range = 0 .. 20, output = listprocedure)

>

H := proc (t) local pos, vel, kinetic, potential; pos := rhs(sol(t)[2]); vel := rhs(sol(t)[3]); kinetic := (1/2)*m*vel^2; potential := (1/2)*k*pos^2+(1/4)*beta*pos^4; return kinetic+potential end proc; timeRange := 0 .. 10; plot(H(t), t = timeRange, 0 .. .1, title = "Hamiltonian Energy vs Time", labels = ["Time (t)", "Energy (H)"])

proc (t) local pos, vel, kinetic, potential; pos := rhs(sol(t)[2]); vel := rhs(sol(t)[3]); kinetic := (1/2)*m*vel^2; potential := (1/2)*k*pos^2+(1/4)*beta*pos^4; return kinetic+potential end proc

In this non-linear mass-spring-damper system (without damping), the Hamiltonian is essentially the total energy of the system, as there is no energy loss.

•	The Hamiltonian remains relatively constant over time, as expected in an ideal undamped system, indicating energy conservation.

•	Small numerical variations may occur due to the non-linear terms and the precision of numerical methods used for solving the differential equations.

>

eq1 := diff(x(t), t) = v(t); eq2 := diff(v(t), t) = -k*x(t)/m-beta*x(t)^3/m

(16)

>

initial_conditions := [[x(0) = -.104, v(0) = 0]]; with(DEtools); DEplot([eq1, eq2], [x(t), v(t)], t = 0 .. 1, x = -.3 .. .3, v = -2 .. 2, initial_conditions, title = "Phase Portrait for Nonlinear Oscillator", axes = boxed, labels = ["Displacement x", "Velocity v"], linecolor = blue)