Items tagged with model

Feed

 

Physics

 

 

Maple provides a state-of-the-art environment for algebraic and tensorial computations in Physics, with emphasis on ensuring that the computational experience is as natural as possible.

 

The theme of the Physics project for Maple 2017 has been the consolidation of the functionality introduced in previous releases, together with significant enhancements and new functionality in General Relativity, in connection with classification of solutions to Einstein's equations and tensor representations to work in an embedded 3D curved space - a new ThreePlusOne  package. This package is relevant in numerical relativity and a Hamiltonian formulation of gravity. The developments also include first steps in connection with computational representations for all the objects entering the Standard Model in particle physics.

Classification of solutions to Einstein's equations and the Tetrads package

 

In Maple 2016, the digitizing of the database of solutions to Einstein's equations  was finished, added to the standard Maple library, with all the metrics from "Stephani, H.; Kramer, D.; MacCallum, M.; Hoenselaers, C.; and Herlt, E., Exact Solutions to Einstein's Field Equations". These metrics can be loaded to work with them, or change them, or searched using g_  (the Physics command representing the spacetime metric that also sets the metric to your choice in one go) or using the command DifferentialGeometry:-Library:-MetricSearch .


In Maple 2017, the Physics:-Tetrads  package has been vastly improved and extended, now including new commands like PetrovType  and SegreType  to classify these metrics, and the TransformTetrad  now has an option canonicalform to automatically derive a transformation and put the tetrad in canonical form (reorientation of the axis of the local system of references), a relevant step in resolving the equivalence between two metrics.

Examples

 

Petrov and Segre types, tetrads in canonical form

   

Equivalence for Schwarzschild metric (spherical and Kruskal coordinates)

 

Formulation of the problem (remove mixed coordinates)

   

Solving the Equivalence

   

The ThreePlusOne (3 + 1) new Maple 2017 Physics package

 

ThreePlusOne , is a package to cast Einstein's equations in a 3+1 form, that is, representing spacetime as a stack of nonintersecting 3-hypersurfaces Σ. This 3+1 description is key in the Hamiltonian formulation of gravity as well as in the study of gravitational waves, black holes, neutron stars, and in general to study the evolution of physical system in general relativity by running numerical simulations as traditional initial value (Cauchy) problems. ThreePlusOne includes computational representations for the spatial metric gamma[i, j] that is induced by g[mu, nu] on the 3-dimensional hypersurfaces, and the related covariant derivative, Christoffel symbols and Ricci and Riemann tensors, the Lapse, Shift, Unit normal and Time vectors and Extrinsic curvature related to the ADM equations.

 

The following is a list of the available commands:

 

ADMEquations

Christoffel3

D3_

ExtrinsicCurvature

gamma3_

Lapse

Ricci3

Riemann3

Shift

TimeVector

UnitNormalVector

 

 

The other four related new Physics  commands:

 

• 

Decompose , to decompose 4D tensorial expressions (free and/or contracted indices) into the space and time parts.

• 

gamma_ , representing the three-dimensional metric tensor, with which the element of spatial distance is defined as  `#mrow(msup(mi("dl"),mrow(mo("⁢"),mn("2"))),mo("="),msub(mi("γ",fontstyle = "normal"),mrow(mi("i"),mo(","),mi("j"))),mo("⁢"),msup(mi("dx"),mi("i")),mo("⁢"),msup(mi("dx"),mi("j")))`.

• 

