Applications, Examples and Libraries

Share your work here

Thanks to @salim-barzani, I became interested in the Laplace Adomian Decompositon method to solve partial differential equations. It produces a sum of analytical components that converge to the exact solution. Probably it is not that efficient for numerical solutions, but the possibility of finding an exact solution to nonlinear PDEs by finding a formula for the components and therefore the infinite sum is interesting.

The version here covers many of the common forms of PDEs, but is still a work in progress. The method for finding the Adomian polynomials could be improved by using one of the reccurence formulas in the literature. Extension to more PDE forms, ODEs etc are possible and I will attempt some of these before uploading an improved version to the application centre. In the meantime, feedback about bugs, usability etc. are appreciated.

Laplace Adomian Decomposition method to solve partial differential equations.
 D.A. Harrington, June 2025, v 1.16.
 Procedure LAD is in the startup code region.

Arguments:

1. 

Partial differential equation in one "time" variable and several other variables (called spatial variables here). Specified as an equation, or an expression implicitly equal to zero.
Comprising: (1) one time derivative term of order m (denoted L), (2) linear terms that are derivatives in the spatial variables (R), (3) nonlinear terms with derivatives in the spatial variables (N). Together, (1)-(3) is a multivariate polynomial in the dependent variable, its derivatives, and any parameters. (4) inhomogenous terms in the time and spatial variables that do not contain the dependent variable or its derivatives (G). The general form is L = R+N-G.
Any symbolic parameters are assumed to be real. Any numeric coefficients or parameters should be of type algebraic. Floating point values will be converted to rationals with convert(value, rational).

2. 

Function to solve for, e.g., u(x, y, t).

3. 

m initial condition(s) (at time zero). For m > 1, these must be in a list, and are the values of the function and its successive derivatives with respect to time, evaluated at time zero.

4. 

name of time variable.

5. 

number of iterations.

6. 

(optional) order for a fractional derivative, specified as fracorder = name or fracorder = numeric*value (floating point values are converted to type rational). The permissible values alpha are related to the value of m: m-1 < alpha and alpha <= m, i.e., an mth order time derivative is specified in the pde and m initial conditions are provided, as well as the value of alpha in the correct range. A symbolic alpha is assumed to be in this range.

infolevel[LAD] may be set to different values to print out additional information as follows. Greater numbers include the information from smaller values.

1. 

The nonlinear expansion variables (those in the Adomian polynomials).

2. 

The different parts of the the pde (L, R, N, G) in jet notation.

3. 

Progress is indicated, as the time at each iteration.

4. 

Values of the components of the solution as they are produced (one per iteration).

restart

It may be useful to load the Physics package if derivatives contain abs. LAD converts these to a form without abs but with the conjugate of the dependent variable. The Physics package affects the internal simplifications and may affect the form of the output or the running time.

PDEtools:-declare(U(x, t), quiet)

Example from A-M. Wazwaz, Applied Mathematics and Computation 111 (2000) 53. Sec. 4.

pde := diff(U(x, t), t)+U(x, t)^2*(diff(U(x, t), x)); inx := 3*x

diff(U(x, t), t)+U(x, t)^2*(diff(U(x, t), x))

3*x

infolevel[LAD] := 2

approx := LAD(pde, U(x, t), inx, t, 10)

LAD: L = U[t]; R = 0; N = -U^2*U[x]; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: [U, U[x]]

175692092397588*t^10*x^11-5650915252554*t^9*x^10+184670433090*t^8*x^9-6155681103*t^7*x^8+210450636*t^6*x^7-7440174*t^5*x^6+275562*t^4*x^5-10935*t^3*x^4+486*t^2*x^3-27*t*x^2+3*x

By extrapolating to an infinite number of terms, we can find an exact solution. (Wazwaz calculates the coefficient of x^4*t^3 incorrectly and deduces an incorrect exact solution.) Multiplying by t gives a series in y = x*t for which guessgf can use the coefficients to find the exact solution.

algsubs(x*t = y, expand(approx*t)); [seq(coeff(%, y, i), i = 0 .. degree(%, y))]; gfun:-guessgf(%, y)[1]; exact := (eval(%, y = x*t))/t

175692092397588*y^11-5650915252554*y^10+184670433090*y^9-6155681103*y^8+210450636*y^7-7440174*y^6+275562*y^5-10935*y^4+486*y^3-27*y^2+3*y

[0, 3, -27, 486, -10935, 275562, -7440174, 210450636, -6155681103, 184670433090, -5650915252554, 175692092397588]

6*y/(1+(1+36*y)^(1/2))

6*x/(1+(36*t*x+1)^(1/2))

Check

pdetest(U(x, t) = exact, [pde, U(x, 0) = inx])

[0, 0]

An example with inhomogeneous terms (but no nonlinear part). Ex. 3 from R. Shah et al, Entropy 21 (2019) 335.

pde := diff(U(x, t), t)+diff(U(x, t), `$`(x, 3)) = -sin(Pi*x)*sin(t)-Pi^3*cos(Pi*x)*cos(t); inx := sin(Pi*x); exact := sin(Pi*x)*cos(t); pdetest(U(x, t) = exact, [pde, U(x, 0) = inx])

diff(U(x, t), t)+diff(diff(diff(U(x, t), x), x), x) = -sin(Pi*x)*sin(t)-Pi^3*cos(Pi*x)*cos(t)

sin(Pi*x)

sin(Pi*x)*cos(t)

[0, 0]

approx := LAD(pde, U(x, t), inx, t, 10)

LAD: L = U[t]; R = -U[x,x,x]; N = 0; G = -sin(Pi*x)*sin(t)-Pi^3*cos(Pi*x)*cos(t).

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: []

sin(Pi*x)-cos(Pi*x)*Pi^3*sin(t)+sin(Pi*x)*(-1+cos(t))-Pi^3*(Pi^3*(-1+cos(t))*sin(Pi*x)-cos(Pi*x)*sin(t))+(Pi^3*(-sin(t)+t)*cos(Pi*x)+sin(Pi*x)*(-1+cos(t)))*Pi^6-(1/2)*Pi^9*(Pi^3*(t^2+2*cos(t)-2)*sin(Pi*x)+2*cos(Pi*x)*(-sin(t)+t))-(1/6)*Pi^12*(Pi^3*(t^3+6*sin(t)-6*t)*cos(Pi*x)-3*sin(Pi*x)*(t^2+2*cos(t)-2))+(1/24)*(Pi^3*(t^4-12*t^2-24*cos(t)+24)*sin(Pi*x)+4*cos(Pi*x)*(t^3+6*sin(t)-6*t))*Pi^15+(1/120)*Pi^18*(Pi^3*(t^5-20*t^3-120*sin(t)+120*t)*cos(Pi*x)-5*sin(Pi*x)*(t^4-12*t^2-24*cos(t)+24))-(1/720)*Pi^21*(Pi^3*(t^6-30*t^4+360*t^2+720*cos(t)-720)*sin(Pi*x)+6*cos(Pi*x)*(t^5-20*t^3-120*sin(t)+120*t))-(1/5040)*(Pi^3*(t^7-42*t^5+840*t^3+5040*sin(t)-5040*t)*cos(Pi*x)-7*sin(Pi*x)*(t^6-30*t^4+360*t^2+720*cos(t)-720))*Pi^24+(1/40320)*Pi^27*(Pi^3*(t^8-56*t^6+1680*t^4-20160*t^2-40320*cos(t)+40320)*sin(Pi*x)+8*cos(Pi*x)*(t^7-42*t^5+840*t^3+5040*sin(t)-5040*t))+(1/362880)*Pi^30*(Pi^3*(t^9-72*t^7+3024*t^5-60480*t^3-362880*sin(t)+362880*t)*cos(Pi*x)-9*sin(Pi*x)*(t^8-56*t^6+1680*t^4-20160*t^2-40320*cos(t)+40320))

Expand the exact solution as a series in t to the same order and verify that it is the same as above.

series(exact-approx, t, 11, oterm = false)

0

An example with a complex solution

pde := I*(diff(U(x, t), t))+diff(U(x, t), `$`(x, 2))+2*abs(U(x, t))^2*U(x, t); inx := sech(x); exact := sech(x)*exp(I*t); `assuming`([simplify(pdetest(U(x, t) = exact, [pde, U(x, 0) = inx]))], [real])

I*(diff(U(x, t), t))+diff(diff(U(x, t), x), x)+2*abs(U(x, t))^2*U(x, t)

sech(x)

sech(x)*exp(I*t)

[0, 0]

approx := LAD(pde, U(x, t), inx, t, 10)

LAD: L = I*U[t]; R = -U[x,x]; N = -2*U^2*C; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: [C, U]

sech(x)+I*sech(x)*t-(1/2)*sech(x)*t^2-((1/6)*I)*sech(x)*t^3+(1/24)*sech(x)*t^4+((1/120)*I)*sech(x)*t^5-(1/720)*sech(x)*t^6-((1/5040)*I)*sech(x)*t^7+(1/40320)*sech(x)*t^8+((1/362880)*I)*sech(x)*t^9-(1/3628800)*sech(x)*t^10

The series for exp(I*t) is apparent.

collect(approx, sech); series(approx-exact, t, 11, oterm = false)

(1+I*t-(1/2)*t^2-((1/6)*I)*t^3+(1/24)*t^4+((1/120)*I)*t^5-(1/720)*t^6-((1/5040)*I)*t^7+(1/40320)*t^8+((1/362880)*I)*t^9-(1/3628800)*t^10)*sech(x)

0

An example with fractional order differentiation wrt t of order alpha. For m-1 < alpha and alpha <= m, the derivative must be entered with order m and there must be a list of m initial conditions:f(x, 0), (D[2](f))(x, 0), (D[2, 2](f))(x, 0), () .. ()

 Ex. 1 from R. Shah et al, Entropy 21 (2019) 335

pde := diff(U(x, t), t) = -2*(diff(U(x, t), x))-(diff(U(x, t), `$`(x, 3))); inx := sin(x)

diff(U(x, t), t) = -2*(diff(U(x, t), x))-(diff(diff(diff(U(x, t), x), x), x))

sin(x)

Symbolic order is accepted.

approx := LAD(pde, U(x, t), inx, t, 6, fracorder = alpha)

LAD: L = U[t]; R = -2*U[x]-U[x,x,x]; N = 0; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: []

sin(x)-cos(x)*t^alpha/GAMMA(1+alpha)-sin(x)*t^(2*alpha)/GAMMA(1+2*alpha)+cos(x)*t^(3*alpha)/GAMMA(1+3*alpha)+sin(x)*t^(4*alpha)/GAMMA(1+4*alpha)-cos(x)*t^(5*alpha)/GAMMA(1+5*alpha)-sin(x)*t^(6*alpha)/GAMMA(1+6*alpha)

Giving a numerical value for the order will generally be more efficient. For alpha = 1/2, we can find an exact solution.

approx := collect(LAD(pde, U(x, t), inx, t, 20, fracorder = 1/2), [sin, cos])

LAD: L = U[t]; R = -2*U[x]-U[x,x,x]; N = 0; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: []

(1-t+(1/2)*t^2-(1/6)*t^3+(1/24)*t^4-(1/120)*t^5+(1/720)*t^6-(1/5040)*t^7+(1/40320)*t^8-(1/362880)*t^9+(1/3628800)*t^10)*sin(x)+(-2*t^(1/2)/Pi^(1/2)+(4/3)*t^(3/2)/Pi^(1/2)-(8/15)*t^(5/2)/Pi^(1/2)+(16/105)*t^(7/2)/Pi^(1/2)-(32/945)*t^(9/2)/Pi^(1/2)+(64/10395)*t^(11/2)/Pi^(1/2)-(128/135135)*t^(13/2)/Pi^(1/2)+(256/2027025)*t^(15/2)/Pi^(1/2)-(512/34459425)*t^(17/2)/Pi^(1/2)+(1024/654729075)*t^(19/2)/Pi^(1/2))*cos(x)

The first series is exp(-t), the second series is not so obvious, After some scaling, the second series may be passed to guessgf

ser2 := `assuming`([simplify(expand(sqrt(Pi)*op(2, approx)/(sqrt(t)*cos(x))))], [t > 0]); [seq(coeff(ser2, t, i), i = 0 .. 9)]; ser2ex := gfun:-guessgf(%, t)[1]

-2+(4/3)*t-(8/15)*t^2+(16/105)*t^3-(32/945)*t^4+(64/10395)*t^5-(128/135135)*t^6+(256/2027025)*t^7-(512/34459425)*t^8+(1024/654729075)*t^9

-exp(-t)*Pi^(1/2)*erfi(t^(1/2))/t^(1/2)

So the full solution is

exact := exp(-t)*sin(x)+ser2ex*sqrt(t)*cos(x)/sqrt(Pi)

exp(-t)*sin(x)-exp(-t)*erfi(t^(1/2))*cos(x)

Check. fracdiff's default method returns an integral, method = laplace fails, but method = series works with cancellation of all terms but the "left over" last half-integral-order one.

Ord := 20; fracdiff(exact, t, 1/2, method = series, methodoptions = [order = Ord])-series(eval(rhs(pde), U(x, t) = exact), t, Ord)

20

-(1048576/319830986772877770815625)*sin(x)*t^(39/2)/Pi^(1/2)+O(t^(39/2))-O(t^20)

A second-order example (linear)

pde := diff(U(x, t), t, t) = -c^2*(diff(U(x, t), `$`(x, 2))); inx := [exp(I*x/c), exp(I*x/c)]; exact := exp(t+I*x/c)

diff(diff(U(x, t), t), t) = -c^2*(diff(diff(U(x, t), x), x))

[exp(I*x/c), exp(I*x/c)]

exp(t+I*x/c)

pdetest(U(x, t) = exact, [pde, U(x, 0) = inx[1], (D[2](U))(x, 0) = inx[2]])

[0, 0, 0]

approx := LAD(pde, U(x, t), inx, t, 6)

LAD: L = U[t,t]; R = -c^2*U[x,x]; N = 0; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: []

exp(I*x/c)*(1+t)+(1/6)*exp(I*x/c)*t^2*(3+t)+(1/120)*exp(I*x/c)*t^4*(t+5)+(1/5040)*exp(I*x/c)*t^6*(7+t)+(1/362880)*exp(I*x/c)*t^8*(9+t)+(1/39916800)*exp(I*x/c)*t^10*(t+11)+(1/6227020800)*exp(I*x/c)*t^12*(13+t)

series(exact-approx, t, 13, oterm = false)

0

Fractional order example

approx := collect(LAD(pde, U(x, t), inx, t, 10, fracorder = 3/2), [exp, t])

LAD: L = U[t,t]; R = -c^2*U[x,x]; N = 0; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: []

(1+(1/20922789888000)*t^16+(1/1307674368000)*t^15+(32768/6190283353629375)*t^(29/2)/Pi^(1/2)+(16384/213458046676875)*t^(27/2)/Pi^(1/2)+(1/6227020800)*t^13+(1/479001600)*t^12+(4096/316234143225)*t^(23/2)/Pi^(1/2)+(2048/13749310575)*t^(21/2)/Pi^(1/2)+(1/3628800)*t^10+(1/362880)*t^9+(512/34459425)*t^(17/2)/Pi^(1/2)+(256/2027025)*t^(15/2)/Pi^(1/2)+(1/5040)*t^7+(1/720)*t^6+(64/10395)*t^(11/2)/Pi^(1/2)+(32/945)*t^(9/2)/Pi^(1/2)+(1/24)*t^4+(1/6)*t^3+(8/15)*t^(5/2)/Pi^(1/2)+(4/3)*t^(3/2)/Pi^(1/2)+t)*exp(I*x/c)

There appear to be two separate series in t, with integer and half-integer powers.

t_series := approx/exp(I*x/c); t1, t2 := selectremove(proc (w) options operator, arrow; (degree(w, t))::integer end proc, t_series)

1+(1/20922789888000)*t^16+(1/1307674368000)*t^15+(1/6227020800)*t^13+(1/479001600)*t^12+(1/3628800)*t^10+(1/362880)*t^9+(1/5040)*t^7+(1/720)*t^6+(1/24)*t^4+(1/6)*t^3+t, (32768/6190283353629375)*t^(29/2)/Pi^(1/2)+(16384/213458046676875)*t^(27/2)/Pi^(1/2)+(4096/316234143225)*t^(23/2)/Pi^(1/2)+(2048/13749310575)*t^(21/2)/Pi^(1/2)+(512/34459425)*t^(17/2)/Pi^(1/2)+(256/2027025)*t^(15/2)/Pi^(1/2)+(64/10395)*t^(11/2)/Pi^(1/2)+(32/945)*t^(9/2)/Pi^(1/2)+(8/15)*t^(5/2)/Pi^(1/2)+(4/3)*t^(3/2)/Pi^(1/2)

The integer-powers series is the series for exp(t) with every third term missing: missing powers 2, 5, 8, 11, ...

t1sum := simplify(exp(t)-(sum(t^(3*i+2)/factorial(3*i+2), i = 0 .. infinity))); series(%, t, 17)

(2/3)*exp(t)+(1/3)*(cos((1/2)*3^(1/2)*t)+3^(1/2)*sin((1/2)*3^(1/2)*t))*exp(-(1/2)*t)

series(1+t+(1/6)*t^3+(1/24)*t^4+(1/720)*t^6+(1/5040)*t^7+(1/362880)*t^9+(1/3628800)*t^10+(1/479001600)*t^12+(1/6227020800)*t^13+(1/1307674368000)*t^15+(1/20922789888000)*t^16+O(t^17),t,17)

And for the other one, the following scaled series shows coefficients with powers of two divided by double factorials, every third term missing: missing powers 2,5,8,11,...

expand(sqrt(Pi)*t2/(4*t^(3/2))); all := add(2^i*t^i/doublefactorial(2*(i+1)+1), i = 0 .. 9); missing := add(eval(2^i*t^i/doublefactorial(2*(i+1)+1), i = 3*j+2), j = 0 .. 3)

(8192/6190283353629375)*t^13+(4096/213458046676875)*t^12+(1024/316234143225)*t^10+(512/13749310575)*t^9+(128/34459425)*t^7+(64/2027025)*t^6+(16/10395)*t^4+(8/945)*t^3+(2/15)*t+1/3

1/3+(2/15)*t+(4/105)*t^2+(8/945)*t^3+(16/10395)*t^4+(32/135135)*t^5+(64/2027025)*t^6+(128/34459425)*t^7+(256/654729075)*t^8+(512/13749310575)*t^9

(4/105)*t^2+(32/135135)*t^5+(256/654729075)*t^8+(2048/7905853580625)*t^11

t2sum := 4*t^(3/2)*(sum(2^i*t^i/doublefactorial(2*(i+1)+1), i = 0 .. infinity)-(sum(eval(2^i*t^i/doublefactorial(2*(i+1)+1), i = 3*j+2), j = 0 .. infinity)))/sqrt(Pi)

