ecterrab

7691 Reputation

20 Badges

14 years, 160 days

MaplePrimes Activity


These are answers submitted by ecterrab

Take a look at the help page of PDEtools, the table of symmetry commands. You can attempt the computation of symmetries of any sort.

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

I just received the following from Katherina von Bulow.

 

Again, the problem:

pde := convert(diff(T(x, y, t), x, x)+diff(T(x, y, t), y, y) = (diff(T(x, y, t), t))/(0.44e-2), rational, exact); conds := T(0, y, t) = 1, T(1, y, t) = 1, T(x, 0, t) = 1, T(x, 1, t) = 1, T(x, y, 0) = 0

diff(diff(T(x, y, t), x), x)+diff(diff(T(x, y, t), y), y) = (2272727273/10000000)*(diff(T(x, y, t), t))

 

T(0, y, t) = 1, T(1, y, t) = 1, T(x, 0, t) = 1, T(x, 1, t) = 1, T(x, y, 0) = 0

(1)

The following change of variables (rather easy ..) suffices to solve the problem as an infinite series, with no composed differential operators around

tr := T(x, y, t) = U(x, y, t)+1

T(x, y, t) = U(x, y, t)+1

(2)

itr := isolate(tr, U(x, y, t))

U(x, y, t) = T(x, y, t)-1

(3)

PDEtools:-dchange(tr, [pde, conds], [U(x, y, t)])

[diff(diff(U(x, y, t), x), x)+diff(diff(U(x, y, t), y), y) = (2272727273/10000000)*(diff(U(x, y, t), t)), U(0, y, t)+1 = 1, U(1, y, t)+1 = 1, U(x, 0, t)+1 = 1, U(x, 1, t)+1 = 1, U(x, y, 0)+1 = 0]

(4)

Call now,

pdsolve([diff(diff(U(x, y, t), x), x)+diff(diff(U(x, y, t), y), y) = (2272727273/10000000)*(diff(U(x, y, t), t)), U(0, y, t)+1 = 1, U(1, y, t)+1 = 1, U(x, 0, t)+1 = 1, U(x, 1, t)+1 = 1, U(x, y, 0)+1 = 0])

U(x, y, t) = Sum(Sum(4*(-(-1)^(n1+n)+(-1)^n1+(-1)^n-1)*sin(n*Pi*x)*sin(n1*Pi*y)*exp(-(10000000/2272727273)*t*Pi^2*(n^2+n1^2))/(n*Pi^2*n1), n = 1 .. infinity), n1 = 1 .. infinity)

(5)

Change variables back and isolate T(x, y, t)

PDEtools:-dchange(itr, U(x, y, t) = Sum(Sum(4*(-(-1)^(n1+n)+(-1)^n1+(-1)^n-1)*sin(n*Pi*x)*sin(n1*Pi*y)*exp(-(10000000/2272727273)*t*Pi^2*(n^2+n1^2))/(n*Pi^2*n1), n = 1 .. infinity), n1 = 1 .. infinity), [T(x, y, t)])

T(x, y, t)-1 = Sum(Sum(4*(-(-1)^(n1+n)+(-1)^n1+(-1)^n-1)*sin(n*Pi*x)*sin(n1*Pi*y)*exp(-(10000000/2272727273)*t*Pi^2*(n^2+n1^2))/(n*Pi^2*n1), n = 1 .. infinity), n1 = 1 .. infinity)

(6)

isolate(T(x, y, t)-1 = Sum(Sum(4*(-(-1)^(n1+n)+(-1)^n1+(-1)^n-1)*sin(n*Pi*x)*sin(n1*Pi*y)*exp(-(10000000/2272727273)*t*Pi^2*(n^2+n1^2))/(n*Pi^2*n1), n = 1 .. infinity), n1 = 1 .. infinity), T(x, y, t))

T(x, y, t) = Sum(Sum(4*(-(-1)^(n1+n)+(-1)^n1+(-1)^n-1)*sin(n*Pi*x)*sin(n1*Pi*y)*exp(-(10000000/2272727273)*t*Pi^2*(n^2+n1^2))/(n*Pi^2*n1), n = 1 .. infinity), n1 = 1 .. infinity)+1

(7)

Besides that this approach solves the problem in a more useful way, the lines above also mean there is something else to adjust in the internal series routines. Although sometimes completely non-obvious, the internal routines are supposed to unveil the need of a transformation of this kind for examples like this one.

``


 

Download PDE_solution_via_change_of_variables.mw

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

Hi Larry
Good catch. The problem is fixed for everybody within the Physics Updates v.453 and higher.