Redefine , to redefine the coordinates and the spacetime metric according to changes in the signature from any of the four possible signatures(− + + +), (+ − − −), (+ + + −) and ((− + + +) to any of the other ones.

• 

EnergyMomentum , is a computational representation for the energy-momentum tensor entering Einstein's equations as well as their 3+1 form, the ADMEquations .

 

Examples

 

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

`Default differentiation variables for d_, D_ and dAlembertian are: `*{X = (x, y, z, t)}

 

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

 

[coordinatesystems = {X}]

(2.1.1)

with(ThreePlusOne)

`Setting lowercaselatin_is letters to represent space indices `

 

0, "%1 is not a command in the %2 package", ThreePlusOne, Physics

 

`Changing the signature of spacetime from `(`- - - +`)*` to `(`+ + + -`)*` in order to match the signature customarily used in the ADM formalism`

 

[ADMEquations, Christoffel3, D3_, ExtrinsicCurvature, Lapse, Ricci3, Riemann3, Shift, TimeVector, UnitNormalVector, gamma3_]

(2.1.2)

Note the different color for gamma[mu, nu], now a 4D tensor representing the metric of a generic 3-dimensional hypersurface induced by the 4D spacetime metric g[mu, nu]. All the ThreePlusOne tensors are displayed in black to distinguish them of the corresponding 4D or 3D tensors. The particular hypersurface gamma[mu, nu] operates is parameterized by the Lapse  alpha and the Shift  beta[mu].

The induced metric gamma[mu, nu]is defined in terms of the UnitNormalVector  n[mu] and the 4D metric g[mu, nu] as

gamma3_[definition]

Physics:-ThreePlusOne:-gamma3_[mu, nu] = Physics:-ThreePlusOne:-UnitNormalVector[mu]*Physics:-ThreePlusOne:-UnitNormalVector[nu]+Physics:-g_[mu, nu]

(2.1.3)

where n[mu] is defined in terms of the Lapse  alpha and the derivative of a scalar function t that can be interpreted as a global time function

UnitNormalVector[definition]

Physics:-ThreePlusOne:-UnitNormalVector[mu] = -Physics:-ThreePlusOne:-Lapse*Physics:-D_[mu](t)

(2.1.4)

The TimeVector  is defined in terms of the Lapse  alpha and the Shift  beta[mu] and this vector  n[mu] as

TimeVector[definition]

Physics:-ThreePlusOne:-TimeVector[mu] = Physics:-ThreePlusOne:-Lapse*Physics:-ThreePlusOne:-UnitNormalVector[mu]+Physics:-ThreePlusOne:-Shift[mu]

(2.1.5)

The ExtrinsicCurvature  is defined in terms of the LieDerivative  of  gamma[mu, nu]

ExtrinsicCurvature[definition]

Physics:-ThreePlusOne:-ExtrinsicCurvature[mu, nu] = -(1/2)*Physics:-LieDerivative[Physics:-ThreePlusOne:-UnitNormalVector](Physics:-ThreePlusOne:-gamma3_[mu, nu])

(2.1.6)

The metric gamma[mu, nu]is also a projection tensor in that it projects 4D tensors into the 3D hypersurface Σ. The definition for any 4D tensor that is also a 3D tensor in Σ, can thus be written directly by contracting their indices with gamma[mu, nu]. In the case of Christoffel3 , Ricci3  and Riemann3,  these tensors can be defined by replacing the 4D metric g[mu, nu] by gamma[mu, nu] and the 4D Christoffel symbols GAMMA[mu, nu, alpha] by the ThreePlusOne GAMMA[mu, nu, alpha] in the definitions of the corresponding 4D tensors. So, for instance

Christoffel3[definition]

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

(2.1.7)

Ricci3[definition]

Physics:-ThreePlusOne:-Ricci3[mu, nu] = Physics:-d_[alpha](Physics:-ThreePlusOne:-Christoffel3[`~alpha`, mu, nu], [X])-Physics:-d_[nu](Physics:-ThreePlusOne:-Christoffel3[`~alpha`, mu, alpha], [X])+Physics:-ThreePlusOne:-Christoffel3[`~beta`, mu, nu]*Physics:-ThreePlusOne:-Christoffel3[`~alpha`, beta, alpha]-Physics:-ThreePlusOne:-Christoffel3[`~beta`, mu, alpha]*Physics:-ThreePlusOne:-Christoffel3[`~alpha`, nu, beta]

(2.1.8)

Riemann3[definition]

Physics:-ThreePlusOne:-Riemann3[mu, nu, alpha, beta] = Physics:-g_[mu, lambda]*(Physics:-d_[alpha](Physics:-ThreePlusOne:-Christoffel3[`~lambda`, nu, beta], [X])-Physics:-d_[beta](Physics:-ThreePlusOne:-Christoffel3[`~lambda`, nu, alpha], [X])+Physics:-ThreePlusOne:-Christoffel3[`~lambda`, upsilon, alpha]*Physics:-ThreePlusOne:-Christoffel3[`~upsilon`, nu, beta]-Physics:-ThreePlusOne:-Christoffel3[`~lambda`, upsilon, beta]*Physics:-ThreePlusOne:-Christoffel3[`~upsilon`, nu, alpha])

(2.1.9)

When working with the ADM formalism, the line element of an arbitrary spacetime metric can be expressed in terms of the differentials of the coordinates dx^mu, the Lapse , the Shift  and the spatial components of the 3D metric gamma3_ . From this line element one can derive the relation between the Lapse , the spatial part of the Shift , the spatial part of the gamma3_  metric and the g[0, j] components of the 4D spacetime metric.

For this purpose, define a tensor representing the differentials of the coordinates and an alias  dt = `#msup(mi("dx"),mn("0"))`

Define(dx[mu])

`Defined objects with tensor properties`

 

{Physics:-ThreePlusOne:-D3_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-ThreePlusOne:-Ricci3[mu, nu], Physics:-ThreePlusOne:-Shift[mu], Physics:-d_[mu], dx[mu], Physics:-g_[mu, nu], Physics:-ThreePlusOne:-gamma3_[mu, nu], Physics:-gamma_[i, j], Physics:-ThreePlusOne:-Christoffel3[mu, nu, alpha], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-ThreePlusOne:-Riemann3[mu, nu, alpha, beta], Physics:-ThreePlusOne:-TimeVector[mu], Physics:-ThreePlusOne:-ExtrinsicCurvature[mu, nu], Physics:-ThreePlusOne:-UnitNormalVector[mu], Physics:-SpaceTimeVector[mu](X)}

(2.1.10)

"alias(dt = dx[~0]):"

The expression for the line element in terms of the Lapse  and Shift   is (see [2], eq.(2.123))

ds^2 = (-Lapse^2+Shift[i]^2)*dt^2+2*Shift[i]*dt*dx[`~i`]+gamma_[i, j]*dx[`~i`]*dx[`~j`]

ds^2 = (-Physics:-ThreePlusOne:-Lapse^2+Physics:-ThreePlusOne:-Shift[i]*Physics:-ThreePlusOne:-Shift[`~i`])*dt^2+2*Physics:-ThreePlusOne:-Shift[i]*dt*dx[`~i`]+Physics:-gamma_[i, j]*dx[`~i`]*dx[`~j`]

(2.1.11)

Compare this expression with the 3+1 decomposition of the line element in an arbitrary system. To avoid the automatic evaluation of the metric components, work with the inert form of the metric %g_

ds^2 = %g_[mu, nu]*dx[`~mu`]*dx[`~nu`]

ds^2 = %g_[mu, nu]*dx[`~mu`]*dx[`~nu`]

(2.1.12)

Decompose(ds^2 = %g_[mu, nu]*dx[`~mu`]*dx[`~nu`])

ds^2 = %g_[0, 0]*dt^2+%g_[0, j]*dt*dx[`~j`]+%g_[i, 0]*dt*dx[`~i`]+%g_[i, j]*dx[`~i`]*dx[`~j`]

(2.1.13)

The second and third terms on the right-hand side are equal

op(2, rhs(ds^2 = dt^2*%g_[0, 0]+dt*%g_[0, j]*dx[`~j`]+dt*%g_[i, 0]*dx[`~i`]+%g_[i, j]*dx[`~i`]*dx[`~j`])) = op(3, rhs(ds^2 = dt^2*%g_[0, 0]+dt*%g_[0, j]*dx[`~j`]+dt*%g_[i, 0]*dx[`~i`]+%g_[i, j]*dx[`~i`]*dx[`~j`]))

%g_[0, j]*dt*dx[`~j`] = %g_[i, 0]*dt*dx[`~i`]

(2.1.14)

subs(%g_[0, j]*dt*dx[`~j`] = %g_[i, 0]*dt*dx[`~i`], ds^2 = dt^2*%g_[0, 0]+dt*%g_[0, j]*dx[`~j`]+dt*%g_[i, 0]*dx[`~i`]+%g_[i, j]*dx[`~i`]*dx[`~j`])

ds^2 = %g_[0, 0]*dt^2+2*%g_[i, 0]*dt*dx[`~i`]+%g_[i, j]*dx[`~i`]*dx[`~j`]

(2.1.15)

Taking the difference between this expression and the one in terms of the Lapse  and Shift  we get

simplify((ds^2 = dt^2*%g_[0, 0]+2*dt*%g_[i, 0]*dx[`~i`]+%g_[i, j]*dx[`~i`]*dx[`~j`])-(ds^2 = (-Physics:-ThreePlusOne:-Lapse^2+Physics:-ThreePlusOne:-Shift[i]*Physics:-ThreePlusOne:-Shift[`~i`])*dt^2+2*Physics:-ThreePlusOne:-Shift[i]*dt*dx[`~i`]+Physics:-gamma_[i, j]*dx[`~i`]*dx[`~j`]))

0 = (Physics:-ThreePlusOne:-Lapse^2-Physics:-ThreePlusOne:-Shift[i]*Physics:-ThreePlusOne:-Shift[`~i`]+%g_[0, 0])*dt^2+2*dx[`~i`]*(%g_[i, 0]-Physics:-ThreePlusOne:-Shift[i])*dt-dx[`~i`]*dx[`~j`]*(Physics:-gamma_[i, j]-%g_[i, j])

(2.1.16)

Taking coefficients, we get equations for the Shift , the Lapse  and the spatial components of the metric gamma3_

eq[1] := coeff(coeff(rhs(0 = (Physics:-ThreePlusOne:-Lapse^2-Physics:-ThreePlusOne:-Shift[i]*Physics:-ThreePlusOne:-Shift[`~i`]+%g_[0, 0])*dt^2+2*dx[`~i`]*(%g_[i, 0]-Physics:-ThreePlusOne:-Shift[i])*dt-dx[`~i`]*dx[`~j`]*(Physics:-gamma_[i, j]-%g_[i, j])), dt), dx[`~i`]) = 0

2*%g_[i, 0]-2*Physics:-ThreePlusOne:-Shift[i] = 0

(2.1.17)

eq[2] := coeff(rhs(0 = (Physics:-ThreePlusOne:-Lapse^2-Physics:-ThreePlusOne:-Shift[i]*Physics:-ThreePlusOne:-Shift[`~i`]+%g_[0, 0])*dt^2+2*dx[`~i`]*(%g_[i, 0]-Physics:-ThreePlusOne:-Shift[i])*dt-dx[`~i`]*dx[`~j`]*(Physics:-gamma_[i, j]-%g_[i, j])), dt^2)

Physics:-ThreePlusOne:-Lapse^2-Physics:-ThreePlusOne:-Shift[i]*Physics:-ThreePlusOne:-Shift[`~i`]+%g_[0, 0]

(2.1.18)

eq[3] := coeff(coeff(rhs(0 = (Physics:-ThreePlusOne:-Lapse^2-Physics:-ThreePlusOne:-Shift[i]*Physics:-ThreePlusOne:-Shift[`~i`]+%g_[0, 0])*dt^2+2*dx[`~i`]*(%g_[i, 0]-Physics:-ThreePlusOne:-Shift[i])*dt-dx[`~i`]*dx[`~j`]*(Physics:-gamma_[i, j]-%g_[i, j])), dx[`~i`]), dx[`~j`]) = 0

-Physics:-gamma_[i, j]+%g_[i, j] = 0

(2.1.19)

Using these equations, these quantities can all be expressed in terms of the time and space components of the 4D metric g[0, 0] and g[i, j]

isolate(eq[1], Shift[i])

Physics:-ThreePlusOne:-Shift[i] = %g_[i, 0]

(2.1.20)

isolate(eq[2], Lapse^2)

Physics:-ThreePlusOne:-Lapse^2 = Physics:-ThreePlusOne:-Shift[i]*Physics:-ThreePlusOne:-Shift[`~i`]-%g_[0, 0]

(2.1.21)

isolate(eq[3], gamma_[i, j])

Physics:-gamma_[i, j] = %g_[i, j]

(2.1.22)

References

 
  

[1] Landau, L.D., and Lifshitz, E.M. The Classical Theory of Fields, Course of Theoretical Physics Volume 2, fourth revised English edition. Elsevier, 1975.

  

[2] Alcubierre, M., Introduction to 3+1 Numerical Relativity, International Series of Monographs on Physics 140, Oxford University Press, 2008.

  

[3] Baumgarte, T.W., Shapiro, S.L., Numerical Relativity, Solving Einstein's Equations on a Computer, Cambridge University Press, 2010.

  

[4] Gourgoulhon, E., 3+1 Formalism and Bases of Numerical Relativity, Lecture notes, 2007, https://arxiv.org/pdf/gr-qc/0703035v1.pdf.

  

[5] Arnowitt, R., Dese, S., Misner, C.W., The Dynamics of General Relativity, Chapter 7 in Gravitation: an introduction to current research (Wiley, 1962), https://arxiv.org/pdf/gr-qc/0405109v1.pdf

  

 

Examples: Decompose, gamma_

 

restartwith(Physics)NULL

Setup(mathematicalnotation = true)

[mathematicalnotation = true]

(2.2.1)

Define  now an arbitrary tensor A

Define(A)

`Defined objects with tensor properties`

 

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

(2.2.2)

So A^mu is a 4D tensor with only one free index, where the position of the time-like component is the position of the different sign in the signature, that you can query about via

Setup(signature)

[signature = `- - - +`]

(2.2.3)

To perform a decomposition into space and time, set - for instance - the lowercase latin letters from i to s to represent spaceindices and

Setup(spaceindices = lowercase_is)

[spaceindices = lowercaselatin_is]

(2.2.4)

Accordingly, the 3+1 decomposition of A^mu is

Decompose(A[`~mu`])

Array(%id = 18446744078724512334)

(2.2.5)

The 3+1 decomposition of the inert representation %g_[mu,nu] of the 4D spacetime metric; use the inert representation when you do not want the actual components of the metric appearing in the output

Decompose(%g_[mu, nu])

Matrix(%id = 18446744078724507998)

(2.2.6)

Note the position of the component %g_[0, 0], related to the trailing position of the time-like component in the signature "(- - - +)".

Compare the decomposition of the 4D inert with the decomposition of the 4D active spacetime metric

g[]

g[mu, nu] = (Matrix(4, 4, {(1, 1) = -1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(2.2.7)

Decompose(g_[mu, nu])

Matrix(%id = 18446744078724494270)

(2.2.8)

Note that in general the 3D space part of g[mu, nu] is not equal to the 3D metric gamma[i, j] whose definition includes another term (see [1] Landau & Lifshitz, eq.(84.7)).

gamma_[definition]

Physics:-gamma_[i, j] = -Physics:-g_[i, j]+Physics:-g_[0, i]*Physics:-g_[0, j]/Physics:-g_[0, 0]

(2.2.9)

The 3D space part of -g[`~mu`, `~nu`] is actually equal to the 3D metric "gamma[]^(i,j)"

"gamma_[~,definition];"

Physics:-gamma_[`~i`, `~j`] = -Physics:-g_[`~i`, `~j`]

(2.2.10)

To derive the formula  for the covariant components of the 3D metric, Decompose into 3+1 the identity

%g_[`~alpha`, `~mu`]*%g_[mu, beta] = KroneckerDelta[`~alpha`, beta]

%g_[`~alpha`, `~mu`]*%g_[mu, beta] = Physics:-KroneckerDelta[beta, `~alpha`]

(2.2.11)

To the side, for illustration purposes, these are the 3 + 1 decompositions, first excluding the repeated indices, then excluding the free indices

Eq := Decompose(%g_[`~alpha`, `~mu`]*%g_[mu, beta] = Physics[KroneckerDelta][beta, `~alpha`], repeatedindices = false)

Matrix(%id = 18446744078132963318)

(2.2.12)

Eq := Decompose(%g_[`~alpha`, `~mu`]*%g_[mu, beta] = Physics[KroneckerDelta][beta, `~alpha`], freeindices = false)

%g_[0, beta]*%g_[`~alpha`, `~0`]+%g_[i, beta]*%g_[`~alpha`, `~i`] = Physics:-KroneckerDelta[beta, `~alpha`]

(2.2.13)

Compare with a full decomposition

Eq := Decompose(%g_[`~alpha`, `~mu`]*%g_[mu, beta] = Physics[KroneckerDelta][beta, `~alpha`])

Matrix(%id = 18446744078724489454)

(2.2.14)

Eq is a symmetric matrix of equations involving non-contracted occurrences of `#msup(mi("g"),mrow(mn("0"),mo(","),mn("0")))`, `#msup(mi("g"),mrow(mi("j"),mo(","),mo("0")))` and `#msup(mi("g"),mrow(mi("j"),mo(","),mi("i")))`. Isolate, in Eq[1, 2], `#msup(mi("g"),mrow(mi("j"),mo(","),mo("0")))`, that you input as %g_[~j, ~0], and substitute into Eq[1, 1]

"isolate(Eq[1, 2], `%g_`[~j, ~0]);"

%g_[`~j`, `~0`] = -%g_[i, 0]*%g_[`~j`, `~i`]/%g_[0, 0]

(2.2.15)

subs(%g_[`~j`, `~0`] = -%g_[i, 0]*%g_[`~j`, `~i`]/%g_[0, 0], Eq[1, 1])

-%g_[0, k]*%g_[i, 0]*%g_[`~j`, `~i`]/%g_[0, 0]+%g_[i, k]*%g_[`~j`, `~i`] = Physics:-KroneckerDelta[k, `~j`]

(2.2.16)

Collect `#msup(mi("g"),mrow(mi("j"),mo(","),mi("i")))`, that you input as %g_[~j, ~i]

collect(-%g_[0, k]*%g_[i, 0]*%g_[`~j`, `~i`]/%g_[0, 0]+%g_[i, k]*%g_[`~j`, `~i`] = Physics[KroneckerDelta][k, `~j`], %g_[`~j`, `~i`])

(-%g_[0, k]*%g_[i, 0]/%g_[0, 0]+%g_[i, k])*%g_[`~j`, `~i`] = Physics:-KroneckerDelta[k, `~j`]

(2.2.17)

Since the right-hand side is the identity matrix and, from , `#msup(mi("g"),mrow(mi("i"),mo(","),mi("j")))` = -`#msup(mi("γ",fontstyle = "normal"),mrow(mi("i"),mo(","),mi("j")))`, the expression between parenthesis, multiplied by -1, is the reciprocal of the contravariant 3D metric `#msup(mi("γ",fontstyle = "normal"),mrow(mi("i"),mo(","),mi("j")))`, that is the covariant 3D metric gamma[i, j], in accordance to its definition for the signature `- - - +`

gamma_[definition]

Physics:-gamma_[i, j] = -Physics:-g_[i, j]+Physics:-g_[0, i]*Physics:-g_[0, j]/Physics:-g_[0, 0]

(2.2.18)

NULL

References

 
  

[1] Landau, L.D., and Lifshitz, E.M. The Classical Theory of Fields, Course of Theoretical Physics Volume 2, fourth revised English edition. Elsevier, 1975.

Example: Redefine

   

Tensors in Special and General Relativity

 

A number of relevant changes happened in the tensor routines of the Physics package, towards making the routines pack more functionality both for special and general relativity, as well as working more efficiently and naturally, based on Maple's Physics users' feedback collected during 2016.

New functionality

 
• 

Implement conversions to most of the tensors of general relativity (relevant in connection with functional differentiation)

• 

New setting in the Physics Setup  allows for specifying the cosmologicalconstant and a default tensorsimplifier


 

Download PhysicsMaple2017.mw

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

I modelled a telecamera and line follower sensors in Maple in order to calculate the "error" for my PID controller. Know i would like to share this information with the PID in MapleSim but I have no idea how to do it. Is there a way to import maple libraries into MapleSim?

Thanks!

 

Howdy,

I am trying to create a mathmatical model that shows a predator-prey relationship along with seasonal variations (hibernation). I made two piecewise functions where one is "on" during the winter and the other is "off" during the winter and the model worked well for one unit of time. 

My problem is that I would like the piecewise functions to apply over an arbitrary amount of time units periods without having to set the range within the function itself for 100 units of time. I was thinking that a for loop would work but i'm not sure how I would impliment that.

Thanks in advance!!

Edit: 

Here's a little bit more info on the problem.

I have three O.D.E's in an array..

ode := {r'(t),h'(t),c'(t)}

which are dependent on a list of parameters..

constants := {a,b,...,s,w,c0,h0,r0};

where s and w are my piecewise functions.

I am able to solve the three odes using dsolve for one period of s and w.

Again, I'm really not sure where to start to make my model periodic.

The CIELAB perceptual model of human vision can be used to predict the color of a wavelength in sRGB. I found the CIEDE2000 corrections to hue, chroma, and lighness were required to get the right values in the blue to violet region. The CIEDE2000 results seem pretty close to what I see looking through a diffraction grating considering I have no way of knowing how my eyes are adapted.

Use the new, faster  version. It can run in under 1 minute;

nprocWavelength_Color_CIEDE2000.mw   and  WtoC_data.xlsx

You will need the data for the color matching functions in the WtoC Excel File. The data directory can now be controlled in the mw file.

 

loock i have 3 models and a Matrix
>m3:=importmatrix("loslobos.txt");
>model[1]:=a*x;
>model[2]:=a^2*x+b;
>model[3]:=C-x^2+a;

so i apply a command named DataFit and it gives me 
>g:=Datafit(matrix,model[1]);

and that throws
         g:=[3,[a:=5],33];
so
         >g[1];
3
         >g[2];
5
         >g[3];
33
the variable that i care about is g[1] so every 3 models have differents g[1] so i need to pick the smallest g[1] and then use his model.
Example:=
>g:=Datafit(matrix,model[1]);
g:=[3,[a:=5],33];
>gg:=Datafit(matrix,model[2]);
gg:=[10,[a:=5,b:=10],44];
>G:=Datafit(matrix,model[3]);
G:=[1.1,[C:=5,a:=2],88];

then 

>plots:-display(plot(m3, style = point, symbol = diamond, symbolsize = 9), plot(eval("MODEL THAT SHOULD BE PLOT", THE G[2] OF THE MODEL), x = 0 .. 27, color = black));

The right one here would be 

>plots:-display(plot(m3, style = point, symbol = diamond, symbolsize = 9), plot(eval(model[3], G[2]), x = 0 .. 27, color = black));

so how do i select the model[3] for my plot? i knnow how the eval works so the problem here is to pick the right model please Help i hope that i isn't that Hard thnx
   

The model of fixed-bed adsorption column

Fluid phase:

PDE:= diff(U(x, tau),tau)+ psi*Theta*diff(U(x, tau),x)-(1/Pe)*psi*Theta*diff(U(x, tau),$(x, 2))=-3*psi*xi*(U(x, tau)-Q/K);

 

IBC:={U(x, 0) = 0,U(0, tau) = 1+(1/Pe)*(D[1](U))(0, tau),(D[1](U))(1, tau)=0};

Particle:

PDE:= diff(Q(r, tau), tau) = diff(Q(r, tau), $(r, 2))+(2/r)*diff(Q(r, tau),r);

IBC:={Q(r, 0) = 0,(D[1](Q))(0, tau) = 0,(1/K)*(D[1](Q))(1, tau)=xi*(U-Q(1, tau)/K)};

Pe:=0.01:

psi:=6780:

Theta:=3.0:

xi:=10000:

I will really appreciate your help. Thanks in anticipation.

The following error occurred when I simulate a build-in model, anyone could help me to solve this problem? Thanks first

I want to make the model which moves along the specific direction (translational), So I used the "prismatic joint" to give the direction that I want to enforce. However, "Prismatic joint" only offers the direction along the "X", "Y", and "Z" axes though I want to give other direction.

 

Is there any way to give the specific direction (vector) to make my model move in that way ?

Hi there,

I would like to compute and display the nullclines of a set of ordinary differential equations.

AFAIK, I can compute the nullclines in Maple by defining the equations and solving the system

e.g.:

# Define the equations
eq1 := u(t)*(1-u(t)/kappa)-u(t)*v(t) = 0;
eq2 := g*(u(t)-1)*v(t) = 0;

# Solve the system (i.e. compute the nullclines)
sol := solve({eq1, eq2}, {u(t), v(t)});

However, I am not quite able to imagine how to display them over a dfieldplot or a phaseportrait.

Attached is an example with some differential equations, and their vector field and trajectories: MaplePrimes_Predator_prey_model_nullclines.mw.

It can be use to illustrate how to (compute and) display the nullclines.

 

Thank you,

jon

This application calculates the number of photons reaching a camera sensor for a given exposure. A blackbody model of the sun is generated. The "Sunny 16" rule for exposure is demonstrated. Calculations are done using units.Photon_Exposure_Array.mw

Photon ExposureNULLNULL

Blackbody Model of the Sun

    h := Units:-Standard:-`*`(Units:-Standard:-`*`(0.6626069e-33, Units:-Standard:-`^`(Unit('m'), 2)), Units:-Standard:-`*`(Unit('kg'), Units:-Standard:-`/`(Unit('s')))): 

Plank Constant       

  kb := Units:-Standard:-`*`(Units:-Standard:-`*`(0.1380650e-22, Units:-Standard:-`*`(Units:-Standard:-`^`(Unit('m'), 2), Units:-Standard:-`/`(Units:-Standard:-`^`(Unit('s'), 2)))), Units:-Standard:-`*`(Unit('kg'), Units:-Standard:-`/`(Unit('K')))): 

Boltzman Constant  

c := Units:-Standard:-`*`(0.2997925e9, Units:-Standard:-`*`(Unit('m'), Units:-Standard:-`/`(Unit('s')))):  ``

Light Speed

Rsun := Units:-Standard:-`*`(Units:-Standard:-`*`(6.955, Units:-Standard:-`^`(10, 8)), Unit('m')): ``

Sun Radius  

Re_orb := Units:-Standard:-`*`(Units:-Standard:-`*`(1.496, Units:-Standard:-`^`(10, 11)), Unit('m')): ``

Earth Orbit

Tsun := Units:-Standard:-`*`(5800, Unit('K')): ``

Sun Color Temperature     

 tf_atm := .718: 

Transmission Factor  

 

Sun: Spectral Radiant Exitance to Earth: Spectral Irradiance                   

  "M(lambda):=(2*Pi*h*c^(2))/((lambda)^(5))*1/((e)^((h*c)/(lambda*kb*Tsun))-1)*(Rsun/(Re_orb))^(2)*tf_atm:" NULL

evalf(M(Units:-Standard:-`*`(555, Unit('nm')))) = 1277414308.*Units:-Unit(('kg')/(('m')*('s')^3))"(->)"1.277414308*Units:-Unit(('W')/(('nm')*('m')^2))NULL

Photopic Relative Response VP vs λ

 

csvFile := FileTools[Filename]("/VPhotopic.csv")NULL = "VPhotopic.csv"NULL

VPdata := ImportMatrix(csvFile) = Vector(4, {(1) = ` 471 x 2 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})NULLNULL

 

`λP` := [seq(1 .. 4000)]:

VP := ArrayInterpolation(VPdata, `λP`):             (ArrayInterpolation for x,y data VPdata returns y' for new x data lambdaP)

NULLVParray := [`$`([`λP`[n], VP[n]], n = 1 .. 4000)]:                     

Mearth := [`$`([n, Units:-Standard:-`*`(Units:-Standard:-`*`(M(Units:-Standard:-`*`(n, Unit('nm'))), Unit('nm')), Units:-Standard:-`*`(Units:-Standard:-`^`(Unit('s'), 3), Units:-Standard:-`/`(Unit('kg'))))], n = 1 .. 4000)]:````

``

dualaxisplot(plot([Mearth], lambda = 300 .. 900, style = line, color = [blue], labels = ["λ (nm)", "M (W/nm m^2)"], title = "Spectral Radiant Exitance of the Sun", titlefont = ["ARIAL", 15], legend = [Exitance], size = [800, 300]), plot([VParray], style = line, color = [green], labels = ["λ (nm)", "Relative Response"], legend = [Units:-Standard:-`*`(Units:-Standard:-`*`(Photopic, Relative), Response)]))

 

``

 

 

 

Illuminance in Radiometric and Photometric Units:

E__r := sum(Units:-Standard:-`*`(M(Units:-Standard:-`*`(lambda, Unit('nm'))), Unit('nm')), lambda = 200 .. 4000) = 984.7275549*Units:-Unit(('kg')/('s')^3)"(->)"984.7275549*Units:-Unit(('W')/('m')^2)NULL

NULL

E__po := Units:-Standard:-`*`(Units:-Standard:-`*`(683.002, Units:-Standard:-`*`(Unit('lm'), Units:-Standard:-`/`(Unit('W')))), sum(Units:-Standard:-`*`(Units:-Standard:-`*`(VP[lambda], M(Units:-Standard:-`*`(lambda, Unit('nm')))), Unit('nm')), lambda = 200 .. 4000)) = HFloat(91873.47376063903)*Units:-Unit('lx')NULL

Translation from Illuminance to Luminance for Reflected Light;

 

Object Reflectance          R__o:      

Object Luminance           L__po := proc (R__o) options operator, arrow; R__o*E__po/(Pi*Unit('sr')) end proc:                evalf(L__po(1)) = HFloat(29244.234968360346)*Units:-Unit(('cd')/('m')^2) 

 

Illuminance of a Camera Sensor  Eps applied for time texp determines Luminous Exposure Hp;

Ideal Illuminance is determined by the exposure time texp, effective f-number N and to a less extent the angle to the optical axis θ;

 

• 

H       Luminous Exposure

• 

Eps     Illuminance to the Camera

• 

N                                               Effective F-Number

• 

texp             Exposure Time

• 

θ        Angle to the Optical Axis    

 

E__ps_ideal = Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`/`(4)), L__po), Units:-Standard:-`*`(Units:-Standard:-`^`(cos(theta), 4), Units:-Standard:-`/`(Units:-Standard:-`^`(N, 2)))):

H__p_ideal = Units:-Standard:-`*`(E__ps_ideal, t__exp):

 

The camera meter determines the exposure time texp to balance the object luminance, reflectance and effective f-number. It does this based on an internal constant k and the camera ISO s.

• 

s        ISO Gain (Based on saturation at 3 stops above the average scene luminance)

• 

k       Reflected Light Meter Calibration Constant      k__m := Units:-Standard:-`*`(Units:-Standard:-`*`(12.5, Unit('lx')), Unit('s')):  

                                                                                                  for Nikon, Canon and Sekonic

• 

c        Incident Light Meter Calibration Constant       c__m := Units:-Standard:-`*`(Units:-Standard:-`*`(250, Unit('lx')), Unit('s')):        

                                                                                                  for Sekonic with flat domeNULL

N^2/t__exp = `#mrow(mi("\`E__po\`"),mo("⋅"),mi("s"))`/c__m                        (Incident Light Meter)  NULL 

Units:-Standard:-`*`(Units:-Standard:-`^`(N, 2), Units:-Standard:-`/`(t__exp)) = Units:-Standard:-`*`(`#mrow(mi("\`L__po\`"),mo("⋅"),mi("s"))`, Units:-Standard:-`/`(k__m)):                        (Reflected Light Meter)

NULL

Solve for H in terms of the Camera Meter Constant k and s

 

Es = Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`/`(4)), Lo), Units:-Standard:-`*`(Units:-Standard:-`^`(cos(theta), 4), Units:-Standard:-`/`(Units:-Standard:-`^`(N, 2)))): NULL

t = Units:-Standard:-`*`(Units:-Standard:-`*`(km, Units:-Standard:-`^`(N, 2)), Units:-Standard:-`/`(Units:-Standard:-`*`(Lo, s))):NULL

NULL

NULL

H = Es*t

H = Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`/`(4)), Lo), Units:-Standard:-`*`(Units:-Standard:-`^`(cos(theta), 4), Units:-Standard:-`/`(Units:-Standard:-`^`(N, 2)))), Units:-Standard:-`*`(Units:-Standard:-`*`(km, Units:-Standard:-`^`(N, 2)), Units:-Standard:-`/`(Units:-Standard:-`*`(Lo, s))))"(=)"H = (1/4)*Pi*cos(theta)^4*km/sNULLNULL

 t = H/Es

t = Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`/`(4)), Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`^`(cos(theta), 4), km), Units:-Standard:-`/`(s))), Units:-Standard:-`/`(Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`/`(4)), Lo), Units:-Standard:-`*`(Units:-Standard:-`^`(cos(theta), 4), Units:-Standard:-`/`(Units:-Standard:-`^`(N, 2))))))"(=)"t = km*N^2/(Lo*s)NULLNULL

H__p := proc (s, theta) options operator, arrow; (1/4)*Pi*k__m*cos(theta)^4/s end proc:                                              

  evalf(H__p(100, 0)) = 0.9817477044e-1*Units:-Unit(('cd')*('s')/('m')('radius')^2)"(->)"0.9817477044e-1*Units:-Unit(('lx')*('s'))NULL

 

Note:  Meters are typically set for a scene reflectance 3 stops below 100% or 12.5%.

           

  E__ps := proc (N, R__o, theta) options operator, arrow; (1/4)*Pi*Unit('sr')*R__o*E__po*cos(theta)^4/(Pi*Unit('sr')*N^2) end proc:               

 evalf(E__ps(16, Units:-Standard:-`/`(Units:-Standard:-`^`(2, 3)), 0)) = HFloat(11.215023652421756)*Units:-Unit('lx')                                                                                                   

t__exp_ideal := proc (N, s, R__o) options operator, arrow; H__p(s, theta)/E__ps(N, R__o, theta) end proc:                                     

  evalf(t__exp_ideal(16, 100, Units:-Standard:-`/`(Units:-Standard:-`^`(2, 3)))) = HFloat(0.008753862094289947)*Units:-Unit('s') NULL NULL

 

 