4*t^(3/2)*((1/4)*(Pi^(1/2)*exp(t)-2*t^(1/2)-Pi^(1/2)*erfc(t^(1/2))*exp(t))/t^(3/2)-(4/105)*t^2*hypergeom([1], [3/2, 11/6, 13/6], (1/27)*t^3))/Pi^(1/2)

Putting it together suggests

exact := exp(I*x/c)*simplify(t1sum+t2sum)

(1/3)*exp(I*x/c)*(Pi^(1/2)*exp(-(1/2)*t)*sin((1/2)*3^(1/2)*t)*3^(1/2)+Pi^(1/2)*exp(-(1/2)*t)*cos((1/2)*3^(1/2)*t)+3*Pi^(1/2)*exp(t)*erf(t^(1/2))-(16/35)*t^(7/2)*hypergeom([1], [3/2, 11/6, 13/6], (1/27)*t^3)+2*Pi^(1/2)*exp(t)-6*t^(1/2))/Pi^(1/2)

Check.

Ord := 20; fracdiff(exact, t, 3/2, method = series, methodoptions = [order = Ord+1])-series(eval(rhs(pde), U(x, t) = exact), t, Ord)

20

O(t^(39/2))-(1048576/319830986772877770815625)*exp(I*x/c)*t^(39/2)/Pi^(1/2)-O(t^(41/2))

Numerical calculation of a soliton

pde := I*(diff(U(x, t), t))+diff(U(x, t), `$`(x, 2))+U(x, t)^2*conjugate(U(x, t)); exact := sqrt(2*B)*sech(sqrt(B)*(2*k*t-x-x__0))*exp(I*(k*x+(-k^2+B)*t)); `assuming`([simplify(pdetest(U(x, t) = exact, pde))], [real, B > 0]); inx := eval(exact, t = 0)

I*(diff(U(x, t), t))+diff(diff(U(x, t), x), x)+U(x, t)^2*conjugate(U(x, t))

2^(1/2)*B^(1/2)*sech(B^(1/2)*(2*k*t-x-x__0))*exp(I*(k*x+(-k^2+B)*t))

0

2^(1/2)*B^(1/2)*sech(B^(1/2)*(-x-x__0))*exp(I*k*x)

With numerical parameters it runs faster. Float parameters will be converted to rationals.

params := {B = 2, k = -1/2, x__0 = -2}; inxnum := eval(inx, params); exactnum := eval(exact, params)

2*sech(2^(1/2)*(-x+2))*exp(-((1/2)*I)*x)

2*sech(2^(1/2)*(-t-x+2))*exp(I*(-(1/2)*x+(7/4)*t))

infolevel[LAD] := 3; approx := LAD(pde, U(x, t), inxnum, t, 28)

LAD: L = I*U[t]; R = -U[x,x]; N = -U^2*C; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: [C, U]

LAD: calculating component 1 at time 0.

LAD: calculating component 2 at time .221

[Edited for brevity]

LAD: calculating component 27 at time 90.414

LAD: calculating component 28 at time 103.870

LAD: completed at time 105.309

For larger x and t, the approximation diverges at a point that depends on the number of iterations.

plot3d([abs(exactnum), abs(approx)], t = 0 .. 2, x = -3 .. 3, color = [red, blue], view = [default, default, 0 .. 5], orientation = [-30, 60])

NULL

Download Laplace_Adomian_Decomposition_procedure8b.mw

 

The link below goes to the Proceedings of the Maple 2024 Conference, which includes several articles that will be of interest to the readers of Maple Primes.

There may be one more paper coming in to the proceedings later as per policy; since most things are ready, away we go!

Proceedings of the Maple 2024 Conferenc

In Optical-Illusion--Impossible-Prisms, mcarvalho provides a Maple Learn worksheet to illustrate an intriguing optical illusion.  Here I do that illustration in Maple, and in 3D! 

The graphics below shows 4 bars.  Or is it 3 bars?

Download the worksheet, stare at the graphics for short while to grasp the nature of the optical illusion, and then rotate the 3D graph with the mouse and see the illusion fall apart.

optical-illusion.mw

Geometrical entertainment in the form of rolling without slipping, now inside a torus.
You can also print out the corresponding equations - these graphs are # in the text, but it is better to do this separately from the geometric animation, because textplot3d takes up a lot of resources.

Please consider it not as a Maple program, but simply as an idea for a corresponding algorithm.
for_TORUS_IN_TORUS_for.mw
 

The workbook includes the excel data file but I also included the worksheet and data file separate for versions that can't load workbooks.

 

Average_World_Temperature_workbook-.maple

Average_World_Temperature_since_1850_worksheet.mw

HadCRUT5.zip

Last week, DeepMind announced their new AlphaEvolve software along with a number of new results improving optimization problems in a number of areas. One of those results is a new, smaller, tensor decomposition for multiplying two 4 × 4 matrices in 48 scalar multiplications, improving on the previous best decomposition which needed 49 multiplications (Strassen 1969). While the DeepMind paper was mostly careful about stating this, when writing the introduction, the authors mistakenly/confusingly wrote:

For 56 years, designing an algorithm with rank less than 49 over any field with characteristic 0 was an open problem. AlphaEvolve is the first method to find a rank-48 algorithm to multiply two 4 × 4 complex-valued matrices

This gives the impression that this new result is the way to multiply 4 × 4 matrices over a field with the fewest number of field multiplications. However, it's not even the case that Strassen's 1969 paper gave a faster way to multiply 4 × 4 matrices. In fact, Winograd's 1968 paper on computing inner products faster could be applied to 4 × 4 matrix multiplication to get a formula using only 48 multiplications. Winograd uses a trick 

which changes two multiplications of ab's into one multiplication of ab's and two multiplications involving only a's or only b's. The latter multiplications can be pre-computed and saved when calculating all the innerproducts in a matrix multiplication

    # i=1..4, 4*2 = 8 multiplications
    p[i] := -A[i, 1]*A[i, 2] - A[i, 3]*A[i, 4];
    # 4*2 = 8 multiplications
    q[i] := -B[1, i]*B[2, i] - B[3, i]*B[4, i];
    # i=1..4, j=1..4, 4*4*2 = 32 multiplications
    C[i,j] = p[i] + q[j] 
            + (A[i,1]+B[2,j])*(A[i,2]+B[1,j])
            + (A[i,3]+B[4,j])*(A[i,4]+B[3,j])

It is simple to verify that C[i,j] = A[i,..] . B[..,j] and that the above formula calculates all of the entries of C=A.B with 8+8+32=48 multiplications.

So, if Winograd achieved a formula of 48 multiplications in 1968, why is the DeepMind result still interesting? Well, the relative shortcoming of Winograd's method is that it only works if the entries of A and B commute, while tensor decomposition formulas for matrix multiplication work even over noncommutative rings. That seems very esoteric, but everyone's favorite noncommutative ring is the ring of n × n matrices. So if a formula applies to a matrix of matrices, that means it gives a recursive formula for matrix multiplication - an tensor decomposition formulas do just that.

The original 1969 Strassen result gave a tensor decomposition of 2 × 2 matrix multiplication using only 7 scalar multiplications (rather than the 8 required by the inner product method) and leads to a recursive matrix multiplication algorithm using O(N^log[2](7)) scalar multiplications. Since then, it has been proved that 7 is the smallest rank tensor decompostion for the 2 × 2 case and so researchers have been interested in the 3x3 and larger cases. Alexandre Sedoglavic at the University of Lille catalogs the best tensor decompositions at https://fmm.univ-lille.fr/ with links to a Maple file with an evaluation formula for each.