Details: the problem was not related to  Library:-PerformMatrixOperations - you could reproduce it just entering mt . Dgamma[2] . m . Maple's implementation of Matrix is such that - e.g. - evalb(Matrix(2,[a,b]) = Matrix(2,[a,b])) returns false. The composition of three of these objects put an internal subroutine in trouble. Anyway, it is fixed now. Thanks for posting the problem.

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

Christian, Kambiz1,

A quick look suggests to me that this problem has no exact and closed form solution where the arbitrary constants and functions could be tweaked to satisfy something like  T(0, y, t) = 1, T(1, y, t) = 1, T(x, 0, t) = 1, T(x, 1, t) = 1, T(x, y, 0) = 0.Those conditions are too many and too specific. My impression is that, if a solution to this problem exists, it is an infinite series solution.


One approach is then, always, to investigate first that part of the problem, the set of initial conditions. This is your problem after converting it to an exact form (no floating point numbers here):

pde := convert(diff(T(x, y, t), x, x)+diff(T(x, y, t), y, y) = (diff(T(x, y, t), t))/(0.44e-2), rational, exact)
conds := T(0, y, t) = 1, T(1, y, t) = 1, T(x, 0, t) = 1, T(x, 1, t) = 1, T(x, y, 0) = 0

diff(diff(T(x, y, t), x), x)+diff(diff(T(x, y, t), y), y) = (2272727273/10000000)*(diff(T(x, y, t), t))

 

T(0, y, t) = 1, T(1, y, t) = 1, T(x, 0, t) = 1, T(x, 1, t) = 1, T(x, y, 0) = 0

(1)

The question is: what would be the most general form of initial conditions compatible with the given pde (i.e.: such that the problem is consistent, has a solution regardless of whether we can compute it).

 

The command for that is DEtools[initialdata]

IC := DEtools[initialdata]([isolate(pde, diff(T(x, y, t), x, x))])

table( [( Finite ) = [], ( Infinite ) = [T(x[0], y, t) = _F1(y, t), (D[1](T))(x[0], y, t) = _F2(y, t)] ] )

(2)

So, for this problem to be consistent, the initial conditions depend at most on two arbitrary functions and are of the following form and the problem has an exact solution:

conds_at_x0 := IC[Infinite]

[T(x[0], y, t) = _F1(y, t), (D[1](T))(x[0], y, t) = _F2(y, t)]

(3)

pdsolve([pde, op(conds_at_x0)], T(x, y, t))

T(x, y, t) = _F1(y, t)+Sum((x-x[0])^(2*n)*((proc (U) options operator, arrow; -(diff(diff(U, y), y))+(2272727273/10000000)*(diff(U, t)) end proc)@@n)(_F1(y, t))/factorial(2*n), n = 1 .. infinity)+Sum((x-x[0])^(2*n1+1)*((proc (U) options operator, arrow; -(diff(diff(U, y), y))+(2272727273/10000000)*(diff(U, t)) end proc)@@n1)(_F2(y, t))/factorial(2*n1+1), n1 = 0 .. infinity)

(4)

Voila! The solution got computed using the most abstract method implemented, which requires computing a differential operator, proc (U) options operator, arrow; -(diff(U, y, y))+(2272727273/10000000)*(diff(U, t)) end proc, and then each term of the infinite series is involves its n composition, e.g.:

u := proc (U) options operator, arrow; -(diff(U, y, y))+(2272727273/10000000)*(diff(U, t)) end proc

proc (U) options operator, arrow; -(diff(U, y, y))+(2272727273/10000000)*(diff(U, t)) end proc

(5)

(u@@n)(F(y, t))

(u@@n)(F(y, t))

(6)

eval((u@@n)(F(y, t)), n = 1)

-(diff(diff(F(y, t), y), y))+(2272727273/10000000)*(diff(F(y, t), t))

(7)

eval((u@@n)(F(y, t)), n = 2)

diff(diff(diff(diff(F(y, t), y), y), y), y)-(2272727273/5000000)*(diff(diff(diff(F(y, t), t), y), y))+(5165289257438016529/100000000000000)*(diff(diff(F(y, t), t), t))

(8)

So if you can tweak this general form of the initial conditions (3) to match your problem, then you have a solution here. And then, if an approximate solution (with a finite number of terms) is of use for you, you can simply substitute infinity by a finite number and perform the summation to obtain an exact form.

 