Actual exposure time includes typical lens losses;

 m := Units:-Standard:-`/`(80):``

Magnification  

  T := .9:``

Lens Transmittance

 F := 1.03:``

Lens Flare

V := 1: ``

Vignetting

 

                                                  ``

Total Lens Efficiency

q := Units:-Standard:-`*`(Units:-Standard:-`*`(Units:-Standard:-`*`(T, F), V), Units:-Standard:-`^`(Units:-Standard:-`+`(1, Units:-Standard:-`-`(m)), 2)):                                      evalf(q) = .9039698438NULL

 

Replacing Eps with q*Eps we get the "Sunny 16" relation between exposure time and ISO;  NULL

t__exp := proc (N, s, R__o) options operator, arrow; H__p(s, theta)/(q*E__ps(N, R__o, theta)) end proc:NULL               evalf(t__exp(16, 100, Units:-Standard:-`/`(Units:-Standard:-`^`(2, 3)))) = HFloat(0.009683798806264942)*Units:-Unit('s')NULL

t__exp_alt := proc (N, s, R__o) options operator, arrow; k__m*N^2*Pi/(s*q*R__o*E__po) end proc:                  evalf(t__exp_alt(16, 100, Units:-Standard:-`/`(Units:-Standard:-`^`(2, 3)))) = HFloat(0.00968379880412244)*Units:-Unit('s') 