The previous best 4 × 4 tensor decomposition was based on using the 2 × 2 Strassen decomposition recursively (2 × 2 matrix of 2 × 2 matrices) which led to 7 2 × 2 matrix multiplications each requiring 7 scalar multiplications, for a total of 49. The new DeepMind result reduces that to 48 scalar multiplications, which leads to a recursive algorithm using O(Nlog[4](48)) scalar multiplications: O(N2.7925) vs. O(N2.8074) for Strassen. This is a theoretical improvment over Strassen but in 2024 the best known multiplication algorithm has much lower complexity: O(N2.3713) (see https://arxiv.org/abs/2404.16349). Now, there might be some chance that the DeepMind result could be used in practical implementations, but its reduction in the number of multiplications comes at the cost of many more additions. Before doing any optimizations I counted 1264 additions in the DeepMind tensor, compared to 318 in the unoptimized 4×4 Strassen tensor (which can be optimized to 12*7+4*12=132). Finally, the DeepMind tensor decomposition involves 1/2 and sqrt(-1), so it cannot be used for fields of characteristic 2, or fields without sqrt(-1). Of course, the restriction on characteristic 2 is not a big deal, since the DeepMind team found a 4 × 4 tensor decomposition of rank 47 over GF2 in 2023 (see https://github.com/google-deepmind/alphatensor).

Now if one is only interested in 4 × 4 matrices over a commutative ring, the Winograd result isn't even the best possible. In 1970, Waksman https://ieeexplore.ieee.org/document/1671519 demonstrated an improvement of Winograd's method that used only 46 multiplications as long as the ring allowed division by 2. Waksman's method has since been improved by Rosowski to remove the divisions see https://arxiv.org/abs/1904.07683. Here's that nice compact straight-line program that computes C = A · B in 46 multiplications.

And, here attached are five Maple worksheets that create the explicit formulas for the 4 × 4 matrix methods mentioned in this post, and verify their operation count, and that they are correct.

Strassen_444.mw  Winograd.mw  DM_444.mw Waksman.mw Rosowski.mw

If you are interested in some of the 2022 results, I posted some code to my GitHub account for creating formulas from tensor decompositions, and verifying them on symbolic matrix multiplications: https://github.com/johnpmay/MapleSnippets/tree/main/FastMatrixMultiplication

The clear followup question to all of this is:  does Maple use any of these fast multiplication schemes?  And the answer to that is mostly hidden in low level BLAS code, but generally the answer is: No, the straightforward inner-product scheme for multiplying matrices optimized memory access thus usually ends up being the fastest choice in most cases.

 

For a long time, triggered by a disagreement with one of my teachers, I wanted to demonstrate that Euler equations are not absolutely necessary to reproduce gyroscopic effects. Back then, there were no computer tools like Maple or programming languages with powerful libraries like Python. Propper calculations by hand (combining Newton’s equations and vector calculus) would have required days without guarantee of immediate success. Overall, costs and risks were too high to go into an academic argument with someone in charge of grading students.

Some years ago, I remembered the unfinished discussion with my teacher and simulated with MapleSim the simple gyroscope with two point masses that I had in mind at the time. It took only 10 minutes to demonstrate that I was right. At least I thought so. As I discovered recently when investigating the intermediate axis theorem, MapleSim derives behind the scenes Euler equations. This devalued the demonstration.

This post is about a second and successful attempt of a demonstration with Maple employing Lagrangian mechanics.  A rotating system of three point masses connected by rigid struts is used. The animation below from the attached Maple worksheet exactly reproduces a simulation of an equivalent T-shaped structure of three identical masses presented here.

 

Lagrangian Mechanics

The worksheet uses Lagrangian mechanics to derive equations of motion. Only translational energy terms are used in the Lagrangian to prevent Euler equations from being derived. To account for the bound motion of the three point masses, geometric constraints with Lagrange multipliers were added to the Lagrangian L of the system. This lead to a modified Lagrangian Lthat can be used with dedicated Maple commands to derive with little effort a set of Lagrange’s equations of the first kind and the corresponding constraints

Ein Bild, das Text, Schrift, weiß, Reihe enthält.

KI-generierte Inhalte können fehlerhaft sein.

     Ein Bild, das Schrift, Handschrift, Text, Typografie enthält.

KI-generierte Inhalte können fehlerhaft sein.

(source https://en.wikipedia.org/wiki/Lagrangian_mechanics)

For the system of three point masses the above equations lead to 9 coupled second-order ordinary differential equations (ODEs) and 3 algebraic constraints 

(Maple output from the command Physics,LagrangeEquations)

where xi, yi and zi are the components of the position vectors i of the masses 1 to 3 and the li,j are the constraints between the masses i and j, and b and h are the base and the height of the triangle.

The 12 equations together are also referred to as differential algebraic equations (DAEs). Maple has dedicated solvers for such systems which make implementation easy. The most difficult part is setting the initial conditions for all point masses. In this respect MapleSim is even easier to use since not all initial conditions have to be exactly defined by the user. MapleSim also detects constraints that allow for a simplification of the problem by reducing the number of variables to solve. This leads automatically to 3 instead of 12 equations to be solved. Computational effort is reduced in this way significantly.

 

Newtonian Mechanics

One could argue now that the above is a demonstration with Lagrangian mechanics and not with Newtonian mechanics. To treat the system in a Newtonian way, the masses must be isolated and internal forces acting on each mass via the struts are applied to each mass and effectively become external forces. This leads to 9 ODEs with 27 unknows

Ein Bild, das Text, Schrift, Screenshot enthält.

KI-generierte Inhalte können fehlerhaft sein.

Following actio=reactio for each of the (massless) struts reduces the number of unknows to 18

Ein Bild, das Text, Handschrift, Schrift enthält.

KI-generierte Inhalte können fehlerhaft sein.

To solve for the 18 unknows, 9 more relations are required. 3 algebraic constraints that keep the distance between the masses constant have already been listed in the previous section. 6 further algebraic constraints can be established from the fact that the force vectors point towards the opposing mass (see also below).
The effort to solve this system of equations will be even greater but with the benefit of having information about the internal forces in the system.
Before making this effort, it is advisable to take a closer look at the equations of motion derived so far.

 

Forces and the "mysterious" Lagrange Multipliers

Rearranging equations of motion from Lagrangian mechanics to

Ein Bild, das Text, Schrift, Screenshot, Electric Blue (Farbe) enthält.

KI-generierte Inhalte können fehlerhaft sein.

and comparing this to the equations of motion from Newtonian mechanics yields in vector notation 

or more general for the forces

Ein Bild, das Schrift, Handschrift, Reihe, Typografie enthält.

KI-generierte Inhalte können fehlerhaft sein.

where the i are the position vectors of the individual masses and the  are the constraint forces between them.  

In the case of a triangle formed by struts all internal forces must act in the direction of the edges of the triangle. If they would not act this way, opposing pairs of  and   would create a torque around the struts which would lead to an infinite angular acceleration of the massless struts. The above equation confirms this reasoning: The internal forces act in the direction of the difference vectors between the position vectors of the masses (which describe the edges) and scale with lij.

The beauty of the Lagrange multipliers in this case is that they hit 3 birds (three components of the vectors ) with one stone. This reduces the number of unknowns.

However, the Lagrange multipliers are somehow mysterious because they do not represent a physical quantity, but they can be used to calculate meaningful and correct physical results.

What makes them even more mysterious is the fact that positions constraints can be expressed in different ways. In the above example the square of the distance between the masses is kept constant instead of the distance. There are many more possibilities to formulate the constraint of constant distance and each of them results in different multipliers lij with different units. In principle they should all work equivalently but might not all be usable with dedicated solvers.

According to the above equation, the internal forces in a strut scale with the Lagrange multipliers and the length of the strut. During the back and forth flip of the triangle in the above animation the forces vary which can be appreciated from the lij in the following plot. 

Ein Bild, das Diagramm, Reihe enthält.

KI-generierte Inhalte können fehlerhaft sein.

Some observations:

  • At the start, only the strut between the red and the green masses is tensed by centrifugal forces. This would have been intuitively expected.
  • At the start, the broken symmetry in the initial conditions is already visible by the imbalance of the forces in the two other struts.
  • At no time two forces are zero
  • There is never compression in two struts at the same time. The existence of compression forces renders any attempt to replace the struts by cables useless
  • Plotted together in 3D the Lagrange multipliers describe a seemingly perfect circle.

Ein Bild, das Diagramm, Reihe, Kreis, Origami enthält.

KI-generierte Inhalte können fehlerhaft sein.

The last observation in particular shows that both the Lagrange multipliers and the IAT still hold some secrets that need to be clarified:

  • Is it an exact circle? How to prove that?
  • Does a circle appear for all initial conditions and geometries (obtuse isosceles triangles) when the intermediate axis theorem manifests?
  • What does the circle represent? Is it a kind of an invariant between the Lagrange multipliers that allows the calculation of one force when the two others are known?
  • Is there an analytical representation of the circle as a function of the geometry and the initial conditions?
  • What determines the center, the radius and the orientation of the circle?

 

Conclusion

Overall, the approach of adding non-physical terms that are zero to a physical quantity (the difference between kinetic and potential energy) to derive something meaningful is not obvious at all and underlines the genius of Lagrange. For a system of three bound masses it could be used to calculate internal forces as opposed to the more common use of calculating external constraint forces. Beyond the lack of fully satisfying intuitive explanations, the IAT still offers unanswered scientific questions.

 

PS.: 

  • The exercise was a nice cross-check of MapleSim and Maple.
  • dsolve and odeplot are awsome

 

IAT_without_Euler_equations.mw

My maplet (author: Xavier Cormier) "sudokuvision_multithread3" allows you to generate, solve, and play nxm region sudoku puzzles. It is possible to import lists of puzzles in .csv, .sdm, and .sdk format. It loads a grid in line format. It saves puzzles in .txt and .csv format to import them on Android (Andoku3, OpenSudoku; you must re-save the file in .sdm) and on iOS (SmartSudoku, Good Sudoku). It saves game puzzles (in numbers or colors) and their solution as a .gif image.
Two possible interfaces.
Multithreading is used to generate sudoku puzzles in parallel on multiple threads: we use sequences for the random number and the number of sudokus for each thread:
e.g., [4521,7458]=s1 for the random number and [2,2]=s2 for the number of sudokus. We then generate 2[i] sudokus for the ith thread. with s1[i] as a random integer.
The maplet is compatible with Maple 2021 but not with Maple 2015.

Go to my Maple Cloud link :

https://maple.cloud/app/6283420531032064/sudoku+multithread3?key=8DBF1E6D71EF423B9BE883C27FF0FD892C14B4FCD25B4920A29AFE4E4B4363B8

My maplet (author :Xavier Cormier) studies the partitioning of a set and the decomposition of an integer into a sum.
It can be used in the case of a dice roll.

It determines:
- the number of times each face is obtained for n rolls for an m-sided die.
- the decomposition of an integer nx into the sum of mx integers.
- the partitioning of a set into u non-empty, disjoint pairwise parts.
- the partitions of an integer as well as the dice rolls having a given sum and the probability (equiprobability) of obtaining this sum.

Go to my link on Maple Cloud :

https://maple.cloud/app/6553582941372416/decomposition+d%27un+entier+en+somme+et+partage+d%27une+liste+en+parties+-+app+pour+lanc%C3%A9s+de+d%C3%A9s?key=8533B69313624ADCBE181600B15A1F4F5DDBBFAE86694B479CE492515459A2C9

 

Hello everyone,

I am creating this post to begin a thread where I will share a series of worksheets on important topics in Complex Analysis, written as part of my notes for my classes. Complex_Analysis_Notes.pdf

The planned sections include:

  • Section 1: Infinite Series

  • Section 2: Power Series

  • Section 3: The Radius of Convergence of a Power Series

  • Section 4: The Riemann Zeta Function and the Riemann Hypothesis

  • Section 5: The Prime Number Theorem

Each worksheet will include calculations, plots, and examples using Maple to illustrate key ideas.

I plan to upload one worksheet every week to keep a steady pace and allow time for discussion and feedback between posts.

I hope this thread will be helpful both for learning and for deeper exploration.
Feel free to comment, suggest improvements, or ask questions as I post the materials.

Thank you!

restart; interface(imaginaryunit = 'I'); z := I*(1/3); S_N := proc (n) options operator, arrow; sum(z^k, k = 0 .. n) end proc; limit(S_N(n), n = infinity); S_N(10); S_N(100); S_N(1000); with(plots); points := [seq([Re(evalf(S_N(n))), Im(evalf(S_N(n)))], n = 0 .. 50)]; pointplot(points, connect = true, symbol = solidcircle, symbolsize = 10, color = blue, labels = ["Re", "Im"])

proc (n) options operator, arrow; sum(z^k, k = 0 .. n) end proc

 

9/10+(3/10)*I

 

53144/59049+(5905/19683)*I

 

 

restart; interface(imaginaryunit = 'I'); z := I*(1/2); S_N := proc (n) options operator, arrow; sum(z^k, k = 0 .. n) end proc; limit(S_N(n), n = infinity); S_N(10); S_N(100); S_N(1000); with(plots); points := [seq([Re(evalf(S_N(n))), Im(evalf(S_N(n)))], n = 0 .. 50)]; pointplot(points, connect = true, symbol = solidcircle, symbolsize = 10, color = blue, labels = ["Re", "Im"])

proc (n) options operator, arrow; sum(z^k, k = 0 .. n) end proc

 

4/5+(2/5)*I

 

819/1024+(205/512)*I

 

 

NULL

restart; with(plots); interface(imaginaryunit = 'I'); S := proc (N) local n; sum(((1/2)*I)^n, n = 0 .. N) end proc; fullsum := sum(((1/2)*I)^n, n = 0 .. infinity); realpts := [seq([n, Re(S(n))], n = 0 .. 30)]; imagpts := [seq([n, Im(S(n))], n = 0 .. 30)]; limit(Re(S(n)), n = infinity); limit(Im(S(n)), n = infinity); horiz_reallimit := plot(4/5, k = 0 .. 30, color = black, linestyle = 2, thickness = 2); horiz_imaglimit := plot(2/5, k = 0 .. 30, color = black, linestyle = 2, thickness = 2); plots[display]([pointplot(realpts, symbol = solidcircle, style = pointline, color = blue, labels = ["n", "Value"], legend = "Real Part"), pointplot(imagpts, symbol = solidbox, style = pointline, color = red, labels = ["n", "Value"], legend = "Imaginary Part"), horiz_reallimit, horiz_imaglimit], axes = boxed, labels = ["n", "Partial Sum Value"])

4/5+(2/5)*I

 

4/5

 

2/5

 

 

restart; with(plots); interface(imaginaryunit = 'I'); H := proc (N) local n; sum(1/n, n = 1 .. N) end proc; limit(H(n), n = infinity); limit(Re(H(n)), n = infinity); limit(Im(H(n)), n = infinity); harmonic_pts := [seq([n, H(n)], n = 1 .. 100)]; harmonic_plot := pointplot(harmonic_pts, symbol = solidcircle, style = pointline, color = blue, labels = ["n", "Partial Sum Value"], axes = boxed)

infinity

 

infinity

 

0

 

 

restart; with(plots); interface(imaginaryunit = 'I'); S := proc (N) local n; sum(I^k/k, k = 1 .. N) end proc; realpts := [seq([n, Re(S(n))], n = 1 .. 100)]; imagpts := [seq([n, Im(S(n))], n = 1 .. 100)]; complex_pts := [seq([Re(S(n)), Im(S(n))], n = 1 .. 100)]; S_infinite := sum(I^k/k, k = 1 .. infinity); Re(S_infinite); Im(S_infinite); horiz_reallimit := plot(-(1/2)*ln(2), k = 0 .. 100, color = black, linestyle = 2, thickness = 2); horiz_imaglimit := plot((1/4)*Pi, k = 0 .. 100, color = black, linestyle = 2, thickness = 2); real_plot := pointplot(realpts, symbol = solidcircle, style = pointline, color = blue, labels = ["n", "Partial Sum Value"], axes = boxed, legend = "Real Part"); imag_plot := pointplot(imagpts, symbol = solidbox, style = pointline, color = red, labels = ["n", "Partial Sum Value"], axes = boxed, legend = "Imaginary Part"); plots[display]([real_plot, horiz_reallimit, imag_plot, horiz_imaglimit]); plots[pointplot](complex_pts, symbol = solidcircle, style = pointline, color = blue, axes = boxed, labels = ["Re", "Im"])

-(1/2)*ln(2)+((1/4)*I)*Pi

 

-(1/2)*ln(2)

 

(1/4)*Pi

 

 

 

restart; with(plots); interface(imaginaryunit = 'I'); S := proc (N) local n; sum((-(2/3)*I)^n, n = 0 .. N) end proc; fullsum := sum((-2*I*(1/3))^n, n = 0 .. infinity); realpts := [seq([n, Re(S(n))], n = 0 .. 30)]; imagpts := [seq([n, Im(S(n))], n = 0 .. 30)]; limit(Re(S(n)), n = infinity); limit(Im(S(n)), n = infinity); horiz_reallimit := plot(9/13, k = 0 .. 30, color = black, linestyle = 2, thickness = 2); horiz_imaglimit := plot(-6/13, k = 0 .. 30, color = black, linestyle = 2, thickness = 2); plots[display]([pointplot(realpts, symbol = solidcircle, style = pointline, color = blue, labels = ["n", "Value"], legend = "Real Part"), pointplot(imagpts, symbol = solidbox, style = pointline, color = red, labels = ["n", "Value"], legend = "Imaginary Part"), horiz_reallimit, horiz_imaglimit], axes = boxed, labels = ["n", "Partial Sum Value"])

9/13-(6/13)*I

 

9/13

 

-6/13

 

 

restart; with(plots); interface(imaginaryunit = 'I'); S := proc (N) local n; sum((I*Pi)^n, n = 0 .. N) end proc; realpts := [seq([n, Re(S(n))], n = 0 .. 100)]; imagpts := [seq([n, Im(S(n))], n = 0 .. 100)]; complex_pts := [seq([Re(S(n)), Im(S(n))], n = 0 .. 100)]; limit(S(N), N = infinity); limit(Re(S(n)), n = infinity); limit(Im(S(n)), n = infinity); real_plot := pointplot(realpts, symbol = solidcircle, style = pointline, color = blue, labels = ["n", "Partial Sum (Real Part)"], title = "Real Part of Partial Sums of (Pi i)^n", axes = boxed); imag_plot := pointplot(imagpts, symbol = solidbox, style = pointline, color = red, labels = ["n", "Partial Sum (Imaginary Part)"], title = "Imaginary Part of Partial Sums of (Pi i)^n", axes = boxed); complex_plot := pointplot(complex_pts, symbol = solidcircle, style = pointline, color = blue, labels = ["Re", "Im"], title = "Partial Sums in Complex Plane (Pi i)^n", axes = boxed)

undefined

 

undefined

 

undefined

 

 

 

 

restart; with(plots); interface(imaginaryunit = 'I'); S := proc (N) local n; sum(2*I^k/k, k = 1 .. N) end proc; realpts := [seq([n, Re(S(n))], n = 1 .. 100)]; imagpts := [seq([n, Im(S(n))], n = 1 .. 100)]; complex_pts := [seq([Re(S(n)), Im(S(n))], n = 1 .. 100)]; S_infinite := sum(2*I^k/k, k = 1 .. infinity); Re(S_infinite); Im(S_infinite); horiz_reallimit := plot(-ln(2), k = 0 .. 100, color = black, linestyle = 2, thickness = 2); horiz_imaglimit := plot((1/2)*Pi, k = 0 .. 100, color = black, linestyle = 2, thickness = 2); real_plot := pointplot(realpts, symbol = solidcircle, style = pointline, color = blue, labels = ["n", "Partial Sum Value"], axes = boxed, legend = "Real Part"); imag_plot := pointplot(imagpts, symbol = solidbox, style = pointline, color = red, labels = ["n", "Partial Sum Value"], axes = boxed, legend = "Imaginary Part"); plots[display]([real_plot, horiz_reallimit, imag_plot, horiz_imaglimit]); plots[pointplot](complex_pts, symbol = solidcircle, style = pointline, color = blue, axes = boxed, labels = ["Re", "Im"])

-ln(2)+((1/2)*I)*Pi

 

-ln(2)

 

(1/2)*Pi

 

 

 

restart; with(plots); interface(imaginaryunit = 'I'); S := proc (N) local n; add(exp(Pi*I*n)/n, n = 1 .. N) end proc; realpts := [seq([n, Re(S(n))], n = 1 .. 100)]; imagpts := [seq([n, Im(S(n))], n = 1 .. 100)]; complex_pts := [seq([Re(S(n)), Im(S(n))], n = 1 .. 100)]; S_infinite := sum(exp(Pi*I*n)/n, n = 1 .. infinity); limit_Re := Re(S_infinite); limit_Im := Im(S_infinite); limit_Re; limit_Im; real_plot := pointplot(realpts, symbol = solidcircle, style = pointline, color = blue, labels = ["n", "Partial Sum Value"], title = "Real Part of Partial Sums", axes = boxed); imag_plot := pointplot(imagpts, symbol = solidbox, style = pointline, color = red, labels = ["n", "Partial Sum Value"], title = "Imaginary Part of Partial Sums", axes = boxed); complex_plot := pointplot(complex_pts, symbol = solidcircle, style = pointline, color = blue, labels = ["Re", "Im"], title = "Partial Sums in Complex Plane", axes = boxed); plots[display]([real_plot, imag_plot]); plots[display](complex_plot)

-ln(2)

 

-ln(2)

 

0

 

-ln(2)

 

0

 

 

 

restart; with(plots); interface(imaginaryunit = 'I'); S := proc (N) local n; add(exp(2*Pi*I*n), n = 0 .. N) end proc; realpts := [seq([n, Re(S(n))], n = 0 .. 100)]; imagpts := [seq([n, Im(S(n))], n = 0 .. 100)]; complex_pts := [seq([Re(S(n)), Im(S(n))], n = 0 .. 100)]; S_infinite := sum(exp(2*Pi*I*n), n = 1 .. infinity); limit_Re := Re(S_infinite); limit_Im := Im(S_infinite); real_plot := pointplot(realpts, symbol = solidcircle, style = pointline, color = blue, labels = ["n", "Partial Sum Value"], title = "Real Part of Partial Sums", axes = boxed); imag_plot := pointplot(imagpts, symbol = solidbox, style = pointline, color = red, labels = ["n", "Partial Sum Value"], title = "Imaginary Part of Partial Sums", axes = boxed); complex_plot := pointplot(complex_pts, symbol = solidcircle, style = pointline, color = blue, labels = ["Re", "Im"], title = "Partial Sums in Complex Plane", axes = boxed); plots[display]([real_plot, imag_plot]); plots[display](complex_plot)

infinity

 

infinity

 

0

 

 

 
 

``

Download infinite_series.mw

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 1st, 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 A[mu](x, y, z, t), there are then four coordinates, one for each of the values of the index mu, and there are four parameters x, y, z, t.

 

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

 

with(Physics)

Setup(mathematicalnotation = true, coordinates = cartesian)

[coordinatesystems = {X}, mathematicalnotation = true]

(1)

The lambda*Phi^4 model in classical field theory and corresponding field equations, as in previous releases

CompactDisplay(Phi(X))

Phi(x, y, z, t)*`will now be displayed as`*Phi

(2)

L := (1/2)*d_[mu](Phi(X))*d_[mu](Phi(X))-(1/2)*m^2*Phi(X)^2+(1/4)*lambda*Phi(X)^4

(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

LagrangeEquations(L, Phi)

Phi(X)^3*lambda-Phi(X)*m^2-Physics:-dAlembertian(Phi(X), [X]) = 0

(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 2/sqrt(-%g_[determinant]) (in flat spacetimes sqrt(-%g_[determinant]) = 1). 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

LagrangeEquations(L, EnergyMomentum[mu, nu])

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 T[mu, nu], you can use the Define  command

Define(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]))

`Defined objects with tensor properties`

 

{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 T[mu, nu]

EnergyMomentum[definition]

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)

Library:-IsTensorialSymmetric(EnergyMomentum[mu, nu])

true

(8)

EnergyMomentum[]

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

g_[arb]

_______________________________________________________

 

`Setting `*lowercaselatin_is*` letters to represent `*space*` indices`

 

`The arbitrary metric in coordinates `*[x, y, z, t]

 

`Signature: `(`- - - +`)

 

_______________________________________________________

 

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

L := sqrt(-%g_[determinant])*L

(-%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

LagrangeEquations(L, Phi)

Phi(X)^3*lambda-Phi(X)*m^2-Physics:-D_[kappa](Physics:-d_[`~kappa`](Phi(X), [X]), [X]) = 0

(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

LagrangeEquations(L, EnergyMomentum[mu, nu])

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 g[mu, nu] 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 g[mu, nu], and as in the case of electrodynamics the parameters are the spacetime coordinates X^alpha. 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

L := sqrt(-%g_[determinant])*Ricci[alpha, `~alpha`]

(-%g_[determinant])^(1/2)*Physics:-Ricci[alpha, `~alpha`]

(14)

Einstein's equations in vacuum:

LagrangeEquations(L, g_[mu, nu])

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

(15)

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

Einstein[definition]

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

(16)

The Lagrangian L 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

L__C := convert(L, Christoffel)

(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*(Physics:-d_[nu](Physics:-Christoffel[`~nu`, alpha, lambda], [X])-Physics:-d_[lambda](Physics:-Christoffel[`~nu`, alpha, nu], [X])+Physics:-Christoffel[`~beta`, alpha, lambda]*Physics:-Christoffel[`~nu`, beta, nu]-Physics:-Christoffel[`~beta`, alpha, nu]*Physics:-Christoffel[`~nu`, beta, lambda])

(17)

Recalling the definition

Christoffel[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 L[C] the two terms containing derivatives of Christoffel symbols contain second order derivatives of g[mu, nu]. Now, it is always possible to add a total spacetime derivative to L[C] 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 L[C], obtain a Lagrangian involving only 1st order derivatives. The total derivative, expressed using the inert `&PartialD;` command to see it before the differentiation operation is performed, is

TD := %d_[alpha](g_[`~mu`, `~nu`]*sqrt(-%g_[determinant])*(g_[`~alpha`, mu]*Christoffel[`~beta`, nu, beta]-Christoffel[`~alpha`, mu, nu]))

%d_[alpha](Physics:-g_[`~mu`, `~nu`]*(-%g_[determinant])^(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-Christoffel[`~beta`, beta, nu]-Physics:-Christoffel[`~alpha`, mu, nu]))

(19)

Adding this term to L[C], performing the `&PartialD;` differentiation operation and simplifying we get

L__1 := L__C+TD

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

(20)

L__1 := eval(L__1, %d_ = d_)

(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*(Physics:-d_[nu](Physics:-Christoffel[`~nu`, alpha, lambda], [X])-Physics:-d_[lambda](Physics:-Christoffel[`~nu`, alpha, nu], [X])+Physics:-Christoffel[`~beta`, alpha, lambda]*Physics:-Christoffel[`~nu`, beta, nu]-Physics:-Christoffel[`~beta`, alpha, nu]*Physics:-Christoffel[`~nu`, beta, lambda])+Physics:-d_[alpha](Physics:-g_[`~mu`, `~nu`], [X])*(-%g_[determinant])^(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-Christoffel[`~beta`, beta, nu]-Physics:-Christoffel[`~alpha`, mu, nu])-(1/2)*Physics:-g_[`~mu`, `~nu`]*(Physics:-g_[mu, `~alpha`]*Physics:-Christoffel[`~beta`, beta, nu]-Physics:-Christoffel[`~alpha`, mu, nu])*%g_[determinant]*Physics:-g_[`~kappa`, `~lambda`]*Physics:-d_[alpha](Physics:-g_[kappa, lambda], [X])/(-%g_[determinant])^(1/2)+Physics:-g_[`~mu`, `~nu`]*(-%g_[determinant])^(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-d_[alpha](Physics:-Christoffel[`~beta`, beta, nu], [X])-Physics:-d_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X]))

(21)

L__1 := Simplify(L__1)

(Physics:-Christoffel[alpha, beta, kappa]*Physics:-Christoffel[`~beta`, `~alpha`, `~kappa`]-Physics:-Christoffel[alpha, beta, `~alpha`]*Physics:-Christoffel[`~beta`, kappa, `~kappa`])*(-%g_[determinant])^(1/2)

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

LagrangeEquations(L__1, g_[mu, nu])

-(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 `&equiv;`(L[1], (Physics[Christoffel][alpha, beta, kappa]*Physics[Christoffel][`~beta`, `~alpha`, `~kappa`]-Physics[Christoffel][alpha, beta, `~alpha`]*Physics[Christoffel][`~beta`, kappa, `~kappa`])*(-%g_[determinant])^(1/2)) is no just L[C] ≡ (17) after discarding its two terms involving derivatives of Christoffel symbols. To verify this, split L[C] into the terms containing or not derivatives of Christoffel

L__22, L__11 := selectremove(has, expand(L__C), d_)

(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*Physics:-d_[nu](Physics:-Christoffel[`~nu`, alpha, lambda], [X])-(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*Physics:-d_[lambda](Physics:-Christoffel[`~nu`, alpha, nu], [X]), (-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*Physics:-Christoffel[`~beta`, alpha, lambda]*Physics:-Christoffel[`~nu`, beta, nu]-(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*Physics:-Christoffel[`~beta`, alpha, nu]*Physics:-Christoffel[`~nu`, beta, lambda]

(24)

Comparing, the total derivative TD≡ (19) is not just -L[22], but

TD = -L__22-2*L__11

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

(25)

Things like these, TD = -L__22-2*L__11, 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 `&PartialD;` and simplify to see the equality is true

(lhs-rhs)(%d_[alpha](Physics[g_][`~mu`, `~nu`]*(-%g_[determinant])^(1/2)*(Physics[g_][mu, `~alpha`]*Physics[Christoffel][`~beta`, beta, nu]-Physics[Christoffel][`~alpha`, mu, nu])) = -(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[d_][nu](Physics[Christoffel][`~nu`, alpha, lambda], [X])+(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[d_][lambda](Physics[Christoffel][`~nu`, alpha, nu], [X])-2*(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[Christoffel][`~beta`, alpha, lambda]*Physics[Christoffel][`~nu`, beta, nu]+2*(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[Christoffel][`~beta`, alpha, nu]*Physics[Christoffel][`~nu`, beta, lambda])

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

(26)

eval(%d_[alpha](Physics[g_][`~mu`, `~nu`]*(-%g_[determinant])^(1/2)*(Physics[g_][mu, `~alpha`]*Physics[Christoffel][`~beta`, beta, nu]-Physics[Christoffel][`~alpha`, mu, nu]))+(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[d_][nu](Physics[Christoffel][`~nu`, alpha, lambda], [X])-(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[d_][lambda](Physics[Christoffel][`~nu`, alpha, nu], [X])+2*(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[Christoffel][`~beta`, alpha, lambda]*Physics[Christoffel][`~nu`, beta, nu]-2*(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[Christoffel][`~beta`, alpha, nu]*Physics[Christoffel][`~nu`, beta, lambda], %d_ = d_)

(-%g_[determinant])^(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-Christoffel[`~beta`, beta, nu]-Physics:-Christoffel[`~alpha`, mu, nu])*Physics:-d_[alpha](Physics:-g_[`~mu`, `~nu`], [X])-(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-Christoffel[`~beta`, beta, nu]-Physics:-Christoffel[`~alpha`, mu, nu])*Physics:-g_[`~mu`, `~nu`]*%g_[determinant]*Physics:-g_[`~kappa`, `~lambda`]*Physics:-d_[alpha](Physics:-g_[kappa, lambda], [X])/(-%g_[determinant])^(1/2)+(-%g_[determinant])^(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-d_[alpha](Physics:-Christoffel[`~beta`, beta, nu], [X])-Physics:-d_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X]))*Physics:-g_[`~mu`, `~nu`]+(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*Physics:-d_[nu](Physics:-Christoffel[`~nu`, alpha, lambda], [X])-(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*Physics:-d_[lambda](Physics:-Christoffel[`~nu`, alpha, nu], [X])+2*(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*Physics:-Christoffel[`~beta`, alpha, lambda]*Physics:-Christoffel[`~nu`, beta, nu]-2*(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*Physics:-Christoffel[`~beta`, alpha, nu]*Physics:-Christoffel[`~nu`, beta, lambda]

(27)

Simplify((-%g_[determinant])^(1/2)*(Physics[g_][mu, `~alpha`]*Physics[Christoffel][`~beta`, beta, nu]-Physics[Christoffel][`~alpha`, mu, nu])*Physics[d_][alpha](Physics[g_][`~mu`, `~nu`], [X])-(1/2)*(Physics[g_][mu, `~alpha`]*Physics[Christoffel][`~beta`, beta, nu]-Physics[Christoffel][`~alpha`, mu, nu])*Physics[g_][`~mu`, `~nu`]*%g_[determinant]*Physics[g_][`~kappa`, `~lambda`]*Physics[d_][alpha](Physics[g_][kappa, lambda], [X])/(-%g_[determinant])^(1/2)+(-%g_[determinant])^(1/2)*(Physics[g_][mu, `~alpha`]*Physics[d_][alpha](Physics[Christoffel][`~beta`, beta, nu], [X])-Physics[d_][alpha](Physics[Christoffel][`~alpha`, mu, nu], [X]))*Physics[g_][`~mu`, `~nu`]+(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[d_][nu](Physics[Christoffel][`~nu`, alpha, lambda], [X])-(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[d_][lambda](Physics[Christoffel][`~nu`, alpha, nu], [X])+2*(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[Christoffel][`~beta`, alpha, lambda]*Physics[Christoffel][`~nu`, beta, nu]-2*(-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*Physics[Christoffel][`~beta`, alpha, nu]*Physics[Christoffel][`~nu`, beta, lambda])

0

(28)

That said, it is also true that TD = -L[22]-2*L[11] results in the Lagrangian L[1] = -L[11], and since the equations of movement don't depend on the sign of the Lagrangian, for this Lagrangian `&equiv;`(L[C], (-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*(Physics[d_][nu](Physics[Christoffel][`~nu`, alpha, lambda], [X])-Physics[d_][lambda](Physics[Christoffel][`~nu`, alpha, nu], [X])+Physics[Christoffel][`~beta`, alpha, lambda]*Physics[Christoffel][`~nu`, beta, nu]-Physics[Christoffel][`~beta`, alpha, nu]*Physics[Christoffel][`~nu`, beta, lambda])) adding the term TD happens to be equivalent to just discarding the terms of L__C 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,

S := Intc(L, X)

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 ":"

EE__unsimplified := Fundiff(S, g_[alpha, beta]) = 0

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

EEC := Simplify(EE__unsimplified)

(1/4)*(2*Physics:-Christoffel[chi, iota, kappa]*Physics:-Christoffel[`~iota`, `~chi`, `~kappa`]-2*Physics:-Christoffel[chi, iota, `~chi`]*Physics:-Christoffel[`~iota`, kappa, `~kappa`]-2*Physics:-D_[iota](Physics:-Christoffel[chi, `~chi`, `~iota`], [X])+2*Physics:-D_[chi](Physics:-Christoffel[`~chi`, iota, `~iota`], [X]))*Physics:-g_[`~alpha`, `~beta`]+(1/4)*(Physics:-Christoffel[`~alpha`, chi, `~beta`]+Physics:-Christoffel[`~beta`, chi, `~alpha`])*Physics:-Christoffel[`~chi`, iota, `~iota`]-(1/4)*Physics:-Christoffel[`~beta`, chi, iota]*Physics:-Christoffel[`~chi`, `~alpha`, `~iota`]-(1/2)*Physics:-Christoffel[chi, iota, `~alpha`]*Physics:-Christoffel[`~iota`, `~beta`, `~chi`]+(1/2)*Physics:-Christoffel[chi, `~alpha`, `~beta`]*Physics:-Christoffel[iota, `~chi`, `~iota`]-(1/4)*Physics:-Christoffel[`~alpha`, chi, iota]*Physics:-Christoffel[`~chi`, `~beta`, `~iota`]+(1/4)*Physics:-D_[chi](Physics:-Christoffel[`~alpha`, `~beta`, `~chi`], [X])+(1/4)*Physics:-D_[chi](Physics:-Christoffel[`~beta`, `~alpha`, `~chi`], [X])-(1/2)*Physics:-D_[chi](Physics:-Christoffel[`~chi`, `~alpha`, `~beta`], [X])-(1/4)*Physics:-D_[`~alpha`](Physics:-Christoffel[`~beta`, chi, `~chi`], [X])+(1/2)*Physics:-D_[`~beta`](Physics:-Christoffel[chi, `~alpha`, `~chi`], [X])-(1/4)*Physics:-D_[`~beta`](Physics:-Christoffel[`~alpha`, chi, `~chi`], [X]) = 0

(30)

In this result, we see`&dtri;` 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,

"`&dtri;`[chi](GAMMA[]^(alpha,beta,chi)) :"

% = expand(%)

Physics:-D_[chi](Physics:-Christoffel[`~alpha`, `~beta`, `~chi`], [X]) = Physics:-d_[chi](Physics:-Christoffel[`~alpha`, `~beta`, `~chi`], [X])+Physics:-Christoffel[`~alpha`, chi, mu]*Physics:-Christoffel[`~mu`, `~beta`, `~chi`]+Physics:-Christoffel[`~beta`, chi, mu]*Physics:-Christoffel[`~alpha`, `~chi`, `~mu`]+Physics:-Christoffel[`~chi`, chi, mu]*Physics:-Christoffel[`~alpha`, `~beta`, `~mu`]

(31)

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

Ricci[definition]

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

(32)

Rewrite the noncovariant derivatives `&PartialD;` in terms of`&dtri;` derivatives using the computational representation (31), simplify and isolate one of them

convert(Physics[Ricci][mu, nu] = Physics[d_][alpha](Physics[Christoffel][`~alpha`, mu, nu], [X])-Physics[d_][nu](Physics[Christoffel][`~alpha`, mu, alpha], [X])+Physics[Christoffel][`~beta`, mu, nu]*Physics[Christoffel][`~alpha`, beta, alpha]-Physics[Christoffel][`~beta`, mu, alpha]*Physics[Christoffel][`~alpha`, nu, beta], D_)

Physics:-Ricci[mu, nu] = Physics:-D_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X])-Physics:-Christoffel[`~alpha`, alpha, kappa]*Physics:-Christoffel[`~kappa`, mu, nu]+Physics:-Christoffel[`~kappa`, alpha, mu]*Physics:-Christoffel[`~alpha`, kappa, nu]+Physics:-Christoffel[`~kappa`, alpha, nu]*Physics:-Christoffel[`~alpha`, mu, kappa]-Physics:-D_[nu](Physics:-Christoffel[`~alpha`, alpha, mu], [X])-Physics:-Christoffel[`~lambda`, mu, nu]*Physics:-Christoffel[`~alpha`, alpha, lambda]+Physics:-Christoffel[`~beta`, mu, nu]*Physics:-Christoffel[`~alpha`, alpha, beta]-Physics:-Christoffel[`~beta`, alpha, mu]*Physics:-Christoffel[`~alpha`, beta, nu]

(33)

Simplify(Physics[Ricci][mu, nu] = D_[alpha](Physics[Christoffel][`~alpha`, mu, nu], [X])-Physics[Christoffel][`~alpha`, alpha, kappa]*Physics[Christoffel][`~kappa`, mu, nu]+Physics[Christoffel][`~kappa`, alpha, mu]*Physics[Christoffel][`~alpha`, kappa, nu]+Physics[Christoffel][`~kappa`, alpha, nu]*Physics[Christoffel][`~alpha`, mu, kappa]-D_[nu](Physics[Christoffel][`~alpha`, alpha, mu], [X])-Physics[Christoffel][`~lambda`, mu, nu]*Physics[Christoffel][`~alpha`, alpha, lambda]+Physics[Christoffel][`~beta`, mu, nu]*Physics[Christoffel][`~alpha`, alpha, beta]-Physics[Christoffel][`~beta`, alpha, mu]*Physics[Christoffel][`~alpha`, beta, nu])

Physics:-Ricci[mu, nu] = Physics:-Christoffel[alpha, beta, mu]*Physics:-Christoffel[`~beta`, nu, `~alpha`]-Physics:-Christoffel[beta, mu, nu]*Physics:-Christoffel[alpha, `~alpha`, `~beta`]+Physics:-D_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X])-Physics:-D_[nu](Physics:-Christoffel[alpha, mu, `~alpha`], [X])

(34)

C_to_Ricci := isolate(Physics[Ricci][mu, nu] = Physics[Christoffel][alpha, beta, mu]*Physics[Christoffel][`~beta`, nu, `~alpha`]-Physics[Christoffel][beta, mu, nu]*Physics[Christoffel][alpha, `~alpha`, `~beta`]+D_[alpha](Physics[Christoffel][`~alpha`, mu, nu], [X])-D_[nu](Physics[Christoffel][alpha, mu, `~alpha`], [X]), D_[alpha](Christoffel[`~alpha`, mu, nu]))

Physics:-D_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X]) = -Physics:-Christoffel[alpha, beta, mu]*Physics:-Christoffel[`~beta`, nu, `~alpha`]+Physics:-Christoffel[beta, mu, nu]*Physics:-Christoffel[alpha, `~alpha`, `~beta`]+Physics:-Ricci[mu, nu]+Physics:-D_[nu](Physics:-Christoffel[alpha, mu, `~alpha`], [X])

(35)

Analogously, derive an expression to rewrite`&dtri;` derivatives of Christoffel symbols using the Riemann  tensor

Riemann[`~alpha`, beta, mu, nu, definition]

Physics:-Riemann[`~alpha`, beta, mu, nu] = Physics:-d_[mu](Physics:-Christoffel[`~alpha`, beta, nu], [X])-Physics:-d_[nu](Physics:-Christoffel[`~alpha`, beta, mu], [X])+Physics:-Christoffel[`~alpha`, upsilon, mu]*Physics:-Christoffel[`~upsilon`, beta, nu]-Physics:-Christoffel[`~alpha`, upsilon, nu]*Physics:-Christoffel[`~upsilon`, beta, mu]

(36)

convert(Physics[Riemann][`~alpha`, beta, mu, nu] = Physics[d_][mu](Physics[Christoffel][`~alpha`, beta, nu], [X])-Physics[d_][nu](Physics[Christoffel][`~alpha`, beta, mu], [X])+Physics[Christoffel][`~alpha`, upsilon, mu]*Physics[Christoffel][`~upsilon`, beta, nu]-Physics[Christoffel][`~alpha`, upsilon, nu]*Physics[Christoffel][`~upsilon`, beta, mu], D_)

Physics:-Riemann[`~alpha`, beta, mu, nu] = Physics:-D_[mu](Physics:-Christoffel[`~alpha`, beta, nu], [X])+Physics:-Christoffel[`~kappa`, mu, nu]*Physics:-Christoffel[`~alpha`, beta, kappa]-Physics:-Christoffel[`~alpha`, kappa, mu]*Physics:-Christoffel[`~kappa`, beta, nu]+Physics:-Christoffel[`~kappa`, beta, mu]*Physics:-Christoffel[`~alpha`, kappa, nu]-Physics:-D_[nu](Physics:-Christoffel[`~alpha`, beta, mu], [X])-Physics:-Christoffel[`~lambda`, mu, nu]*Physics:-Christoffel[`~alpha`, beta, lambda]-Physics:-Christoffel[`~lambda`, beta, nu]*Physics:-Christoffel[`~alpha`, lambda, mu]+Physics:-Christoffel[`~alpha`, lambda, nu]*Physics:-Christoffel[`~lambda`, beta, mu]+Physics:-Christoffel[`~alpha`, mu, upsilon]*Physics:-Christoffel[`~upsilon`, beta, nu]-Physics:-Christoffel[`~alpha`, nu, upsilon]*Physics:-Christoffel[`~upsilon`, beta, mu]

(37)

Simplify(Physics[Riemann][`~alpha`, beta, mu, nu] = D_[mu](Physics[Christoffel][`~alpha`, beta, nu], [X])+Physics[Christoffel][`~kappa`, mu, nu]*Physics[Christoffel][`~alpha`, beta, kappa]-Physics[Christoffel][`~alpha`, kappa, mu]*Physics[Christoffel][`~kappa`, beta, nu]+Physics[Christoffel][`~kappa`, beta, mu]*Physics[Christoffel][`~alpha`, kappa, nu]-D_[nu](Physics[Christoffel][`~alpha`, beta, mu], [X])-Physics[Christoffel][`~lambda`, mu, nu]*Physics[Christoffel][`~alpha`, beta, lambda]-Physics[Christoffel][`~lambda`, beta, nu]*Physics[Christoffel][`~alpha`, lambda, mu]+Physics[Christoffel][`~alpha`, lambda, nu]*Physics[Christoffel][`~lambda`, beta, mu]+Physics[Christoffel][`~alpha`, mu, upsilon]*Physics[Christoffel][`~upsilon`, beta, nu]-Physics[Christoffel][`~alpha`, nu, upsilon]*Physics[Christoffel][`~upsilon`, beta, mu])

Physics:-Riemann[`~alpha`, beta, mu, nu] = -Physics:-Christoffel[`~alpha`, kappa, mu]*Physics:-Christoffel[`~kappa`, beta, nu]+Physics:-Christoffel[`~kappa`, beta, mu]*Physics:-Christoffel[`~alpha`, kappa, nu]+Physics:-D_[mu](Physics:-Christoffel[`~alpha`, beta, nu], [X])-Physics:-D_[nu](Physics:-Christoffel[`~alpha`, beta, mu], [X])

(38)

C_to_Riemann := isolate(Physics[Riemann][`~alpha`, beta, mu, nu] = -Physics[Christoffel][`~alpha`, kappa, mu]*Physics[Christoffel][`~kappa`, beta, nu]+Physics[Christoffel][`~kappa`, beta, mu]*Physics[Christoffel][`~alpha`, kappa, nu]+D_[mu](Physics[Christoffel][`~alpha`, beta, nu], [X])-D_[nu](Physics[Christoffel][`~alpha`, beta, mu], [X]), D_[mu](Christoffel[`~alpha`, beta, nu]))

Physics:-D_[mu](Physics:-Christoffel[`~alpha`, beta, nu], [X]) = Physics:-Christoffel[`~alpha`, kappa, mu]*Physics:-Christoffel[`~kappa`, beta, nu]-Physics:-Christoffel[`~kappa`, beta, mu]*Physics:-Christoffel[`~alpha`, kappa, nu]+Physics:-Riemann[`~alpha`, beta, mu, nu]+Physics:-D_[nu](Physics:-Christoffel[`~alpha`, beta, mu], [X])

(39)

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

Substitute(C_to_Riemann, C_to_Ricci, EEC)

-(1/4)*Physics:-Christoffel[`~beta`, psi, `~alpha`]*Physics:-Christoffel[`~psi`, omega, `~omega`]+(1/4)*Physics:-Christoffel[`~beta`, psi, `~omega`]*Physics:-Christoffel[`~psi`, omega, `~alpha`]+(1/4)*Physics:-Christoffel[`~alpha`, lambda, mu]*Physics:-Christoffel[`~lambda`, `~beta`, `~mu`]-(1/4)*Physics:-Christoffel[`~alpha`, lambda, `~mu`]*Physics:-Christoffel[`~lambda`, mu, `~beta`]+(1/4)*Physics:-Christoffel[`~beta`, nu, sigma]*Physics:-Christoffel[`~sigma`, `~alpha`, `~nu`]-(1/4)*Physics:-Christoffel[`~beta`, sigma, `~nu`]*Physics:-Christoffel[`~sigma`, nu, `~alpha`]+(1/2)*Physics:-Christoffel[omicron, zeta, `~beta`]*Physics:-Christoffel[`~zeta`, `~alpha`, `~omicron`]-(1/2)*Physics:-Christoffel[omicron, zeta, `~omicron`]*Physics:-Christoffel[`~zeta`, `~alpha`, `~beta`]-Physics:-Ricci[`~alpha`, `~beta`]-(1/2)*Physics:-Christoffel[`~tau`, `~alpha`, `~beta`]*Physics:-Christoffel[`~upsilon`, tau, upsilon]-(1/2)*Physics:-Christoffel[chi, iota, `~alpha`]*Physics:-Christoffel[`~iota`, `~beta`, `~chi`]+(1/4)*(Physics:-Christoffel[`~alpha`, chi, `~beta`]+Physics:-Christoffel[`~beta`, chi, `~alpha`])*Physics:-Christoffel[`~chi`, iota, `~iota`]+(1/2)*Physics:-Christoffel[`~upsilon`, tau, `~beta`]*Physics:-Christoffel[`~tau`, upsilon, `~alpha`]-(1/4)*Physics:-Christoffel[`~beta`, chi, iota]*Physics:-Christoffel[`~chi`, `~alpha`, `~iota`]+(1/2)*Physics:-Christoffel[chi, `~alpha`, `~beta`]*Physics:-Christoffel[iota, `~chi`, `~iota`]-(1/4)*Physics:-Christoffel[`~alpha`, chi, iota]*Physics:-Christoffel[`~chi`, `~beta`, `~iota`]-(1/4)*Physics:-D_[`~rho1`](Physics:-Christoffel[`~alpha`, rho1, `~beta`], [X])+(1/4)*Physics:-D_[`~mu`](Physics:-Christoffel[`~alpha`, mu, `~beta`], [X])+(1/4)*Physics:-D_[`~nu`](Physics:-Christoffel[`~beta`, nu, `~alpha`], [X])-(1/2)*Physics:-D_[`~beta`](Physics:-Christoffel[`~upsilon`, upsilon, `~alpha`], [X])-(1/4)*Physics:-D_[`~omega`](Physics:-Christoffel[`~beta`, omega, `~alpha`], [X])+(1/2)*Physics:-D_[`~beta`](Physics:-Christoffel[alpha9, `~alpha`, `~alpha9`], [X])-(1/4)*Physics:-Christoffel[`~alpha`, rho, `~beta`]*Physics:-Christoffel[`~rho`, rho1, `~rho1`]+(1/4)*Physics:-Christoffel[`~alpha`, rho, `~rho1`]*Physics:-Christoffel[`~rho`, rho1, `~beta`]+(1/2)*Physics:-Christoffel[alpha10, `~alpha`, `~beta`]*Physics:-Christoffel[alpha9, `~alpha10`, `~alpha9`]-(1/2)*Physics:-Christoffel[alpha9, alpha10, `~alpha`]*Physics:-Christoffel[`~alpha10`, `~alpha9`, `~beta`]+(1/4)*(2*Physics:-Christoffel[chi, iota, kappa]*Physics:-Christoffel[`~iota`, `~chi`, `~kappa`]-2*Physics:-Christoffel[chi, iota, `~chi`]*Physics:-Christoffel[`~iota`, kappa, `~kappa`]-2*Physics:-Christoffel[alpha4, alpha5, alpha6]*Physics:-Christoffel[`~alpha6`, `~alpha4`, `~alpha5`]+2*Physics:-Christoffel[alpha4, alpha6, `~alpha5`]*Physics:-Christoffel[`~alpha6`, alpha5, `~alpha4`]-2*Physics:-D_[`~alpha5`](Physics:-Christoffel[alpha4, alpha5, `~alpha4`], [X])+2*Physics:-Christoffel[`~alpha1`, alpha1, alpha3]*Physics:-Christoffel[`~alpha3`, alpha2, `~alpha2`]-2*Physics:-Christoffel[`~alpha1`, alpha3, `~alpha2`]*Physics:-Christoffel[`~alpha3`, alpha1, alpha2]+2*Physics:-Ricci[alpha2, `~alpha2`]+2*Physics:-D_[`~alpha2`](Physics:-Christoffel[`~alpha1`, alpha1, alpha2], [X]))*Physics:-g_[`~alpha`, `~beta`] = 0

(40)

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

Simplify(-(1/4)*Physics[Christoffel][`~beta`, psi, `~alpha`]*Physics[Christoffel][`~psi`, omega, `~omega`]+(1/4)*Physics[Christoffel][`~beta`, psi, `~omega`]*Physics[Christoffel][`~psi`, omega, `~alpha`]+(1/4)*Physics[Christoffel][`~alpha`, lambda, mu]*Physics[Christoffel][`~lambda`, `~beta`, `~mu`]-(1/4)*Physics[Christoffel][`~alpha`, lambda, `~mu`]*Physics[Christoffel][`~lambda`, mu, `~beta`]+(1/4)*Physics[Christoffel][`~beta`, nu, sigma]*Physics[Christoffel][`~sigma`, `~alpha`, `~nu`]-(1/4)*Physics[Christoffel][`~beta`, sigma, `~nu`]*Physics[Christoffel][`~sigma`, nu, `~alpha`]+(1/2)*Physics[Christoffel][omicron, zeta, `~beta`]*Physics[Christoffel][`~zeta`, `~alpha`, `~omicron`]-(1/2)*Physics[Christoffel][omicron, zeta, `~omicron`]*Physics[Christoffel][`~zeta`, `~alpha`, `~beta`]-(1/2)*Physics[Christoffel][`~tau`, `~alpha`, `~beta`]*Physics[Christoffel][`~upsilon`, tau, upsilon]-(1/2)*Physics[Christoffel][chi, iota, `~alpha`]*Physics[Christoffel][`~iota`, `~beta`, `~chi`]+(1/4)*(Physics[Christoffel][`~alpha`, chi, `~beta`]+Physics[Christoffel][`~beta`, chi, `~alpha`])*Physics[Christoffel][`~chi`, iota, `~iota`]+(1/2)*Physics[Christoffel][`~upsilon`, tau, `~beta`]*Physics[Christoffel][`~tau`, upsilon, `~alpha`]-(1/4)*Physics[Christoffel][`~beta`, chi, iota]*Physics[Christoffel][`~chi`, `~alpha`, `~iota`]+(1/2)*Physics[Christoffel][chi, `~alpha`, `~beta`]*Physics[Christoffel][iota, `~chi`, `~iota`]-(1/4)*Physics[Christoffel][`~alpha`, chi, iota]*Physics[Christoffel][`~chi`, `~beta`, `~iota`]-(1/4)*Physics[Christoffel][`~alpha`, rho, `~beta`]*Physics[Christoffel][`~rho`, rho1, `~rho1`]+(1/4)*Physics[Christoffel][`~alpha`, rho, `~rho1`]*Physics[Christoffel][`~rho`, rho1, `~beta`]+(1/2)*Physics[Christoffel][alpha10, `~alpha`, `~beta`]*Physics[Christoffel][alpha9, `~alpha10`, `~alpha9`]-(1/2)*Physics[Christoffel][alpha9, alpha10, `~alpha`]*Physics[Christoffel][`~alpha10`, `~alpha9`, `~beta`]+(1/4)*(2*Physics[Christoffel][chi, iota, kappa]*Physics[Christoffel][`~iota`, `~chi`, `~kappa`]-2*Physics[Christoffel][chi, iota, `~chi`]*Physics[Christoffel][`~iota`, kappa, `~kappa`]-2*Physics[Christoffel][alpha4, alpha5, alpha6]*Physics[Christoffel][`~alpha6`, `~alpha4`, `~alpha5`]+2*Physics[Christoffel][alpha4, alpha6, `~alpha5`]*Physics[Christoffel][`~alpha6`, alpha5, `~alpha4`]-2*D_[`~alpha5`](Physics[Christoffel][alpha4, alpha5, `~alpha4`], [X])+2*Physics[Christoffel][`~alpha1`, alpha1, alpha3]*Physics[Christoffel][`~alpha3`, alpha2, `~alpha2`]-2*Physics[Christoffel][`~alpha1`, alpha3, `~alpha2`]*Physics[Christoffel][`~alpha3`, alpha1, alpha2]+2*Physics[Ricci][alpha2, `~alpha2`]+2*D_[`~alpha2`](Physics[Christoffel][`~alpha1`, alpha1, alpha2], [X]))*Physics[g_][`~alpha`, `~beta`]-(1/4)*D_[`~rho1`](Physics[Christoffel][`~alpha`, rho1, `~beta`], [X])+(1/4)*D_[`~mu`](Physics[Christoffel][`~alpha`, mu, `~beta`], [X])+(1/4)*D_[`~nu`](Physics[Christoffel][`~beta`, nu, `~alpha`], [X])-(1/2)*D_[`~beta`](Physics[Christoffel][`~upsilon`, upsilon, `~alpha`], [X])-(1/4)*D_[`~omega`](Physics[Christoffel][`~beta`, omega, `~alpha`], [X])+(1/2)*D_[`~beta`](Physics[Christoffel][alpha9, `~alpha`, `~alpha9`], [X])-Physics[Ricci][`~alpha`, `~beta`] = 0)

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

`Systems of spacetime coordinates are:`*{X = (x, y, z, t)}

 

_______________________________________________________

 

[coordinatesystems = {X}]

(42)

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

g_[]

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 h[mu, nu], 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 delta[i, j](X) by delta[i, j])

h[mu, nu] = Matrix(4, proc (i, j) options operator, arrow; delta[i, j](X) end proc, shape = symmetric)

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

(44)

In the above it is understood that `&ll;`(`&delta;__i,j`, 1), so that quadratic or higher powers of it or its derivatives can be approximated to 0 (discarded). Define the components of `h__&mu;,&nu;` accordingly

"Define(?)"

`Defined objects with tensor properties`

 

{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 `&eta;__&mu;,&nu;` representing the unperturbed Minkowski metric

"eta[mu,nu] = rhs(?)"

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

(46)

"Define(?)"

`Defined objects with tensor properties`

 

{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

g_[mu, nu] = eta[mu, nu]+h[mu, nu]

Physics:-g_[mu, nu] = eta[mu, nu]+h[mu, nu]

(48)

Make this be the definition of the metric

Define(Physics[g_][mu, nu] = eta[mu, nu]+h[mu, nu])

_______________________________________________________

 

`Coordinates: `[x, y, z, t]*`. Signature: `(`- - - +`)

 

_______________________________________________________

 

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

 

_______________________________________________________

 

`Setting `*lowercaselatin_is*` letters to represent `*space*` indices`

 

`Defined objects with tensor properties`

 

{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

Ricci[definition]

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

(50)

Christoffel[`~alpha`, mu, nu, definition]

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)

Substitute(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])), Physics[Ricci][mu, nu] = Physics[d_][alpha](Physics[Christoffel][`~alpha`, mu, nu], [X])-Physics[d_][nu](Physics[Christoffel][`~alpha`, mu, alpha], [X])+Physics[Christoffel][`~beta`, mu, nu]*Physics[Christoffel][`~alpha`, beta, alpha]-Physics[Christoffel][`~beta`, mu, alpha]*Physics[Christoffel][`~alpha`, nu, beta])

Physics:-Ricci[mu, nu] = Physics:-d_[alpha]((1/2)*Physics:-g_[`~alpha`, `~kappa`]*(Physics:-d_[nu](Physics:-g_[kappa, mu], [X])+Physics:-d_[mu](Physics:-g_[kappa, nu], [X])-Physics:-d_[kappa](Physics:-g_[mu, nu], [X])), [X])-Physics:-d_[nu]((1/2)*Physics:-g_[`~alpha`, `~tau`]*(Physics:-d_[mu](Physics:-g_[tau, alpha], [X])+Physics:-d_[alpha](Physics:-g_[tau, mu], [X])-Physics:-d_[tau](Physics:-g_[alpha, mu], [X])), [X])+(1/4)*Physics:-g_[`~beta`, `~iota`]*(Physics:-d_[nu](Physics:-g_[iota, mu], [X])+Physics:-d_[mu](Physics:-g_[iota, nu], [X])-Physics:-d_[iota](Physics:-g_[mu, nu], [X]))*Physics:-g_[`~alpha`, `~lambda`]*(Physics:-d_[beta](Physics:-g_[lambda, alpha], [X])+Physics:-d_[alpha](Physics:-g_[lambda, beta], [X])-Physics:-d_[lambda](Physics:-g_[alpha, beta], [X]))-(1/4)*Physics:-g_[`~beta`, `~omega`]*(Physics:-d_[mu](Physics:-g_[omega, alpha], [X])+Physics:-d_[alpha](Physics:-g_[omega, mu], [X])-Physics:-d_[omega](Physics:-g_[alpha, mu], [X]))*Physics:-g_[`~alpha`, `~chi`]*(Physics:-d_[nu](Physics:-g_[chi, beta], [X])+Physics:-d_[beta](Physics:-g_[chi, nu], [X])-Physics:-d_[chi](Physics:-g_[beta, nu], [X]))

(52)

Introducing (48)g[mu, nu] = eta[mu, nu]+h[mu, nu], and also the inert form of the Ricci tensor to facilitate simplification some steps below,

Substitute(Physics[g_][mu, nu] = eta[mu, nu]+h[mu, nu], Ricci = %Ricci, Physics[Ricci][mu, nu] = Physics[d_][alpha]((1/2)*Physics[g_][`~alpha`, `~kappa`]*(Physics[d_][nu](Physics[g_][kappa, mu], [X])+Physics[d_][mu](Physics[g_][kappa, nu], [X])-Physics[d_][kappa](Physics[g_][mu, nu], [X])), [X])-Physics[d_][nu]((1/2)*Physics[g_][`~alpha`, `~tau`]*(Physics[d_][mu](Physics[g_][tau, alpha], [X])+Physics[d_][alpha](Physics[g_][tau, mu], [X])-Physics[d_][tau](Physics[g_][alpha, mu], [X])), [X])+(1/4)*Physics[g_][`~beta`, `~iota`]*(Physics[d_][nu](Physics[g_][iota, mu], [X])+Physics[d_][mu](Physics[g_][iota, nu], [X])-Physics[d_][iota](Physics[g_][mu, nu], [X]))*Physics[g_][`~alpha`, `~lambda`]*(Physics[d_][beta](Physics[g_][lambda, alpha], [X])+Physics[d_][alpha](Physics[g_][lambda, beta], [X])-Physics[d_][lambda](Physics[g_][alpha, beta], [X]))-(1/4)*Physics[g_][`~beta`, `~omega`]*(Physics[d_][mu](Physics[g_][omega, alpha], [X])+Physics[d_][alpha](Physics[g_][omega, mu], [X])-Physics[d_][omega](Physics[g_][alpha, mu], [X]))*Physics[g_][`~alpha`, `~chi`]*(Physics[d_][nu](Physics[g_][chi, beta], [X])+Physics[d_][beta](Physics[g_][chi, nu], [X])-Physics[d_][chi](Physics[g_][beta, nu], [X])))

%Ricci[mu, nu] = (1/2)*Physics:-d_[alpha](eta[`~alpha`, `~kappa`]+h[`~alpha`, `~kappa`], [X])*(Physics:-d_[nu](eta[kappa, mu]+h[kappa, mu], [X])+Physics:-d_[mu](eta[kappa, nu]+h[kappa, nu], [X])-Physics:-d_[kappa](eta[mu, nu]+h[mu, nu], [X]))+(1/2)*(eta[`~alpha`, `~kappa`]+h[`~alpha`, `~kappa`])*(Physics:-d_[alpha](Physics:-d_[nu](eta[kappa, mu]+h[kappa, mu], [X]), [X])+Physics:-d_[alpha](Physics:-d_[mu](eta[kappa, nu]+h[kappa, nu], [X]), [X])-Physics:-d_[alpha](Physics:-d_[kappa](eta[mu, nu]+h[mu, nu], [X]), [X]))-(1/2)*Physics:-d_[nu](eta[`~alpha`, `~tau`]+h[`~alpha`, `~tau`], [X])*(Physics:-d_[mu](eta[alpha, tau]+h[alpha, tau], [X])+Physics:-d_[alpha](eta[mu, tau]+h[mu, tau], [X])-Physics:-d_[tau](eta[alpha, mu]+h[alpha, mu], [X]))-(1/2)*(eta[`~alpha`, `~tau`]+h[`~alpha`, `~tau`])*(Physics:-d_[mu](Physics:-d_[nu](eta[alpha, tau]+h[alpha, tau], [X]), [X])+Physics:-d_[alpha](Physics:-d_[nu](eta[mu, tau]+h[mu, tau], [X]), [X])-Physics:-d_[nu](Physics:-d_[tau](eta[alpha, mu]+h[alpha, mu], [X]), [X]))+(1/4)*(eta[`~beta`, `~iota`]+h[`~beta`, `~iota`])*(Physics:-d_[nu](eta[iota, mu]+h[iota, mu], [X])+Physics:-d_[mu](eta[iota, nu]+h[iota, nu], [X])-Physics:-d_[iota](eta[mu, nu]+h[mu, nu], [X]))*(eta[`~alpha`, `~lambda`]+h[`~alpha`, `~lambda`])*(Physics:-d_[beta](eta[alpha, lambda]+h[alpha, lambda], [X])+Physics:-d_[alpha](eta[beta, lambda]+h[beta, lambda], [X])-Physics:-d_[lambda](eta[alpha, beta]+h[alpha, beta], [X]))-(1/4)*(eta[`~beta`, `~omega`]+h[`~beta`, `~omega`])*(Physics:-d_[mu](eta[alpha, omega]+h[alpha, omega], [X])+Physics:-d_[alpha](eta[mu, omega]+h[mu, omega], [X])-Physics:-d_[omega](eta[alpha, mu]+h[alpha, mu], [X]))*(eta[`~alpha`, `~chi`]+h[`~alpha`, `~chi`])*(Physics:-d_[nu](eta[beta, chi]+h[beta, chi], [X])+Physics:-d_[beta](eta[chi, nu]+h[chi, nu], [X])-Physics:-d_[chi](eta[beta, nu]+h[beta, nu], [X]))

(53)

This expression contains several terms quadratic in the small perturbation `h__&mu;,&nu;` 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

Library:-Linearize(%Ricci[mu, nu] = (1/2)*Physics[d_][alpha](eta[`~alpha`, `~kappa`]+h[`~alpha`, `~kappa`], [X])*(Physics[d_][nu](eta[kappa, mu]+h[kappa, mu], [X])+Physics[d_][mu](eta[kappa, nu]+h[kappa, nu], [X])-Physics[d_][kappa](eta[mu, nu]+h[mu, nu], [X]))+(1/2)*(eta[`~alpha`, `~kappa`]+h[`~alpha`, `~kappa`])*(Physics[d_][alpha](Physics[d_][nu](eta[kappa, mu]+h[kappa, mu], [X]), [X])+Physics[d_][alpha](Physics[d_][mu](eta[kappa, nu]+h[kappa, nu], [X]), [X])-Physics[d_][alpha](Physics[d_][kappa](eta[mu, nu]+h[mu, nu], [X]), [X]))-(1/2)*Physics[d_][nu](eta[`~alpha`, `~tau`]+h[`~alpha`, `~tau`], [X])*(Physics[d_][mu](eta[alpha, tau]+h[alpha, tau], [X])+Physics[d_][alpha](eta[mu, tau]+h[mu, tau], [X])-Physics[d_][tau](eta[alpha, mu]+h[alpha, mu], [X]))-(1/2)*(eta[`~alpha`, `~tau`]+h[`~alpha`, `~tau`])*(Physics[d_][mu](Physics[d_][nu](eta[alpha, tau]+h[alpha, tau], [X]), [X])+Physics[d_][alpha](Physics[d_][nu](eta[mu, tau]+h[mu, tau], [X]), [X])-Physics[d_][nu](Physics[d_][tau](eta[alpha, mu]+h[alpha, mu], [X]), [X]))+(1/4)*(eta[`~beta`, `~iota`]+h[`~beta`, `~iota`])*(Physics[d_][nu](eta[iota, mu]+h[iota, mu], [X])+Physics[d_][mu](eta[iota, nu]+h[iota, nu], [X])-Physics[d_][iota](eta[mu, nu]+h[mu, nu], [X]))*(eta[`~alpha`, `~lambda`]+h[`~alpha`, `~lambda`])*(Physics[d_][beta](eta[alpha, lambda]+h[alpha, lambda], [X])+Physics[d_][alpha](eta[beta, lambda]+h[beta, lambda], [X])-Physics[d_][lambda](eta[alpha, beta]+h[alpha, beta], [X]))-(1/4)*(eta[`~beta`, `~omega`]+h[`~beta`, `~omega`])*(Physics[d_][mu](eta[alpha, omega]+h[alpha, omega], [X])+Physics[d_][alpha](eta[mu, omega]+h[mu, omega], [X])-Physics[d_][omega](eta[alpha, mu]+h[alpha, mu], [X]))*(eta[`~alpha`, `~chi`]+h[`~alpha`, `~chi`])*(Physics[d_][nu](eta[beta, chi]+h[beta, chi], [X])+Physics[d_][beta](eta[chi, nu]+h[chi, nu], [X])-Physics[d_][chi](eta[beta, nu]+h[beta, nu], [X])), h)

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

(54)

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

g_[min]

_______________________________________________________

 

`The Minkowski metric in coordinates `*[x, y, z, t]

 

`Signature: `(`- - - +`)

 

_______________________________________________________

 

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

(55)

To proceed simplifying, replace in the expression (54) for the Ricci tensor the intermediate Minkowski `&eta;__&mu;,&nu;`by `#msub(mi("g"),mrow(mi("&mu;",fontstyle = "normal"),mo("&comma;"),mi("&nu;",fontstyle = "normal")))`

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

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

(56)

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

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

%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

Einstein[definition]

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

(58)

Compute the trace R[alpha, `~alpha`] directly from the linearized form (57) of the Ricci tensor,

g_[mu, nu]*(%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]))

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

Simplify(%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`])

%Ricci[nu, `~nu`] = -Physics:-dAlembertian(h[alpha, `~alpha`], [X])+Physics:-d_[alpha](Physics:-d_[tau](h[`~alpha`, `~tau`], [X]), [X])

(60)

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

(%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]))-(1/2)*g_[mu, nu]*(%Ricci[nu, `~nu`] = -Physics[dAlembertian](h[alpha, `~alpha`], [X])+Physics[d_][alpha](Physics[d_][tau](h[`~alpha`, `~tau`], [X]), [X]))

%Ricci[mu, nu]-(1/2)*Physics:-g_[mu, nu]*%Ricci[alpha, `~alpha`] = -(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])-(1/2)*Physics:-g_[mu, nu]*(-Physics:-dAlembertian(h[alpha, `~alpha`], [X])+Physics:-d_[alpha](Physics:-d_[tau](h[`~alpha`, `~tau`], [X]), [X]))

(61)

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


You can now redefine the general h[mu, nu] 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

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

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

(62)

Substitute(Physics[d_][mu](h[`~mu`, nu], [X]) = (1/2)*Physics[d_][nu](h[alpha, `~alpha`], [X]), %Ricci[mu, nu]-(1/2)*Physics[g_][mu, nu]*%Ricci[alpha, `~alpha`] = -(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])-(1/2)*Physics[g_][mu, nu]*(-Physics[dAlembertian](h[alpha, `~alpha`], [X])+Physics[d_][alpha](Physics[d_][tau](h[`~alpha`, `~tau`], [X]), [X])))

%Ricci[mu, nu]-(1/2)*Physics:-g_[mu, nu]*%Ricci[alpha, `~alpha`] = -(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]((1/2)*Physics:-d_[mu](h[lambda, `~lambda`], [X]), [X])+(1/2)*Physics:-d_[mu]((1/2)*Physics:-d_[nu](h[kappa, `~kappa`], [X]), [X])-(1/2)*Physics:-g_[mu, nu]*(-Physics:-dAlembertian(h[alpha, `~alpha`], [X])+Physics:-d_[alpha]((1/2)*Physics:-d_[`~alpha`](h[beta, `~beta`], [X]), [X]))

(63)

Simplify(%Ricci[mu, nu]-(1/2)*Physics[g_][mu, nu]*%Ricci[alpha, `~alpha`] = -(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]((1/2)*Physics[d_][mu](h[lambda, `~lambda`], [X]), [X])+(1/2)*Physics[d_][mu]((1/2)*Physics[d_][nu](h[kappa, `~kappa`], [X]), [X])-(1/2)*Physics[g_][mu, nu]*(-Physics[dAlembertian](h[alpha, `~alpha`], [X])+Physics[d_][alpha]((1/2)*Physics[d_][`~alpha`](h[beta, `~beta`], [X]), [X])))

%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]

_______________________________________________________

 

`Systems of spacetime coordinates are:`*{X = (r, theta, phi, t)}

 

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

 

`Setting `*lowercaselatin_is*` letters to represent `*space*` indices`

 

`The Schwarzschild metric in coordinates `*[r, theta, phi, t]

 

`Parameters: `[m]

 

`Signature: `(`- - - +`)

 

_______________________________________________________

 

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

(65)

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

Define(T[mu])

`Defined objects with tensor properties`

 

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

Define(R[mu], relativeweight = 1)

`Defined objects with tensor properties`

 

{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 r, theta, phi, t to  rho, theta, phi, t where

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

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

(68)

The transformed components of T[mu] and R[mu] are, respectively,

TransformCoordinates(TR, T[mu], [rho, theta, phi, t])

Vector[column](%id = 36893488151820263660)

(69)

TransformCoordinates(TR, R[mu], [rho, theta, phi, t])

Vector[column](%id = 36893488151823155676)

(70)

where, when comparing both results, we see that the transformed components for R[mu] are all multiplied by J^n with n = 1 and J is the determinant of the transformation:

J__matrix := simplify(VectorCalculus:-Jacobian([rhs(TR), theta, phi, t], [rho, theta, phi, t]))

Matrix(%id = 36893488151964220700)

(71)

J = LinearAlgebra:-Determinant(J__matrix)

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 T[mu] and R[mu] used above,

Library:-GetRelativeWeight(T[mu])

0

(73)

Library:-GetRelativeWeight(R[mu])

1

(74)

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

Library:-GetRelativeWeight(R[`~mu`])

1

(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

Library:-GetRelativeWeight(LeviCivita[alpha, beta, mu, nu])

-1

(76)

"Library:-GetRelativeWeight(LeviCivita[~alpha, ~beta, ~mu,~nu])   "

1

(77)

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

R[mu]^2

R[mu]*R[`~mu`]

(78)

Library:-GetRelativeWeight(R[mu]*R[`~mu`])

2

(79)

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

1/R[mu]^2

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

(80)

Library:-GetRelativeWeight(1/(R[mu]*R[`~mu`]))

-2

(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.

"R[~mu] + LeviCivita[~alpha, ~beta, ~mu,~nu]*T[alpha] T[beta] T[nu]"

T[alpha]*T[beta]*T[nu]*Physics:-LeviCivita[`~alpha`, `~beta`, `~mu`, `~nu`]+R[`~mu`]

(82)

Library:-GetRelativeWeight(T[alpha]*T[beta]*T[nu]*Physics[LeviCivita][`~alpha`, `~beta`, `~mu`, `~nu`]+R[`~mu`])

1

(83)

The relative weight of any determinant is always equal to 2

%g_[determinant]

%g_[determinant]

(84)

Library:-GetRelativeWeight(%g_[determinant])

2

(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.

g__det := %g_[:-determinant]

%g_[determinant]

(86)

Library:-GetRelativeTerm(g__det, mu)

-2*Physics:-Christoffel[`~nu`, mu, nu]*%g_[determinant]

(87)
  

Consequently,

(%D_[mu] = D_[mu])(g__det)

%D_[mu](%g_[determinant]) = 0

(88)
  

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

convert(%D_[mu](%g_[determinant]) = 0, d_)

%d_[mu](%g_[determinant], [X])-2*Physics:-Christoffel[`~alpha`, alpha, mu]*%g_[determinant] = 0

(89)
  

evaluate the inert %d_

factor(eval(%d_[mu](%g_[determinant], [X])-2*Physics[Christoffel][`~alpha`, alpha, mu]*%g_[determinant] = 0, %d_ = d_))

%g_[determinant]*(Physics:-g_[`~alpha`, `~nu`]*Physics:-d_[mu](Physics:-g_[alpha, nu], [X])-2*Physics:-Christoffel[`~alpha`, alpha, mu]) = 0

(90)
  

The factor in parentheses is equal to `#mrow(msubsup(mi("g"),none(),mrow(mi("&alpha;",fontstyle = "normal"),mo("&comma;"),mi("&nu;",fontstyle = "normal"))),msub(mi("&dtri;"),mi("&mu;",fontstyle = "normal")),mo("&ApplyFunction;"),mfenced(msub(mi("g"),mrow(mi("&alpha;",fontstyle = "normal"),mo("&comma;"),mi("&nu;",fontstyle = "normal")))))`, where the covariant derivative of the metric is equal to zero, so

Simplify(%g_[determinant]*(Physics[g_][`~alpha`, `~nu`]*Physics[d_][mu](Physics[g_][alpha, nu], [X])-2*Physics[Christoffel][`~alpha`, alpha, mu]) = 0)

0 = 0

(91)
  

Consider the covariant derivative of T[mu] and R[mu] defined in (66) and (67)

Library:-GetRelativeWeight(T[mu])

0

(92)

Library:-GetRelativeWeight(R[mu])

1

(93)

The corresponding covariant derivatives

(%D_[mu] = D_[mu])(T[mu](X))

%D_[mu](T[`~mu`](X)) = Physics:-D_[mu](T[`~mu`](X), [X])

(94)

expand(%D_[mu](T[`~mu`](X)) = D_[mu](T[`~mu`](X), [X]))

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

(%D_[mu] = D_[mu])(R[mu](X))

%D_[mu](R[`~mu`](X)) = Physics:-D_[mu](R[`~mu`](X), [X])

(96)

expand(%D_[mu](R[`~mu`](X)) = D_[mu](R[`~mu`](X), [X]))

%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

Library:-GetRelativeTerm(R[`~mu`](X), mu)

-Physics:-Christoffel[`~nu`, mu, nu]*R[`~mu`](X)

(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

Christoffel[mu, alpha, beta] = Library:-ConvertToF(Christoffel[mu, alpha, beta], g_)

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 A[mu] representing the 4D electromagnetic potential as a function of the coordinates X and F[mu, nu] representing the electromagnetic field tensors

Define(A[mu] = A[mu](X), quiet)

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

Define(F[mu, nu] = d_[mu](A[nu])-d_[nu](A[mu]))

`Defined objects with tensor properties`

 

{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 A[mu]

F[mu, nu] = Library:-ConvertToF(F[mu, nu], A)

F[mu, nu] = Physics:-d_[mu](A[nu], [X])-Physics:-d_[nu](A[mu], [X])

(103)

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

F[definition]

F[mu, nu] = Physics:-d_[mu](A[nu], [X])-Physics:-d_[nu](A[mu], [X])

(104)

The rewriting, however, works also with tensorial expressions

F[mu, nu]*A[mu]*A[nu]

F[mu, nu]*A[`~mu`]*A[`~nu`]

(105)

Library:-ConvertToF(F[mu, nu]*A[`~mu`]*A[`~nu`], A)

(Physics:-d_[mu](A[nu], [X])-Physics:-d_[nu](A[mu], [X]))*A[`~mu`]*A[`~nu`]

(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

NULL


 

Download New_in_Physics_2025.mw

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

Computing Einstein's equations using variational principles

Edgardo S. Cheb-Terrab

Freddy Baudine

 

One of the biggest challenges for a computer algebra system is to compute Einstein's equations from first principles, by equating to zero the functional derivative of the Action for gravity, without using human-shortcuts or tricks. Developments during 2024, appearing as new in Maple 2025, broke through that barrier and now the LagrangeEquations  Physics command, introduced in 2023, can perform that elusive computation of Einstein's equations, not using tabulated cases, properly handling several (traditional or not) alternative ways of presenting the Lagrangian (the integrand in the Action), taking advantage of the functional differentiation  capabilities of the Physics  package .

 

The computation can be performed in one call to Physics:-LagrangeEquations, or in steps interactively using the Physics:-Fundiff  and Physics:-Simplify  commands that include newly implemented capabilities for simplifying tensorial expressions in curved spacetimes. This is an exciting breakthrough/milestone in Computer Algebra, also in an area out of reach of neural network AIs.

 

Einstein's equations are a system of second order nonlinear coupled partial differential equations for the 10 components of the spacetime metric g[mu, nu]

with(Physics); Setup(coordinates = cartesian, metric = arbitrary)

`Systems of spacetime coordinates are:`*{X = (x, y, z, t)}

 

_______________________________________________________

 

`Setting `*lowercaselatin_is*` letters to represent `*space*` indices`

 

`The arbitrary metric in coordinates `*[x, y, z, t]

 

`Signature: `(`- - - +`)

 

_______________________________________________________

 

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

 

_______________________________________________________; "_noterminate"

(1)

In the Lagrangian formulation, the coordinates of the problem are the 10 components of the metric,

CompactDisplay(g_[])

f__1(x, y, z, t)*`will now be displayed as`*f__1

 

f__10(x, y, z, t)*`will now be displayed as`*f__10

 

f__2(x, y, z, t)*`will now be displayed as`*f__2

 

f__3(x, y, z, t)*`will now be displayed as`*f__3

 

f__4(x, y, z, t)*`will now be displayed as`*f__4

 

f__5(x, y, z, t)*`will now be displayed as`*f__5

 

f__6(x, y, z, t)*`will now be displayed as`*f__6

 

f__7(x, y, z, t)*`will now be displayed as`*f__7

 

f__8(x, y, z, t)*`will now be displayed as`*f__8

 

f__9(x, y, z, t)*`will now be displayed as`*f__9

(2)

and the parameters of the variational problem are the spacetime coordinates X^alpha.

 

The traditional Lagrangian

 

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 and the determinant of the metric by

L := sqrt(-%g_[determinant])*Ricci[alpha, `~alpha`]

(-%g_[determinant])^(1/2)*Physics:-Ricci[alpha, `~alpha`]

(3)

What was almost a dream in previous years, Einstein's equations can now be computed in one simple instruction as the Lagrange equations for this Lagrangian, in the traditional compact form, taking the components of the metric tensor as the coordinates (for an interactive, step by step computation, see further below)

LagrangeEquations(L, g_[mu, nu])

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

(4)

The tensorial equation computed is also the definition of the Einstein  tensor

Einstein[definition]

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

(5)

A Lagrangian depending only on first derivatives of the metric

 

The Lagrangian L used to compute Einstein's equations (4)  contains first and second derivatives of the metric; the latter make the computation significantly more complicated. The second order derivatives, however, can be removed from the formulation. To see that, rewrite L in terms of Christoffel  symbols

L__C := convert(L, Christoffel)

(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*(Physics:-d_[nu](Physics:-Christoffel[`~nu`, alpha, lambda], [X])-Physics:-d_[lambda](Physics:-Christoffel[`~nu`, alpha, nu], [X])+Physics:-Christoffel[`~beta`, alpha, lambda]*Physics:-Christoffel[`~nu`, beta, nu]-Physics:-Christoffel[`~beta`, alpha, nu]*Physics:-Christoffel[`~nu`, beta, lambda])

(6)

Recalling the definition

Christoffel[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])

(7)

in L[C] the two terms containing derivatives of Christoffel symbols contain second order derivatives of g[mu, nu].

 

Now, it is always possible to add a total spacetime derivative to L[C] 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 L[C], obtain a Lagrangian involving only 1st order derivatives. The total derivative, expressed using the inert `&PartialD;` command to see it before the differentiation operation is performed, is

TD := %d_[alpha](g_[`~mu`, `~nu`]*sqrt(-%g_[determinant])*(g_[`~alpha`, mu]*Christoffel[`~beta`, nu, beta]-Christoffel[`~alpha`, mu, nu]))

%d_[alpha](Physics:-g_[`~mu`, `~nu`]*(-%g_[determinant])^(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-Christoffel[`~beta`, beta, nu]-Physics:-Christoffel[`~alpha`, mu, nu]))

(8)

Adding this term to L[C], performing the `&PartialD;` differentiation operation and simplifying we get

L__1 := L__C+TD

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

(9)

L__1 := eval(L__1, %d_ = d_)

(-%g_[determinant])^(1/2)*Physics:-g_[`~alpha`, `~lambda`]*(Physics:-d_[nu](Physics:-Christoffel[`~nu`, alpha, lambda], [X])-Physics:-d_[lambda](Physics:-Christoffel[`~nu`, alpha, nu], [X])+Physics:-Christoffel[`~beta`, alpha, lambda]*Physics:-Christoffel[`~nu`, beta, nu]-Physics:-Christoffel[`~beta`, alpha, nu]*Physics:-Christoffel[`~nu`, beta, lambda])+Physics:-d_[alpha](Physics:-g_[`~mu`, `~nu`], [X])*(-%g_[determinant])^(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-Christoffel[`~beta`, beta, nu]-Physics:-Christoffel[`~alpha`, mu, nu])-(1/2)*Physics:-g_[`~mu`, `~nu`]*(Physics:-g_[mu, `~alpha`]*Physics:-Christoffel[`~beta`, beta, nu]-Physics:-Christoffel[`~alpha`, mu, nu])*%g_[determinant]*Physics:-g_[`~kappa`, `~lambda`]*Physics:-d_[alpha](Physics:-g_[kappa, lambda], [X])/(-%g_[determinant])^(1/2)+Physics:-g_[`~mu`, `~nu`]*(-%g_[determinant])^(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-d_[alpha](Physics:-Christoffel[`~beta`, beta, nu], [X])-Physics:-d_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X]))

(10)

L__1 := Simplify(L__1)

(Physics:-Christoffel[alpha, beta, kappa]*Physics:-Christoffel[`~beta`, `~alpha`, `~kappa`]-Physics:-Christoffel[alpha, beta, `~alpha`]*Physics:-Christoffel[`~beta`, kappa, `~kappa`])*(-%g_[determinant])^(1/2)

(11)

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 (4); this computation can also now be performed in a single call

LagrangeEquations(L__1, g_[mu, nu])

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

(12)

Simplification capabilities developed to cover these computations

 

To illustrate the new Maple 2025 tensorial simplification capabilities note that `&equiv;`(L[1], (Physics[Christoffel][alpha, beta, kappa]*Physics[Christoffel][`~beta`, `~alpha`, `~kappa`]-Physics[Christoffel][alpha, beta, `~alpha`]*Physics[Christoffel][`~beta`, kappa, `~kappa`])*(-%g_[determinant])^(1/2)) is no just L[C] ≡ (6) after discarding its two terms involving derivatives of Christoffel symbols. To verify this, split L[C] into the terms containing or not derivatives of Christoffel

L__22, L__11 := selectremove(has, expand(L__C), d_)

(-%g_[determinant])^(1/2)*Physics:-d_[nu](Physics:-Christoffel[`~nu`, alpha, lambda], [X])*Physics:-g_[`~alpha`, `~lambda`]-(-%g_[determinant])^(1/2)*Physics:-d_[lambda](Physics:-Christoffel[`~nu`, alpha, nu], [X])*Physics:-g_[`~alpha`, `~lambda`], (-%g_[determinant])^(1/2)*Physics:-Christoffel[`~beta`, alpha, lambda]*Physics:-Christoffel[`~nu`, beta, nu]*Physics:-g_[`~alpha`, `~lambda`]-(-%g_[determinant])^(1/2)*Physics:-Christoffel[`~beta`, alpha, nu]*Physics:-Christoffel[`~nu`, beta, lambda]*Physics:-g_[`~alpha`, `~lambda`]

(13)

Comparing, the total derivative TD≡ (8) is not just -L[22], but

TD = -L__22-2*L__11

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

(14)

Things like these, TD = -L__22-2*L__11, 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 `&PartialD;` and simplify to see the equality is true

(lhs-rhs)(%d_[alpha](Physics[g_][`~mu`, `~nu`]*(-%g_[determinant])^(1/2)*(Physics[g_][mu, `~alpha`]*Physics[Christoffel][`~beta`, beta, nu]-Physics[Christoffel][`~alpha`, mu, nu])) = -(-%g_[determinant])^(1/2)*Physics[d_][nu](Physics[Christoffel][`~nu`, alpha, lambda], [X])*Physics[g_][`~alpha`, `~lambda`]+(-%g_[determinant])^(1/2)*Physics[d_][lambda](Physics[Christoffel][`~nu`, alpha, nu], [X])*Physics[g_][`~alpha`, `~lambda`]-2*(-%g_[determinant])^(1/2)*Physics[Christoffel][`~beta`, alpha, lambda]*Physics[Christoffel][`~nu`, beta, nu]*Physics[g_][`~alpha`, `~lambda`]+2*(-%g_[determinant])^(1/2)*Physics[Christoffel][`~beta`, alpha, nu]*Physics[Christoffel][`~nu`, beta, lambda]*Physics[g_][`~alpha`, `~lambda`])

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

(15)

eval(%d_[alpha](Physics[g_][`~mu`, `~nu`]*(-%g_[determinant])^(1/2)*(Physics[g_][mu, `~alpha`]*Physics[Christoffel][`~beta`, beta, nu]-Physics[Christoffel][`~alpha`, mu, nu]))+(-%g_[determinant])^(1/2)*Physics[d_][nu](Physics[Christoffel][`~nu`, alpha, lambda], [X])*Physics[g_][`~alpha`, `~lambda`]-(-%g_[determinant])^(1/2)*Physics[d_][lambda](Physics[Christoffel][`~nu`, alpha, nu], [X])*Physics[g_][`~alpha`, `~lambda`]+2*(-%g_[determinant])^(1/2)*Physics[Christoffel][`~beta`, alpha, lambda]*Physics[Christoffel][`~nu`, beta, nu]*Physics[g_][`~alpha`, `~lambda`]-2*(-%g_[determinant])^(1/2)*Physics[Christoffel][`~beta`, alpha, nu]*Physics[Christoffel][`~nu`, beta, lambda]*Physics[g_][`~alpha`, `~lambda`], %d_ = d_)

(-%g_[determinant])^(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-Christoffel[`~beta`, beta, nu]-Physics:-Christoffel[`~alpha`, mu, nu])*Physics:-d_[alpha](Physics:-g_[`~mu`, `~nu`], [X])-(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-Christoffel[`~beta`, beta, nu]-Physics:-Christoffel[`~alpha`, mu, nu])*Physics:-g_[`~mu`, `~nu`]*%g_[determinant]*Physics:-g_[`~kappa`, `~lambda`]*Physics:-d_[alpha](Physics:-g_[kappa, lambda], [X])/(-%g_[determinant])^(1/2)+(-%g_[determinant])^(1/2)*(Physics:-g_[mu, `~alpha`]*Physics:-d_[alpha](Physics:-Christoffel[`~beta`, beta, nu], [X])-Physics:-d_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X]))*Physics:-g_[`~mu`, `~nu`]+(-%g_[determinant])^(1/2)*Physics:-d_[nu](Physics:-Christoffel[`~nu`, alpha, lambda], [X])*Physics:-g_[`~alpha`, `~lambda`]-(-%g_[determinant])^(1/2)*Physics:-d_[lambda](Physics:-Christoffel[`~nu`, alpha, nu], [X])*Physics:-g_[`~alpha`, `~lambda`]+2*(-%g_[determinant])^(1/2)*Physics:-Christoffel[`~beta`, alpha, lambda]*Physics:-Christoffel[`~nu`, beta, nu]*Physics:-g_[`~alpha`, `~lambda`]-2*(-%g_[determinant])^(1/2)*Physics:-Christoffel[`~beta`, alpha, nu]*Physics:-Christoffel[`~nu`, beta, lambda]*Physics:-g_[`~alpha`, `~lambda`]

(16)

Simplify((-%g_[determinant])^(1/2)*(Physics[g_][mu, `~alpha`]*Physics[Christoffel][`~beta`, beta, nu]-Physics[Christoffel][`~alpha`, mu, nu])*Physics[d_][alpha](Physics[g_][`~mu`, `~nu`], [X])-(1/2)*(Physics[g_][mu, `~alpha`]*Physics[Christoffel][`~beta`, beta, nu]-Physics[Christoffel][`~alpha`, mu, nu])*Physics[g_][`~mu`, `~nu`]*%g_[determinant]*Physics[g_][`~kappa`, `~lambda`]*Physics[d_][alpha](Physics[g_][kappa, lambda], [X])/(-%g_[determinant])^(1/2)+(-%g_[determinant])^(1/2)*(Physics[g_][mu, `~alpha`]*Physics[d_][alpha](Physics[Christoffel][`~beta`, beta, nu], [X])-Physics[d_][alpha](Physics[Christoffel][`~alpha`, mu, nu], [X]))*Physics[g_][`~mu`, `~nu`]+(-%g_[determinant])^(1/2)*Physics[d_][nu](Physics[Christoffel][`~nu`, alpha, lambda], [X])*Physics[g_][`~alpha`, `~lambda`]-(-%g_[determinant])^(1/2)*Physics[d_][lambda](Physics[Christoffel][`~nu`, alpha, nu], [X])*Physics[g_][`~alpha`, `~lambda`]+2*(-%g_[determinant])^(1/2)*Physics[Christoffel][`~beta`, alpha, lambda]*Physics[Christoffel][`~nu`, beta, nu]*Physics[g_][`~alpha`, `~lambda`]-2*(-%g_[determinant])^(1/2)*Physics[Christoffel][`~beta`, alpha, nu]*Physics[Christoffel][`~nu`, beta, lambda]*Physics[g_][`~alpha`, `~lambda`])

0

(17)

That said, it is also true that TD = -L[22]-2*L[11] results in the Lagrangian L[1] = -L[11], and since the equations of movement don't depend on the sign of the Lagrangian, for this Lagrangian `&equiv;`(L[C], (-%g_[determinant])^(1/2)*Physics[g_][`~alpha`, `~lambda`]*(Physics[d_][nu](Physics[Christoffel][`~nu`, alpha, lambda], [X])-Physics[d_][lambda](Physics[Christoffel][`~nu`, alpha, nu], [X])+Physics[Christoffel][`~beta`, alpha, lambda]*Physics[Christoffel][`~nu`, beta, nu]-Physics[Christoffel][`~beta`, alpha, nu]*Physics[Christoffel][`~nu`, beta, lambda])) adding the term TD happens to be equivalent to just discarding the terms of L__C involving derivatives of Christoffel symbols.

 

Derivation step by step using functional differentiation (variational principle)

 

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

S := Intc(L, X)

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

(18)

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

EE__unsimplified := Fundiff(S, g_[alpha, beta]) = 0

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

EEC := Simplify(EE__unsimplified)

(1/4)*(2*Physics:-Christoffel[chi, iota, kappa]*Physics:-Christoffel[`~iota`, `~chi`, `~kappa`]-2*Physics:-Christoffel[chi, iota, `~chi`]*Physics:-Christoffel[`~iota`, kappa, `~kappa`]-2*Physics:-D_[iota](Physics:-Christoffel[chi, `~chi`, `~iota`], [X])+2*Physics:-D_[chi](Physics:-Christoffel[`~chi`, iota, `~iota`], [X]))*Physics:-g_[`~alpha`, `~beta`]+(1/4)*(Physics:-Christoffel[`~alpha`, chi, `~beta`]+Physics:-Christoffel[`~beta`, chi, `~alpha`])*Physics:-Christoffel[`~chi`, iota, `~iota`]-(1/4)*Physics:-Christoffel[`~beta`, chi, iota]*Physics:-Christoffel[`~chi`, `~alpha`, `~iota`]-(1/2)*Physics:-Christoffel[chi, iota, `~alpha`]*Physics:-Christoffel[`~iota`, `~beta`, `~chi`]+(1/2)*Physics:-Christoffel[chi, `~alpha`, `~beta`]*Physics:-Christoffel[iota, `~chi`, `~iota`]-(1/4)*Physics:-Christoffel[`~alpha`, chi, iota]*Physics:-Christoffel[`~chi`, `~beta`, `~iota`]+(1/4)*Physics:-D_[chi](Physics:-Christoffel[`~alpha`, `~beta`, `~chi`], [X])+(1/4)*Physics:-D_[chi](Physics:-Christoffel[`~beta`, `~alpha`, `~chi`], [X])-(1/2)*Physics:-D_[chi](Physics:-Christoffel[`~chi`, `~alpha`, `~beta`], [X])-(1/4)*Physics:-D_[`~alpha`](Physics:-Christoffel[`~beta`, chi, `~chi`], [X])+(1/2)*Physics:-D_[`~beta`](Physics:-Christoffel[chi, `~alpha`, `~chi`], [X])-(1/4)*Physics:-D_[`~beta`](Physics:-Christoffel[`~alpha`, chi, `~chi`], [X]) = 0

(19)

In this result, we see`&dtri;` 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,

"`&dtri;`[chi](GAMMA[]^(alpha,beta,chi)) :"

% = expand(%)

Physics:-D_[chi](Physics:-Christoffel[`~alpha`, `~beta`, `~chi`], [X]) = Physics:-d_[chi](Physics:-Christoffel[`~alpha`, `~beta`, `~chi`], [X])+Physics:-Christoffel[`~alpha`, chi, mu]*Physics:-Christoffel[`~mu`, `~beta`, `~chi`]+Physics:-Christoffel[`~beta`, chi, mu]*Physics:-Christoffel[`~alpha`, `~chi`, `~mu`]+Physics:-Christoffel[`~chi`, chi, mu]*Physics:-Christoffel[`~alpha`, `~beta`, `~mu`]

(20)

With this computational meaning for the`&dtri;` derivatives of Christoffel symbols appearing in (19), rewrite EEC(19) in terms of the Ricci  and Riemann  tensors. For that, consider the definition

Ricci[definition]

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

(21)

Rewrite the noncovariant derivatives `&PartialD;` in terms of`&dtri;` derivatives using the computational representation (20), simplify and isolate one of them

convert(Physics[Ricci][mu, nu] = Physics[d_][alpha](Physics[Christoffel][`~alpha`, mu, nu], [X])-Physics[d_][nu](Physics[Christoffel][`~alpha`, mu, alpha], [X])+Physics[Christoffel][`~beta`, mu, nu]*Physics[Christoffel][`~alpha`, beta, alpha]-Physics[Christoffel][`~beta`, mu, alpha]*Physics[Christoffel][`~alpha`, nu, beta], D_)

Physics:-Ricci[mu, nu] = Physics:-D_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X])-Physics:-Christoffel[`~alpha`, alpha, kappa]*Physics:-Christoffel[`~kappa`, mu, nu]+Physics:-Christoffel[`~kappa`, alpha, mu]*Physics:-Christoffel[`~alpha`, kappa, nu]+Physics:-Christoffel[`~kappa`, alpha, nu]*Physics:-Christoffel[`~alpha`, mu, kappa]-Physics:-D_[nu](Physics:-Christoffel[`~alpha`, alpha, mu], [X])-Physics:-Christoffel[`~lambda`, mu, nu]*Physics:-Christoffel[`~alpha`, alpha, lambda]+Physics:-Christoffel[`~beta`, mu, nu]*Physics:-Christoffel[`~alpha`, alpha, beta]-Physics:-Christoffel[`~beta`, alpha, mu]*Physics:-Christoffel[`~alpha`, beta, nu]

(22)

Simplify(Physics[Ricci][mu, nu] = D_[alpha](Physics[Christoffel][`~alpha`, mu, nu], [X])-Physics[Christoffel][`~alpha`, alpha, kappa]*Physics[Christoffel][`~kappa`, mu, nu]+Physics[Christoffel][`~kappa`, alpha, mu]*Physics[Christoffel][`~alpha`, kappa, nu]+Physics[Christoffel][`~kappa`, alpha, nu]*Physics[Christoffel][`~alpha`, mu, kappa]-D_[nu](Physics[Christoffel][`~alpha`, alpha, mu], [X])-Physics[Christoffel][`~lambda`, mu, nu]*Physics[Christoffel][`~alpha`, alpha, lambda]+Physics[Christoffel][`~beta`, mu, nu]*Physics[Christoffel][`~alpha`, alpha, beta]-Physics[Christoffel][`~beta`, alpha, mu]*Physics[Christoffel][`~alpha`, beta, nu])

Physics:-Ricci[mu, nu] = Physics:-Christoffel[alpha, beta, mu]*Physics:-Christoffel[`~beta`, nu, `~alpha`]-Physics:-Christoffel[beta, mu, nu]*Physics:-Christoffel[alpha, `~alpha`, `~beta`]+Physics:-D_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X])-Physics:-D_[nu](Physics:-Christoffel[alpha, mu, `~alpha`], [X])

(23)

C_to_Ricci := isolate(Physics[Ricci][mu, nu] = Physics[Christoffel][alpha, beta, mu]*Physics[Christoffel][`~beta`, nu, `~alpha`]-Physics[Christoffel][beta, mu, nu]*Physics[Christoffel][alpha, `~alpha`, `~beta`]+D_[alpha](Physics[Christoffel][`~alpha`, mu, nu], [X])-D_[nu](Physics[Christoffel][alpha, mu, `~alpha`], [X]), D_[alpha](Christoffel[`~alpha`, mu, nu]))

Physics:-D_[alpha](Physics:-Christoffel[`~alpha`, mu, nu], [X]) = -Physics:-Christoffel[alpha, beta, mu]*Physics:-Christoffel[`~beta`, nu, `~alpha`]+Physics:-Christoffel[beta, mu, nu]*Physics:-Christoffel[alpha, `~alpha`, `~beta`]+Physics:-Ricci[mu, nu]+Physics:-D_[nu](Physics:-Christoffel[alpha, mu, `~alpha`], [X])

(24)

Analogously, derive an expression to rewrite`&dtri;` derivatives of Christoffel symbols using the Riemann  tensor

Riemann[`~alpha`, beta, mu, nu, definition]

Physics:-Riemann[`~alpha`, beta, mu, nu] = Physics:-d_[mu](Physics:-Christoffel[`~alpha`, beta, nu], [X])-Physics:-d_[nu](Physics:-Christoffel[`~alpha`, beta, mu], [X])+Physics:-Christoffel[`~alpha`, upsilon, mu]*Physics:-Christoffel[`~upsilon`, beta, nu]-Physics:-Christoffel[`~alpha`, upsilon, nu]*Physics:-Christoffel[`~upsilon`, beta, mu]

(25)

convert(Physics[Riemann][`~alpha`, beta, mu, nu] = Physics[d_][mu](Physics[Christoffel][`~alpha`, beta, nu], [X])-Physics[d_][nu](Physics[Christoffel][`~alpha`, beta, mu], [X])+Physics[Christoffel][`~alpha`, upsilon, mu]*Physics[Christoffel][`~upsilon`, beta, nu]-Physics[Christoffel][`~alpha`, upsilon, nu]*Physics[Christoffel][`~upsilon`, beta, mu], D_)

Physics:-Riemann[`~alpha`, beta, mu, nu] = Physics:-D_[mu](Physics:-Christoffel[`~alpha`, beta, nu], [X])+Physics:-Christoffel[`~kappa`, mu, nu]*Physics:-Christoffel[`~alpha`, beta, kappa]-Physics:-Christoffel[`~alpha`, kappa, mu]*Physics:-Christoffel[`~kappa`, beta, nu]+Physics:-Christoffel[`~kappa`, beta, mu]*Physics:-Christoffel[`~alpha`, kappa, nu]-Physics:-D_[nu](Physics:-Christoffel[`~alpha`, beta, mu], [X])-Physics:-Christoffel[`~lambda`, mu, nu]*Physics:-Christoffel[`~alpha`, beta, lambda]-Physics:-Christoffel[`~lambda`, beta, nu]*Physics:-Christoffel[`~alpha`, lambda, mu]+Physics:-Christoffel[`~alpha`, lambda, nu]*Physics:-Christoffel[`~lambda`, beta, mu]+Physics:-Christoffel[`~alpha`, mu, upsilon]*Physics:-Christoffel[`~upsilon`, beta, nu]-Physics:-Christoffel[`~alpha`, nu, upsilon]*Physics:-Christoffel[`~upsilon`, beta, mu]

(26)

Simplify(Physics[Riemann][`~alpha`, beta, mu, nu] = D_[mu](Physics[Christoffel][`~alpha`, beta, nu], [X])+Physics[Christoffel][`~kappa`, mu, nu]*Physics[Christoffel][`~alpha`, beta, kappa]-Physics[Christoffel][`~alpha`, kappa, mu]*Physics[Christoffel][`~kappa`, beta, nu]+Physics[Christoffel][`~kappa`, beta, mu]*Physics[Christoffel][`~alpha`, kappa, nu]-D_[nu](Physics[Christoffel][`~alpha`, beta, mu], [X])-Physics[Christoffel][`~lambda`, mu, nu]*Physics[Christoffel][`~alpha`, beta, lambda]-Physics[Christoffel][`~lambda`, beta, nu]*Physics[Christoffel][`~alpha`, lambda, mu]+Physics[Christoffel][`~alpha`, lambda, nu]*Physics[Christoffel][`~lambda`, beta, mu]+Physics[Christoffel][`~alpha`, mu, upsilon]*Physics[Christoffel][`~upsilon`, beta, nu]-Physics[Christoffel][`~alpha`, nu, upsilon]*Physics[Christoffel][`~upsilon`, beta, mu])

Physics:-Riemann[`~alpha`, beta, mu, nu] = -Physics:-Christoffel[`~alpha`, kappa, mu]*Physics:-Christoffel[`~kappa`, beta, nu]+Physics:-Christoffel[`~kappa`, beta, mu]*Physics:-Christoffel[`~alpha`, kappa, nu]+Physics:-D_[mu](Physics:-Christoffel[`~alpha`, beta, nu], [X])-Physics:-D_[nu](Physics:-Christoffel[`~alpha`, beta, mu], [X])

(27)

C_to_Riemann := isolate(Physics[Riemann][`~alpha`, beta, mu, nu] = -Physics[Christoffel][`~alpha`, kappa, mu]*Physics[Christoffel][`~kappa`, beta, nu]+Physics[Christoffel][`~kappa`, beta, mu]*Physics[Christoffel][`~alpha`, kappa, nu]+D_[mu](Physics[Christoffel][`~alpha`, beta, nu], [X])-D_[nu](Physics[Christoffel][`~alpha`, beta, mu], [X]), D_[mu](Christoffel[`~alpha`, beta, nu]))

Physics:-D_[mu](Physics:-Christoffel[`~alpha`, beta, nu], [X]) = Physics:-Christoffel[`~alpha`, kappa, mu]*Physics:-Christoffel[`~kappa`, beta, nu]-Physics:-Christoffel[`~kappa`, beta, mu]*Physics:-Christoffel[`~alpha`, kappa, nu]+Physics:-Riemann[`~alpha`, beta, mu, nu]+Physics:-D_[nu](Physics:-Christoffel[`~alpha`, beta, mu], [X])

(28)

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

Substitute(C_to_Riemann, C_to_Ricci, EEC)

(1/4)*(2*Physics:-Christoffel[chi, iota, kappa]*Physics:-Christoffel[`~iota`, `~chi`, `~kappa`]-2*Physics:-Christoffel[chi, iota, `~chi`]*Physics:-Christoffel[`~iota`, kappa, `~kappa`]-2*Physics:-Christoffel[alpha4, alpha5, alpha6]*Physics:-Christoffel[`~alpha6`, `~alpha4`, `~alpha5`]+2*Physics:-Christoffel[alpha4, alpha6, `~alpha5`]*Physics:-Christoffel[`~alpha6`, alpha5, `~alpha4`]-2*Physics:-D_[`~alpha5`](Physics:-Christoffel[alpha4, alpha5, `~alpha4`], [X])+2*Physics:-Christoffel[`~alpha1`, alpha1, alpha3]*Physics:-Christoffel[`~alpha3`, alpha2, `~alpha2`]-2*Physics:-Christoffel[`~alpha1`, alpha3, `~alpha2`]*Physics:-Christoffel[`~alpha3`, alpha1, alpha2]+2*Physics:-Ricci[alpha2, `~alpha2`]+2*Physics:-D_[`~alpha2`](Physics:-Christoffel[`~alpha1`, alpha1, alpha2], [X]))*Physics:-g_[`~alpha`, `~beta`]-Physics:-Ricci[`~alpha`, `~beta`]+(1/2)*Physics:-Christoffel[alpha10, `~alpha`, `~beta`]*Physics:-Christoffel[alpha9, `~alpha10`, `~alpha9`]-(1/2)*Physics:-Christoffel[alpha9, alpha10, `~alpha`]*Physics:-Christoffel[`~alpha10`, `~alpha9`, `~beta`]-(1/4)*Physics:-Christoffel[`~alpha`, chi, iota]*Physics:-Christoffel[`~chi`, `~beta`, `~iota`]-(1/4)*Physics:-Christoffel[`~beta`, chi, iota]*Physics:-Christoffel[`~chi`, `~alpha`, `~iota`]+(1/2)*Physics:-Christoffel[chi, `~alpha`, `~beta`]*Physics:-Christoffel[iota, `~chi`, `~iota`]+(1/4)*(Physics:-Christoffel[`~alpha`, chi, `~beta`]+Physics:-Christoffel[`~beta`, chi, `~alpha`])*Physics:-Christoffel[`~chi`, iota, `~iota`]-(1/2)*Physics:-Christoffel[chi, iota, `~alpha`]*Physics:-Christoffel[`~iota`, `~beta`, `~chi`]-(1/4)*Physics:-Christoffel[`~alpha`, rho, `~beta`]*Physics:-Christoffel[`~rho`, rho1, `~rho1`]+(1/4)*Physics:-Christoffel[`~alpha`, rho, `~rho1`]*Physics:-Christoffel[`~rho`, rho1, `~beta`]-(1/2)*Physics:-Christoffel[omicron, zeta, `~omicron`]*Physics:-Christoffel[`~zeta`, `~alpha`, `~beta`]+(1/2)*Physics:-Christoffel[omicron, zeta, `~beta`]*Physics:-Christoffel[`~zeta`, `~alpha`, `~omicron`]-(1/4)*Physics:-Christoffel[`~beta`, psi, `~alpha`]*Physics:-Christoffel[`~psi`, omega, `~omega`]+(1/4)*Physics:-Christoffel[`~beta`, psi, `~omega`]*Physics:-Christoffel[`~psi`, omega, `~alpha`]+(1/2)*Physics:-Christoffel[`~tau`, upsilon, `~alpha`]*Physics:-Christoffel[`~upsilon`, tau, `~beta`]-(1/2)*Physics:-Christoffel[`~tau`, `~alpha`, `~beta`]*Physics:-Christoffel[`~upsilon`, tau, upsilon]+(1/4)*Physics:-Christoffel[`~beta`, nu, sigma]*Physics:-Christoffel[`~sigma`, `~alpha`, `~nu`]-(1/4)*Physics:-Christoffel[`~beta`, sigma, `~nu`]*Physics:-Christoffel[`~sigma`, nu, `~alpha`]-(1/4)*Physics:-Christoffel[`~alpha`, lambda, `~mu`]*Physics:-Christoffel[`~lambda`, mu, `~beta`]+(1/4)*Physics:-Christoffel[`~alpha`, lambda, mu]*Physics:-Christoffel[`~lambda`, `~beta`, `~mu`]-(1/2)*Physics:-D_[`~beta`](Physics:-Christoffel[`~upsilon`, upsilon, `~alpha`], [X])+(1/4)*Physics:-D_[`~nu`](Physics:-Christoffel[`~beta`, nu, `~alpha`], [X])+(1/4)*Physics:-D_[`~mu`](Physics:-Christoffel[`~alpha`, mu, `~beta`], [X])-(1/4)*Physics:-D_[`~omega`](Physics:-Christoffel[`~beta`, omega, `~alpha`], [X])+(1/2)*Physics:-D_[`~beta`](Physics:-Christoffel[alpha9, `~alpha`, `~alpha9`], [X])-(1/4)*Physics:-D_[`~rho1`](Physics:-Christoffel[`~alpha`, rho1, `~beta`], [X]) = 0

(29)

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

Simplify(-(1/2)*D_[`~beta`](Physics[Christoffel][`~upsilon`, upsilon, `~alpha`], [X])+(1/4)*D_[`~nu`](Physics[Christoffel][`~beta`, nu, `~alpha`], [X])+(1/4)*D_[`~mu`](Physics[Christoffel][`~alpha`, mu, `~beta`], [X])-(1/4)*D_[`~omega`](Physics[Christoffel][`~beta`, omega, `~alpha`], [X])+(1/2)*D_[`~beta`](Physics[Christoffel][alpha9, `~alpha`, `~alpha9`], [X])-(1/4)*D_[`~rho1`](Physics[Christoffel][`~alpha`, rho1, `~beta`], [X])-Physics[Ricci][`~alpha`, `~beta`]-(1/2)*Physics[Christoffel][alpha9, alpha10, `~alpha`]*Physics[Christoffel][`~alpha10`, `~alpha9`, `~beta`]-(1/4)*Physics[Christoffel][`~alpha`, chi, iota]*Physics[Christoffel][`~chi`, `~beta`, `~iota`]-(1/4)*Physics[Christoffel][`~beta`, chi, iota]*Physics[Christoffel][`~chi`, `~alpha`, `~iota`]+(1/2)*Physics[Christoffel][chi, `~alpha`, `~beta`]*Physics[Christoffel][iota, `~chi`, `~iota`]+(1/4)*(Physics[Christoffel][`~alpha`, chi, `~beta`]+Physics[Christoffel][`~beta`, chi, `~alpha`])*Physics[Christoffel][`~chi`, iota, `~iota`]-(1/2)*Physics[Christoffel][chi, iota, `~alpha`]*Physics[Christoffel][`~iota`, `~beta`, `~chi`]-(1/4)*Physics[Christoffel][`~alpha`, rho, `~beta`]*Physics[Christoffel][`~rho`, rho1, `~rho1`]+(1/4)*Physics[Christoffel][`~alpha`, rho, `~rho1`]*Physics[Christoffel][`~rho`, rho1, `~beta`]-(1/2)*Physics[Christoffel][omicron, zeta, `~omicron`]*Physics[Christoffel][`~zeta`, `~alpha`, `~beta`]+(1/2)*Physics[Christoffel][omicron, zeta, `~beta`]*Physics[Christoffel][`~zeta`, `~alpha`, `~omicron`]-(1/4)*Physics[Christoffel][`~beta`, psi, `~alpha`]*Physics[Christoffel][`~psi`, omega, `~omega`]+(1/4)*Physics[Christoffel][`~beta`, psi, `~omega`]*Physics[Christoffel][`~psi`, omega, `~alpha`]+(1/2)*Physics[Christoffel][`~tau`, upsilon, `~alpha`]*Physics[Christoffel][`~upsilon`, tau, `~beta`]-(1/2)*Physics[Christoffel][`~tau`, `~alpha`, `~beta`]*Physics[Christoffel][`~upsilon`, tau, upsilon]+(1/4)*Physics[Christoffel][`~beta`, nu, sigma]*Physics[Christoffel][`~sigma`, `~alpha`, `~nu`]-(1/4)*Physics[Christoffel][`~beta`, sigma, `~nu`]*Physics[Christoffel][`~sigma`, nu, `~alpha`]-(1/4)*Physics[Christoffel][`~alpha`, lambda, `~mu`]*Physics[Christoffel][`~lambda`, mu, `~beta`]+(1/4)*Physics[Christoffel][`~alpha`, lambda, mu]*Physics[Christoffel][`~lambda`, `~beta`, `~mu`]+(1/4)*(2*Physics[Christoffel][chi, iota, kappa]*Physics[Christoffel][`~iota`, `~chi`, `~kappa`]-2*Physics[Christoffel][chi, iota, `~chi`]*Physics[Christoffel][`~iota`, kappa, `~kappa`]-2*Physics[Christoffel][alpha4, alpha5, alpha6]*Physics[Christoffel][`~alpha6`, `~alpha4`, `~alpha5`]+2*Physics[Christoffel][alpha4, alpha6, `~alpha5`]*Physics[Christoffel][`~alpha6`, alpha5, `~alpha4`]-2*D_[`~alpha5`](Physics[Christoffel][alpha4, alpha5, `~alpha4`], [X])+2*Physics[Christoffel][`~alpha1`, alpha1, alpha3]*Physics[Christoffel][`~alpha3`, alpha2, `~alpha2`]-2*Physics[Christoffel][`~alpha1`, alpha3, `~alpha2`]*Physics[Christoffel][`~alpha3`, alpha1, alpha2]+2*Physics[Ricci][alpha2, `~alpha2`]+2*D_[`~alpha2`](Physics[Christoffel][`~alpha1`, alpha1, alpha2], [X]))*Physics[g_][`~alpha`, `~beta`]+(1/2)*Physics[Christoffel][alpha10, `~alpha`, `~beta`]*Physics[Christoffel][alpha9, `~alpha10`, `~alpha9`] = 0)

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

(30)

NULL


 

Download Einstein_Equations_From_First_Principles.mw

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

I’m absolutely delighted to announce the launch of Maple 2025!

Although you see a new release every year, new features take anything from a few fast-paced weeks to develop, to months of careful cultivation.

Working on so many features in parallel, each with varying time scales, isn't easy! We have to fastidiously manage and track our work.

So it's easy to lose ourselves in the daily minutiae of software development. To help us maintain perspective, we constantly ask ourselves questions like:

  • What user problem are we solving and how often does this problem occur?
  • Can we validate our proposed solution with preliminary user feedback?
  • Is this a solution to a problem that doesn't exist and will never exist, or are we pre-empting a future need?
  • Are we offering value to our users?

Given the answers, we course-correct to make sure we stay on track for our central mission - to make you happy, and to keep you coming back year-after-year.

With Maple 2025, I think we've smashed that goal. We have many new features that'll appeal to many different types of users - from students, educators and mathematicians, to engineers, scientists and technical professionals

Let me walk you through some of my personal highlights.

It’ll be difficult for anyone to miss this - Maple 2025 has a new interface! It’s a ribbon-based UI that look clean and contemporary, and helps you find and discover tools more quickly than before.

You have large, meaningful icons.

Items are logically grouped.

The ribbons is contextual. If you click on a plot, you get new tabs for interacting with and drawing on the plot.

A new Education tab collects pedagogical resources that were scattered around the interface in prior releases.

This is the biggest visual overhaul to Maple in many years. We hope you like it! 

We also appreciate that changes in look and feel can be divisive. Please rest assured that we will refine and finesse the interface with each successive release; your comments and suggestions are most welcome.

The new interface is available on Windows and Linux, and as a technology preview on Mac.

The right arrow key on my keyboard is wearing out…and it’s all because of Maple. I’m knee deep in Maple nearly every day entering equations, and I’m always using right-arrow to move the cursor. It gets kind of tedious!

This anecdote reflects some investigative work we did. We comprehensively examined our internal library of thousands of Maple worksheets and discovered that these three input patterns are extremely common.

Previously, you’d use the right-arrow key to move the cursor out of the exponential, division or subscript.

Now, in Maple 2025, when you

  • type ^, /, or enter a literal subscript with a double-underscore,
  • followed by a number or symbol
  • and then input another operator (such as +)

the operator is automatically inserted on the baseline (except when y = 1).

Of course, you can also make the cursor return to or stay in the exponent or denominator with a simple keystroke, when that is what is needed.

This is one of those little quality of life refinements that I’m very fond of - it’s a little visual and usability dopamine hit.

The sum command (and its typeset form) now indexes into vectors without you needing to spam unevaluation quotes all over your expression.

Maple 2024

Maple 2025

We’ve been integrating units deeper into the Maple system, release after release. Much of this is driven by our engineering users.

A few releases ago, we made int(numeric) compatible with units. With Maple 2025, you can now numerically differentiate  expressions and procedures that have units.

I’m a grizzled thermodynamics hack, so here’s an example in which I calculate the specific heat capacity of water by differentiating enthalpy with respect to temperature (and then confirm the result with the built-in value):

This is in addition to many other improvements to the units experience.

Although this is a part of Maple that I don’t touch often (my colleague Karishma takes point on the education side), I REALLY wish I’d had this when I was struggling with math.

You can now automatically generate unlimited variants of the same problem for students to solve with the Try Another feature, which has been added to Maple’s Check My Work tools (another feature I really could have used!). This is available for many common math principles, including factorization, simplification, integration and more.

This is just one of the improvements in Maple 2025 for teaching and learning.

 If you’ve ever found yourself going back and forth (and back and forth) between two large, almost identical-looking Maple expressions, trying to figure out how they are different, you’re going to love this one.  ExpressionTools is a new package that lets you compare the differences between two expressions.

I really like the use of color to highlight differences. Less squinting at the screen!

You can now run Maple Flow worksheets from Maple (you don’t need Maple Flow installed to do this). You can send parameters into the Flow worksheet and extract your desired results.

This means you can use the entire flexibility of Maple to analysis and manipulate your Flow worksheet. You could, for example:

  • Attach a Flow worksheet to a Maple workbook and create an interactive application
  • Carry out parameter studies of a Flow worksheet by evaluating it over many parameter sets in Maple
  • Create an Excel interface for a Flow worksheet using the Maple add-in for Excel

Simplify is one of those functions that literally tens of thousands of people use each day. Every time we make an incremental improvement, the cumulative benefits across our entire user base are significant.

We’ve refined simplify in a number of critical ways. For example, simplify now recognizes when exponentials can be profitably converted to hyperbolic trig functions:

The analysis of many scientific phenomena result in Laplace transforms that do not have a symbolic inverse which can be expressed in terms of elementary functions. This includes applications in heat transfer, fluid mechanics, fractional diffusion processes, control systems and electrical transmission.

For example, this monster Laplace transform results from an analysis of voltage on a transmission line:

You can now numerically invert this transform courtesy of an enhancement to inttrans:-invlaplace - a fast quadrature method.

I’ve saved what I think has the most future potential for last.

I’m sure nearly all of you have experimented with the various AI tools. They’re an inevitable part of our present and future, whether we're comfortable with it or not.

This is something we've been mulling over for some time.

  • In Maple 2019, the DeepLearning package made its debut. This package provides tools for machine learning, supporting operations such as classification and regression using neural networks.
  • In Maple 2024, we introduced an AI-powered formula lookup feature.

In Maple 2025, we’re giving you an early-stage technology preview of AI-powered document generation.

You can automatically generate worksheet content by prompting an AI, and then gradually refine the content

If you’re an educator, you might want some content that describes applications of calculus. So you might ask the AI “How do I derive the formula for the area of a circle” by entering your prompt into this text box:

This is the worksheet content that may be returned:

If you’re structural engineer who wants to know how to calculate the hardness of concrete, you might ask the AI: “How do I calculate the compressive strength of slow hardening concrete as a function of time? Use the CEB-FIP Model Code 90. Include a worked example with Maple code”.

This worksheet content that could be generated (note the live Maple code):

We’re labelling AI-generated worksheet content as a technology preview. You might see

  • text that might be misleading (but sounds plausible)
  • code that doesn’t work (but looks plausible)
  • or different results each time you click “Generate Document”

For the moment, I would not rely on AI-generated worksheet content without realistic expectations, a healthy dose of scepticism and a modicum of detached analysis. But AI models are rapidly growing in robustness, and we want to position ourselves to best exploit their future potential. The next few years will be VERY exciting.

We can never cover everything in a short blog post like this. So if you want to know more, head on over to the What’s New pages for Maple 2025!

Hi MaplePrimes,
I made a quick code to calculate Aliquot Sequences.

The code works.

I was inspired by a Numberphile YouTube video.

see

aliquot_sequence_cut_1.mw

aliquot_sequence_cut_1.pdf

This should be helpful.

Regards,

Matt

1 2 3 4 5 6 7 Last Page 1 of 76