Likewise, since conds_at_x0 := [T(x[0], y, t) = _F1(y, t), (D[1](T))(x[0], y, t) = _F2(y, t)] has only one expansion point x = x__0 (and note that for y and t you have arbitrariness, so you can tweak the problem almost at will regarding them), then this problem is also suitable for a direct computation of a truncated series solution; here is up to order 2 but of course you can compute them to arbitrary order - the advantage with respect to the exact infinite sum (4) being that you do not have the composed differential operator:

pdsolve([pde, op(conds_at_x0)], T(x, y, t), series, order = 2)

T(x, y, t) = _F1(0, 0)+(D[2](_F1))(0, 0)*t+(D[1](_F1))(0, 0)*y+_F2(0, 0)*(x-x[0])+(1/2)*(D[2, 2](_F1))(0, 0)*t^2+(D[1, 2](_F1))(0, 0)*y*t+(D[2](_F2))(0, 0)*(x-x[0])*t+(1/2)*(D[1, 1](_F1))(0, 0)*y^2+(D[1](_F2))(0, 0)*(x-x[0])*y+(1/2)*(-(D[1, 1](_F1))(0, 0)+(2272727273/10000000)*(D[2](_F1))(0, 0))*(x-x[0])^2

(9)

Finally, independent of of all the above, while pdsolve is the main guy to compute exact solutions, the fact is that Maple has more PDE solvers, yet to be integrated properly into pdsolve's (large ...) strategy; the main one is PDEtools:-InvariantSolutions ; here are, then, other exact solutions for pde without initial conditions

PDEtools:-InvariantSolutions(pde, dependency = {t, x, y})

T(x, y, t) = _C1/(t*exp((2272727273/40000000)*(x^2+y^2)/t)), T(x, y, t) = _C1+Ei(1, (2272727273/40000000)*(x^2+y^2)/t)*_C2, T(x, y, t) = (_C1*t+_C2*y)/(t^2*exp((2272727273/40000000)*(x^2+y^2)/t)), T(x, y, t) = (_C1*t+_C2*x)/(t^2*exp((2272727273/40000000)*(x^2+y^2)/t)), T(x, y, t) = (_C2*ln(t/(x^2+y^2)^(1/2))+_C1)/(t*exp((2272727273/40000000)*(x^2+y^2)/t)), T(x, y, t) = (2000*(t/x^2)^(1/2)*exp(-(2272727273/40000000)*x^2/t)*22727272730^(1/2)*_C2+2272727273*erf((1/20000)*22727272730^(1/2)/(t/x^2)^(1/2))*Pi^(1/2)*_C2-2272727273*Pi^(1/2)*_C2+_C1)/((t/x^2)^(1/2)*exp((2272727273/40000000)*y^2/t)), T(x, y, t) = (2000*(t/y^2)^(1/2)*exp(-(2272727273/40000000)*y^2/t)*22727272730^(1/2)*_C2+2272727273*erf((1/20000)*22727272730^(1/2)/(t/y^2)^(1/2))*Pi^(1/2)*_C2-2272727273*Pi^(1/2)*_C2+_C1)/((t/y^2)^(1/2)*exp((2272727273/40000000)*x^2/t))

(10)

Check the help page of this command, it has too many options to mention here, and several of them may take you to more exact solutions - all of them of particular type (these are the so-called Lie symmetry invariant solutions) .


 

Download PDE_solution_for_more_general_ICs.mw

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

Hi

First, note that Physics comes with an implementation of Pauli matrices. The command is called Psigma and has a help page. Then note that these matrices satisfy both commutation and also anticommutation relations, and you didn't enter anticommutations for your J[n]. A third comment: if in your (your equation (3)) you replace J by Psigma, the final result S1 is basically what you expected; , but for a weakness: for Simplify(S) you get S1 := -Psigma[1]^2 - Psigma[3]^2 + 2, and then Simplify(S1) returns 0 as expected. I will fix this weakness (it should have returned 0 directly on Simplify(S)) today or tomorrow and post the fix within the next Physics Updates (works with Maple 2019).

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

In general, for questions about DE manipulation, take a look at ?PDEtools, or ?DEtools. Both have alphabetical lists of commands with a one-line description, including links to the corresponding help pages. If you don't find the answer there, take a look at ?PDEtools,Library.

For this particular question, see ?PDEtools:-dcoeffs. The following does what you want with regards to the functions u and vmap(PDEtools:-dcoeffs, [PDEtools:-dcoeffs((lhs - rhs)((3)), u(x, y, z))], v(x, y, z))Of course you get more than just equations on F, P and Q, so to that you may want to asd select(has,%, {F, P, Q}) .

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