• 

The Number of Photons NP Reaching the Sensor Area A;

• 

Circle of confusion for 24x36mm "Full Frame" for 1 arcminute view at twice the diagonal:

                          A__cc := Units:-Standard:-`*`(Units:-Standard:-`*`(Pi, Units:-Standard:-`^`(Units:-Standard:-`*`(12.6, Unit('`μm`')), 2)), Units:-Standard:-`/`(4)):    

     

• 

  Sensor Bandwidth                                          Photopic Response VP

• 

  Exposure Time for Zone 5: Rscene=12.5% , Saturation in Zone 8 Rscene=100%

• 

  Camera ISO differs from Saturation ISO. Typical Saturation ISO is 2300 when the camera is set to 3200. See DxoMark.

 

NULL

The average number of photons for exposure time based on Reflectance of the scene  relative to the metered value:    

Zone 5;   R__meter := R__scene: 

NP := proc (s, R__o, theta) options operator, arrow; (1/4)*t__exp(N, s, R__meter)*A__cc*q*R__scene*cos(theta)^4*(sum(VP[lambda]*M(lambda*Unit('nm'))*Unit('nm')*lambda*Unit('nm')/(h*c), lambda = 200 .. 4000))/N^2 end proc: 

                                                                               evalf(NP(2300, 1, Units:-Standard:-`*`(0, Unit('deg')))) = HFloat(2191.5645712603696)  NULL