Hi Pooyan,
Opening your worksheet I see that your Vr(Rmax(theta, beta), theta, beta) is actually equal to 

Vr(((4.000000000*cos(theta) + sqrt(17 + 8*cos(2*theta)))*(9 + 20*beta/Pi))/9, theta, beta) = 0

The problem being that the first 'independent' variable in Vr(r, theta,beta), that is 'r', cannot depend on theta and beta, so I imagine there is some problem with your formulation, and that is what the error message is telling you: it is unexpected to find {beta, theta} in the first operand of Vr(..., theta, beta).

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

Hi Rouben

I don't recall seeing this problem before ... It is fixed and the fix available to everybody within the Physics Updates v.446 or higher.

Details: at a very low level, Maple does some automatic simplifications: a + 0 is transformed into a, while a + 0.0 is not (for natural reasons), but then 0.0 * _j is automatically transformed into 0.0. Your input 3.0 _i + 0.0 _j thus arrives at `+` as 3.0 _i + 0.0. Independent of that, Vectors is programmed to work with exact vectors, which doesn't prevent it from working with floats, but then an interpretation for the meaning of 3.0 _i + 0.0 is necessary. I coded now to take that as equal to 3.0 _i, so that in the context of adding vectors, 0.0 is the 0 vector.

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

Hi mwahab

Unfortunately, Mapleprime's mechanism for uploading and displaying the contents of a worksheet is not working, so below, you have only the link to this reply.

In brief, what you have done incorrectly is to assume that splitting the determining system into cases concerning k (so, not pde, but its determining system) would take you to the symmetries of each case of splitting pde into cases with respect to k. Although that may work in some cases, it is not assured that it will work in all cases.

In the worksheet attached I show three things: 1) and 2) are concretely different -however equivalent - approaches to split the symmetries into cases with respect to parameters, then 3) is how to tell when the approach you used will work. 

When the Mapleprimes mechanism to upload and display the contents of a worksheet is restored, I will repost this answer with the contents visible. Your question is indeed interesting for a wider audience, I think.

Download How_to_split_symmetries_into_cases.mw

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

 


 

I can't help you with a grtensor issue but Maple already comes with routines more advanced and easier to use than grtensor; what follows then illustrates how you could do your computation using the official-Maple Physics package; for more info you can take a look to its help page, ?Physics

"f(t):=a(t)^(2);"

proc (t) options operator, arrow, function_assign; a(t)^2 end proc

(1)

with(Physics)

Setup(signature = `+---`, coord = cartesian, metric = dt^2-f(t)*(dx^2+dy^2+dz^2))

[coordinatesystems = {X}, metric = {(1, 1) = 1, (2, 2) = -a(t)^2, (3, 3) = -a(t)^2, (4, 4) = -a(t)^2}, signature = `+ - - -`]

(2)

Check the metric

g_[]

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

(3)

The Ricci tensor is Ricci[mu,nu]. Its definition, one component, all the nonzero, and its matrix form

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]

(4)

With compact display, derivatives are displayed indexed (to have them displayed in non-compact form, input OFF for details see CompactDisplay )

Ricci[1, 1]

-3*(diff(diff(a(t), t), t))/a(t)

(5)

Ricci[nonzero]

Physics:-Ricci[mu, nu] = {(1, 1) = -3*(diff(diff(a(t), t), t))/a(t), (2, 2) = 2*(diff(a(t), t))^2+a(t)*(diff(diff(a(t), t), t)), (3, 3) = 2*(diff(a(t), t))^2+a(t)*(diff(diff(a(t), t), t)), (4, 4) = 2*(diff(a(t), t))^2+a(t)*(diff(diff(a(t), t), t))}

(6)

Ricci[]

Physics:-Ricci[mu, nu] = Matrix(%id = 18446744078320134734)

(7)

"Ricci[~]"

Physics:-Ricci[`~mu`, `~nu`] = Matrix(%id = 18446744078195215886)

(8)

Define a tensor as the product of a scalar and a covariant derivative (following your worksheet)

fyr[mu, nu] = fY(t)*Ricci[mu, nu]

fyr[mu, nu] = fY(t)*Physics:-Ricci[mu, nu]

(9)

Define(fyr[mu, nu] = fY(t)*Physics[Ricci][mu, nu])

`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], fyr[mu, nu], Physics:-g_[mu, nu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(10)

Check it out

fyr[definition]

fyr[mu, nu] = fY(t)*Physics:-Ricci[mu, nu]

(11)

fyr[1, 1]

-3*fY(t)*(diff(diff(a(t), t), t))/a(t)

(12)

fyr[nonzero]

fyr[mu, nu] = {(1, 1) = -3*fY(t)*(diff(diff(a(t), t), t))/a(t), (2, 2) = fY(t)*(2*(diff(a(t), t))^2+a(t)*(diff(diff(a(t), t), t))), (3, 3) = fY(t)*(2*(diff(a(t), t))^2+a(t)*(diff(diff(a(t), t), t))), (4, 4) = fY(t)*(2*(diff(a(t), t))^2+a(t)*(diff(diff(a(t), t), t)))}

(13)

fyr[]

fyr[mu, nu] = Matrix(%id = 18446744078328236022)

(14)

"fyr[~]"

fyr[`~mu`, `~nu`] = Matrix(%id = 18446744078328191686)

(15)

The covariant derivative is D_[mu]

D_[mu](fyr[mu, beta])

Physics:-D_[mu](fyr[beta, `~mu`], [X])

(16)

The components of arbitrary tensorial expressions are visible using TensorArray

TensorArray(Physics[D_][mu](fyr[beta, `~mu`], [X]))

Array(%id = 18446744078349157910)

(17)

"simplify(?)"

Array(%id = 18446744078349166342)

(18)

The double contraction

D_[mu](D_[nu](fyr[mu, nu]))

Physics:-D_[mu](Physics:-D_[nu](fyr[`~mu`, `~nu`]), [X])

(19)

Summing over the repeated indices

SumOverRepeatedIndices(Physics[D_][mu](Physics[D_][nu](fyr[`~mu`, `~nu`]), [X]), simplifier = simplify)

(-3*(diff(diff(diff(diff(a(t), t), t), t), t))*fY(t)*a(t)^2+(-6*a(t)^2*(diff(fY(t), t))-9*a(t)*fY(t)*(diff(a(t), t)))*(diff(diff(diff(a(t), t), t), t))-3*fY(t)*(diff(diff(a(t), t), t))^2*a(t)+(-3*a(t)^2*(diff(diff(fY(t), t), t))-9*a(t)*(diff(a(t), t))*(diff(fY(t), t))+15*fY(t)*(diff(a(t), t))^2)*(diff(diff(a(t), t), t))+6*(diff(a(t), t))^3*(diff(fY(t), t)))/a(t)^3

(20)

Simplifying properties (first set indicates symmetires, second set antisimmetries)

Library:-GetTensorSymmetryProperties(fyr)

{[1, 2]}, {}

(21)

LeviCivita[mu, nu, alpha, beta]*fyr[alpha, beta]

Physics:-LeviCivita[alpha, beta, mu, nu]*fyr[alpha, beta]

(22)

Simplify(Physics[LeviCivita][alpha, beta, mu, nu]*fyr[alpha, beta])

0

(23)

Raising and lowering indices

"g_[~mu,~nu]  fyr[nu,alpha]"

Physics:-g_[`~mu`, `~nu`]*fyr[alpha, nu]

(24)

Simplify(Physics[g_][`~mu`, `~nu`]*fyr[alpha, nu])

fyr[alpha, `~mu`]

(25)

Or in one go using the `.` product operator instead of `*`

"g_[~mu,~nu] .  fyr[nu,alpha]"

fyr[alpha, `~mu`]

(26)

I imagine the above suffices for you to develop your computation; if not, could you please give more details (that do not require having grtensor installed).

 

NULL

NULL

Warning: grOptionMetricPath has not been assigned.

Default metric is now fr.

 

grcalc(R(dn, dn))

`CPU Time ` = 0.16e-1

(27)

NULL

grdef("fyr{c d}:=R{c d}*fY(t) ")

Created definition for fyr(dn,dn)

 

 

grcalc(fyr(dn, dn))

`CPU Time ` = 0.

(28)

 

NULLgrdisplay(fyr(dn, dn))

fyr[a]*``[b] = (array( 1 .. 4, 1 .. 4, [( 2, 1 ) = (0), ( 4, 1 ) = (0), ( 1, 4 ) = (0), ( 1, 2 ) = (0), ( 2, 4 ) = (0), ( 3, 2 ) = (0), ( 3, 4 ) = (0), ( 3, 1 ) = (0), ( 4, 3 ) = (0), ( 3, 3 ) = ((2*(diff(a(t), t))^2+a(t)*(diff(diff(a(t), t), t)))*fY(t)), ( 2, 2 ) = ((2*(diff(a(t), t))^2+a(t)*(diff(diff(a(t), t), t)))*fY(t)), ( 2, 3 ) = (0), ( 4, 2 ) = (0), ( 1, 3 ) = (0), ( 4, 4 ) = ((2*(diff(a(t), t))^2+a(t)*(diff(diff(a(t), t), t)))*fY(t)), ( 1, 1 ) = (-3*(diff(diff(a(t), t), t))*fY(t)/a(t))  ] ))

(29)

grcalc(Box[fyr(dn, dn)])

`CPU Time ` = 0.63e-1

(30)

NULL

NULL

grcalc(fyr(dn, dn, cup, cup))

`CPU Time ` = 0.63e-1

(31)

NULLNULL

Parse could not create a functional expression for:

  "Tensor_(rc,[dn,dn,dn],[c,d,e],0)≔Tensor_(R,[dn,dn,cdn],[c,d,e],0)"
Please recheck the definition.

Error, (in grtensor:-grF_strToDef) parse() failed

 

This is generic definition of the generic VcVd (fY·Rab)

grdef("kt2{a b c d}:=fyr{a b;c d}")

This object is already defined. The new definition has been ignored.

 

grdisplay(kt2(dn, dn, up, up))

Created definition for kt2(dn,dn,up,up)

kt2(dn,dn,up,up) has not been calculated.

 

grcalc(kt2(up, up, dn, dn))

`CPU Time ` = 0.

(32)

now i want to define contracted quantity  using contraction with upper metrics g^{cd}

grdef("kt3:=kt2{abcd}*g{^cd}*g{^ab}")

Indices in name: [[], []]

Indices in definition: [[c], []]

Error, (in grtensor:-grdef) lhs/rhs index conflict.

 

NULL 

grdef("kt4{a c}:=kt2{abcd}*g{^bd}")

Indices in name: [[], [a, c]]

Indices in definition: [[b], [a]]

Error, (in grtensor:-grdef) lhs/rhs index conflict.

 

``


 

Download grtensor_v2_(reviewed).mw

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

It suffices to equate to 0 the functional derivative, which you can compute with Physics:-Fundiff. Now, the Maple input line you show is not correct on several counts; for one you are using [] as () (in Maple, [] are used to enter lists of objects), then e^ is probably meant to be exp(). I'd suggest you to try to write the problem correctly, then give a look at the help page ?Physics,Fundiff.

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

Hi WouterSe
I adjusted the solution you posted to be, in 1D Maple input syntax,

sol_Wouter := `ξr`(r,t) = sin(1/2*Pi*r)*cos(1/2*(kappa/Mu)^(1/2)*Pi*t);

Then tried

pdetest(sol_Wouter, [eqn, ic, bc]);

              [      /1     \      /Pi \         ]
              [0, sin|- Pi r| - sin|---|, 0, 0, 0]
              [      \2     /      \2 r/         ]

So the first ic is not satisfied. Indeed, your first ic reads  `ξr`(r, 0) = sin(Pi/(2*r)), not sin(Pi*r/2). Of course, if you change  sin(Pi/(2*r)) by sin(Pi*r/2) in sol_Wouter, then the PDE (your eqn) does not cancel. 

So, the solution you are expecting is not correct for the problem you posted.

Guessing what could be wrong in your post, from the output by pdetest above, if you change the right-hand side of your first ic from sin(Pi/(2*r)) to sin(Pi*r/2), then the (corrected above) solution you posted cancels all of eqn, ic and bc. But then pdsolve also returns the simpler solution:

# change your ic from sin(Pi/(2*r) to sin(Pi*r/2)

ic := `ξr`(r, 0) = sin(Pi*r/2), D[2](`ξr`)(r, 0) = 0:

pdsolve([eqn, ic, bc], Zeta(r, t));

        `ξr`(r,t) = sin(1/2*Pi*r) * cos(1/2*kappa^(1/2)*Pi*t/Mu^(1/2))

In summary: to get from pdsolve the "simple" solution, you need to correct your input (the first ic). For the input you presented, there is no simple solution that I could see, and the infinite sum solution returned by pdsolve is correct.

 

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

Hi

I don't fully understand the formulation you are trying to do - I mean the problem itself, not how to represent it in Maple.

 

Generally speaking, using Physics is probably simpler for this kind of problem, what I understood from your post is that you are working with cartesian coordinates, in two dimensions, with this metric

with(Physics)

Setup(coordinates = cartesian, dimension = 2, metric = -dx^2*t^2+dt^2)

 

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

 

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

 

_______________________________________________________

 

[coordinatesystems = {X}, dimension = 2, metric = {(1, 1) = -t^2, (2, 2) = 1}]

(1)

Check it out

g_[]

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

(2)

and you want to compute the geodesics for this spacetime, ie the solutions of

Geodesics()

[diff(diff(t(tau), tau), tau) = -t(tau)*(diff(x(tau), tau))^2, diff(diff(x(tau), tau), tau) = -2*(diff(x(tau), tau))*(diff(t(tau), tau))/t(tau)]

(3)

You can call dsolve on these equations, or directly:

Geodesics(output = solutions)

{t(tau) = _C2*tau+_C3, x(tau) = _C1}, {t(tau) = _C4*((_C1*_C2+_C1*tau+4)*(_C1*_C2+_C1*tau-4))^(1/2), x(tau) = -(1/2)*ln(_C1*_C2+_C1*tau+4)+(1/2)*ln(_C1*_C2+_C1*tau-4)+_C3}

(4)

Now, what do you mean by "plotting these geodesics in polar coordinates using the time as the radius" and what do you mean by "[the time] normalized to our current time equal to unity". You say "to plot starting at "x = 0, t = 1. X is the polar angle" but there is no X here, only the cartesian x and the image does not suffice to understand what you actually want. And, you know, to plot any of the solutions (4) you'd need initial conditions (not just a starting point [x__0, t__0]) fixing the arbitrary constants present in these solutions.

 

If you could give a clear-cut description of what you want - mathematically speaking, disregard for now the Maple way of representing the problem - it can most probably be done with ease, departing from the input/output above, either using plots:-polaplot  or Physics:-TransformCoordinates  and the commands of the plots  package.


 

Download polar_plot_of_geodesics.mw

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

Hi

I am curious about your question ... No, you cannot define a 'tensor of rank 0' because that is a scalar, and everything in Maple is already a scalar. I mean: from the computational point of view, what distinguishes a tensor from everything else are two things: 1) its behaviour under a transformation of coordinates (only applicable to tensors of rank > 0, see ?Physics:-TransformCoordinates) and 2) the sum rule for repeated indices (also for tensors of rank > 0). So, in Maple, you do not need to define a 'scalar' (tensor of rank 0) because, generally speaking, everything is a scalar unless you indicate it is a tensor (of rank > 0).

Trying to figure out what you intend to do with Theta = varepsilon[mu, ~mu](X) , maybe you want to have Theta as a computational representation for varepsilon[mu, ~mu](X) ... ? You can achieve that also via alias(Theta = varepsilon[mu, ~mu](X) ) or assigning as in  Theta := varepsilon[mu, ~mu](X). Or if you want to compute with a scalar function, e.g. Theta(X), just use Theta(X). If you don't want to have the display showing the (X) dependency, use Physics:-CompactDisplay, or even alias(Theta = Theta(X)).

Or perhaps I am not understanding what you have in mind - if that is the case could you please give more details to frame your question and understand what you intend to do with Theta(X) that would require a tensor (of rank 0) definition?

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

Hi

Thanks to all for the feedback. The sense of community grows with that. The problem is fixed and the fix is available to everybody within the Physics Updates v.430 and higher.

Comments: this problem fixed is not related to the Physics Updates or any of its versions - for instance, without it, you get the same result posted by nm. The problem was actually a combination of two facts. 

  1. int(Sum(0, j = 1 .. infinity), tau = 0 .. t) returns Sum(0, j = 1 .. infinity)*t instead of 0. Although not incorrect, this is what I'd call an inconvenient weakness.
  2. Unfortunately, that was happening in a place where the code checks for 0, and Sum(0, j = 1 .. infinity)*t is syntactically different from 0, an then the flow moved forward in a case where it should have stopped.

This is now tracked in the internal database as an issue in int, and it is always good to detect these cases, no matter how surprising they could seem.

Investigating this problem also helped to detect another place where the flow could be streamlined, so now pdsolve works faster in some examples. But for this problem, the result we get with v.430 is still NULL, not a solution.

I will give a look at the suggestion by Rouben, he always has good ideas; anyway, this is what I think is a general approach for this kind of problem: change variables avoiding symmetric bc's, as shown in the worksheet below.

NOTE AFTER POSTING: The approach explained in the worksheet below is already implemented within the Physics Updates version 431, so that the solution shown below as (6) is the current output for pdsolve([pde, ic, bc]).

Physics:-Version()

`The "Physics Updates" version in the MapleCloud is 430 and is the same as the version installed in this computer, created 2019, September 23, 13:35 hours, found in the directory /Users/ecterrab/maple/toolbox/2019/Physics Updates/lib/`

(1)

pde := diff(u(x, t), `$`(t, 2)) = 4*(diff(u(x, t), `$`(x, 2))); ic := u(x, 0) = 0, (D[2](u))(x, 0) = sin(x)^2; bc := u(-Pi, t) = 0, u(Pi, t) = 0

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

 

u(x, 0) = 0, (D[2](u))(x, 0) = sin(x)^2

 

u(-Pi, t) = 0, u(Pi, t) = 0

(2)

With the Physics Updates v.430 we still get no solution

pdsolve([pde, ic, bc])

So change variables removing the symmetric bcs; ie proc (x) options operator, arrow; xi-Pi end proc

tr := {t = tau, x = xi-Pi, u(x, t) = upsilon(xi, tau)}

itr := solve(tr, {tau, xi, upsilon(xi, tau)})

{tau = t, xi = x+Pi, upsilon(xi, tau) = u(x, t)}

(3)

PDEtools:-dchange(tr, [pde, ic, bc])

[diff(diff(upsilon(xi, tau), tau), tau) = 4*(diff(diff(upsilon(xi, tau), xi), xi)), upsilon(xi, 0) = 0, (D[2](upsilon))(xi, 0) = sin(xi)^2, upsilon(0, tau) = 0, upsilon(2*Pi, tau) = 0]

(4)

Solve this problem

pdsolve([diff(diff(upsilon(xi, tau), tau), tau) = 4*(diff(diff(upsilon(xi, tau), xi), xi)), upsilon(xi, 0) = 0, (D[2](upsilon))(xi, 0) = sin(xi)^2, upsilon(0, tau) = 0, upsilon(2*Pi, tau) = 0])

upsilon(xi, tau) = (1/315)*(315*(Sum(16*sin((1/2)*n*xi)*sin(n*tau)*((-1)^n-1)/(Pi*n^2*(n^2-16)), n = 5 .. infinity))*Pi+672*sin(tau)*sin((1/2)*xi)+160*sin((3/2)*xi)*sin(3*tau))/Pi

(5)

Change variables back

PDEtools:-dchange(itr, upsilon(xi, tau) = (1/315)*(315*(Sum(16*sin((1/2)*n*xi)*sin(n*tau)*((-1)^n-1)/(Pi*n^2*(n^2-16)), n = 5 .. infinity))*Pi+672*sin(tau)*sin((1/2)*xi)+160*sin((3/2)*xi)*sin(3*tau))/Pi)

u(x, t) = (1/315)*(315*(Sum(16*sin((1/2)*n*(x+Pi))*sin(t*n)*((-1)^n-1)/(Pi*n^2*(n^2-16)), n = 5 .. infinity))*Pi+672*sin(t)*cos((1/2)*x)-160*cos((3/2)*x)*sin(3*t))/Pi

(6)

And that is it. Verify this solution

pdetest(u(x, t) = (1/315)*(315*(Sum(16*sin((1/2)*n*(x+Pi))*sin(t*n)*((-1)^n-1)/(Pi*n^2*(n^2-16)), n = 5 .. infinity))*Pi+672*sin(t)*cos((1/2)*x)-160*cos((3/2)*x)*sin(3*t))/Pi, [pde, ic, bc])

[0, 0, (1/210)*(105*Pi*cos(2*x)-105*Pi+3360*(Sum(sin((1/2)*n*(x+Pi))*((-1)^n-1)/(n^3-16*n), n = 5 .. infinity))-320*cos((3/2)*x)+448*cos((1/2)*x))/Pi, 0, 0]

(7)

So we have one ic that needs a plotting test:

zero := [0, 0, (1/210)*(105*Pi*cos(2*x)-105*Pi+3360*(Sum(sin((1/2)*n*(x+Pi))*((-1)^n-1)/(n^3-16*n), n = 5 .. infinity))-320*cos((3/2)*x)+448*cos((1/2)*x))/Pi, 0, 0][3]

(1/210)*(105*Pi*cos(2*x)-105*Pi+3360*(Sum(sin((1/2)*n*(x+Pi))*((-1)^n-1)/(n^3-16*n), n = 5 .. infinity))-320*cos((3/2)*x)+448*cos((1/2)*x))/Pi

(8)

From the bc, the relevant interval is x = -Pi .. Pi. and indeed in that interval zero is numerically sufficiently close to 0 up to Digits (= 10)

plot(eval(zero, [infinity = 1000, Sum = add]), x = -Pi .. Pi)

 

``


 

Download fixed_in_v.440_and_approach_to_solve_the_problem.mw

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

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