Zone 8;       R__meter := Units:-Standard:-`*`(R__scene, Units:-Standard:-`/`(Units:-Standard:-`^`(2, 3))):   NULL

NP__sat := proc (s, theta) options operator, arrow; (1/4)*t__exp(N, s, R__meter)*A__cc*q*R__scene*cos(theta)^4*(sum(VP[lambda]*M(lambda*Unit('nm'))*Unit('nm')*lambda*Unit('nm')/(h*c), lambda = 200 .. 4000))/N^2 end proc:  NULL

                                                                              evalf(NP__sat(2300, Units:-Standard:-`*`(0, Unit('deg')))) = HFloat(17532.516570082957)NULL

NULL

 

Approximate Formula

 

H__sat := proc (s__sat) options operator, arrow; H__p(s__sat, 0)*E__ps(N, 1, 0)/E__ps(N, 1/8, 0) end proc:      

                                                                                       evalf(H__sat(s__sat)) = HFloat(78.53981635)*Units:-Unit(('cd')*('s')/('m')('radius')^2)/s__satNULLNULL

Average Visible Photon Energy

P__e_ave := Units:-Standard:-`*`(Units:-Standard:-`/`(Units:-Standard:-`+`(850, -350)), sum(Units:-Standard:-`*`(Units:-Standard:-`*`(h, c), Units:-Standard:-`/`(Units:-Standard:-`*`(lambda, Unit('nm')))), lambda = 350 .. 850)):                    evalf(P__e_ave) = 0.3533174192e-18*Units:-Unit('J') 

NPtyp := proc (s__sat) options operator, arrow; H__sat(s__sat)*A__cc/(683.002*(Unit('lm')/Unit('W'))*P__e_ave) end proc: 

                               evalf(NPtyp(2300)) = HFloat(17644.363333654386)"(->)"HFloat(17644.363333654386)NULL

NULL

 

Download Photon_Exposure_Array.mw

Hello,

I study the model of crankshaft in MapleSim.

For the moment, i manage to plot different parameters function of the time. For exemple, the angle of the crank function of the time or the displacement of the slider function of the time.

But, i like to plot the displacement of the slider funtion of the angle of the crank.

Is it possible and how can i do this?

Thank you for your help.

Hello,

It's a very easy question but i don't find for the moment.

I want to measure an angle in a revolute joint. For that purpose, i use a probe to measure the angle. The angle is measured in rad. How can i do to change in deg ?

Thank you for your help.

Hello,

It's a very easy question but i don't find for the moment.

I want to measure an angle in a revolute joint. For that purpose, i use a probe to measure the angle. The angle is measured in rad. How can i do to change in deg ?

Thank you for your help.

1 2 3 Page 1 of 3