7374 Reputation

19 Badges

14 years, 77 days

MaplePrimes Activity

These are replies submitted by ecterrab

@J F Ogilvie 

As I already said, the Maple Heun functions include numerical evaluation routines for them, actually sophisticated ones. That these routines can be improved, as everything else, does not change that fact. I also know that you are well aware of the existence of those numerical evaluation routines - you actually reported a couple of problems with them. Writing here in "fake-news" style saying things you know that are not true does not help your point. Anyway I am out of this conversation.

@Jayaprakash J 

Open Maple, and give a look at the help page. To open the help page, input ?PDEtools,dchange.


The following is a worksheet with the input lines you posted (and please the next time upload a worksheet, see the green arrow to upload it in the Mapleprimes window where you write your post)



I am quoting Define so that we see what is what you are defining

('Define')(K[mu, nu] = Matrix(4, {(1, 1) = exp(nu), (2, 2) = -exp(lambda), (3, 3) = -r^2, (4, 4) = -r^2*sin(theta)^2}, shape = symmetric))

Physics:-Define(K[mu, nu] = Matrix(%id = 18446744078249969894))


Release the delayed evaluation of Define


`Defined objects with tensor properties`


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


But at this point the metric is Minkowski, not Schwarzschild:


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


Consequently, all the following (what you posted) are actually correct

TensorArray(K[mu, `~nu`])

Matrix(%id = 18446744078388367718)



K[mu, nu] = Matrix(%id = 18446744078388365662)



K[`~mu`, `~nu`] = Matrix(%id = 18446744078388356990)


Note also that you have this alternative syntax available


K[mu, `~nu`] = Matrix(%id = 18446744078388353974)



K[`~mu`, `~nu`] = Matrix(%id = 18446744078388356990)


Now if you want to work with the Schwarzschild metric, first you need to set it, the simplest way to do that is:




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


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


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


`Parameters: `[m]




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


The definition of `K__μ,ν`remains unchanged


K[mu, nu1] = Matrix(%id = 18446744078249957726)




K[mu, nu] = Matrix(%id = 18446744078249970974)


But of course K[`~mu`, `~nu`]is now different


K[`~mu`, `~nu`] = Matrix(%id = 18446744078249956166)


What Maple is calculating is what it is telling:

"K[~mu,~nu] = g_[~mu,~alpha] g_[~nu,~beta] K[alpha,beta]"

K[`~mu`, `~nu`] = Physics:-g_[`~alpha`, `~mu`]*Physics:-g_[`~beta`, `~nu`]*K[alpha, beta]


TensorArray(K[`~mu`, `~nu`] = Physics[g_][`~alpha`, `~mu`]*Physics[g_][`~beta`, `~nu`]*K[alpha, beta])

Matrix(%id = 18446744078261614942)


To me all this looks OK. So is this what you were expecting? Or otherwise what is what you think would be wrong in these results?



Regarding an update for Maple 2018, you already have the last one. The "Physics Updates" is an interesting project, useful for everybody, but it is a significant amount of extra work, we cannot afford that for more than the current version of Maple.


Edgardo S. Cheb-Terrab
Physics, Differential Equations andMathematicalFunctions, Maplesoft


Both dsolve and odetest are large software projects, not reduced to "do one thing", and if this is a generic conversation don't expect all the details. Generally speaking things are as said, and trust me: I wrote both dsolve and odetest. If you want to know about the internals the two routines you may want to trace are `odetest/implicit` (returns FAIL for your example, because you removed _C1) and `odetest/explicit`. Also, in this same example, even if you don't remove _C1 from the input to odetest, `odetest/implicit` fails in finding 0,

So, what if one approach fails? Naturally, try other things. There are exactly 1800 lines of code distributed into 20 subroutines below odetest, implementing different catch-up strategies accumulated along the years, taking into account feedback from Maple users. In some cases the catch-up approaches work, in other cases they don't. In this your new (third example) the catch-up works.

Summaryzing again: if you want to test with odetest, don't change the conventions of this code and let it do its job. If you change the conventions, don't expect it to do its job properly, even when, depending on the example, it may still succeed.

Anyway, I can't continue in this thread, need to move to other pending things.

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


There is nothing like random behavior with software. If the solution is explicit, like the one you show now in your 'reply', then there is no need to isolate the integration constant, and so you can change _C1 by whatever (not depending on x). The behaviour that called your attention (that you cannot change _C1 by C[1] and expect odetest to handle the problem) happens only with implicit solutions, for the reasons explained (in those cases the testing approach starts by isolating the integration constant, etc.).  

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


As mentioned in a previous thread with you, there is a distinction between essential and non-essential singular solutions. dsolve's default is to return only essential singular solutions, i.e. those that  cannot be obtained from the general solution by specializing (possibly piecewise) integration constants.  On the other hand, casesplit's default is to return all singular solutions unless you explicitly indicate  singsol = essential . When you solve(ode, y(x)) you are splitting (kinda manually) the original problem into two problems not having essential singular solutions.

I realize as well that you are not familiar with essential singular solutions or that distinction - in the previous thread I mentioned you could give a look a the help pages ?DifferentialAlgebra,RosenfeldGroebner, or the old ?diffalg where this distinction between essential and non-essential singular solutions - not mentioned in your posts in this thread - is explained and is at the root of the difference in behaviors you are pointing at.

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


... it might well be that at this point I have lost perspective to tell what is standard knowledge and what is not. It's been too many years working on steroids (using computer algebra) on different subjects.

a) I won't write a book (enough with the Mini-Course on Computer Algebra for Physicists or the new Maple 2019 guide for Tensor Computations). I enjoy telling stories, but not long stories, and English is not my native language, slowing things down.

b) By all means, there exist textbooks on this, whether I remember their names or not; but more interesting than 'singular solutions' is differential algebra itself, and singular solutions appear naturally as a subproduct. The textbook I studied was for me clear enough (my background is in theoretical physics, not mathematics): Differential Algebra - Joseph Fels Ritt. The book is rather small and straightforward.

c) For me, textbooks are at the end; first, there is Google followed by my own exploration on the subject. Through Google, the first link is Wikipedia, and the third link I see points to "Singular solutions of differential equations." The presentation there is for me clear and understandable, although it misses that singular solutions may also be particular cases of general solutions. 

d) You may prefer to read about singular solutions directly in your Maple help pages - see "Analysis of singular solutions" in the help page ?DifferentialAlgebra,RosenfeldGrobner (don't let these names intimidate you). I wrote DifferentialAlgebra and based the presentation on singular solutions you see in that help page after eq.(27) on the one of the old ?diffalg,deprecated (page written by Evelyne Hubert). There you see a discussion of singular solutions that are not particular cases of the general solution.

e) About reading Maple code that is 'way too advanced', I can't help, and in fact I do not recommend that in this case. PDEtools:-casesplit is growing since 1998, and by now does too many things to be understood just from reading it (a ton of functionality - check its help page and that of PDEtools:-dpoyform; its size is 7500+ lines of code, not including the differential elimination engines, e.g. DifferentialThomas itself adds 3900+ lines of code). I frankly do not recommend trying to 'read' something of that size.

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

@nm, @Rouben Rostamian@vv@acer,

Take each of the ODEs discussed in this thread and try

> PDEtools:-casesplit(ODE);

and you will see "the cases," that is: the general and singular cases. As you noticed, the solutions to the singular cases are one differential order less, so their solution does not contain arbitrary constants. The discussion about general versus singular cases is standard in the literature (most textbook on DEs talk about them, Clairaut's equation is one example).

On a separate however related issue, the interesting question is "when are these singular solutions not particular cases of the general solution?" And there is the sub-question: "For ODEs? For PDEs?".

The answer to that question in the case of an ODE was part of the Ph.D. thesis of Evelyn Hubert, on essential singular solutions. An implementation of her algorithm entered the Maple library in 1998, and I coded it in Maple's dsolve to compute only essential singular solutions.

The answer to that same question in the general case, i.e. nonlinear PDE systems, as for instance Einstein's equations for the function components of the spacetime metric in curved spacetimes, got implemented in the Maple library two releases ago within the DifferentialThomas package (for which, unfortunately, I still haven't finished writing all of its help pages).

Summarizing, use PDEtools:-casesplit to see the cases (try its option caseplot to see a graphical representation of all the special cases you are noticing with the ODEs posted here), and remember that singular cases are sometimes particular cases of the general case, but there is now algorithmic ways, implemented in Maple, to compute only the singular cases that are not particular cases of the general case, and dsolve is expected to only return essentially singular cases (modulo bugs, of course).

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

@Carl Love@Rouben Rostamian@Oliveira

This is, again, all about using differential elimination, something implemented in Maple in 1996 (diffalg), then improved with DEtools[rifsimp], DifferentialAlgebra and DifferentialThomas - all differential elimination packages. Exploring the potentiality of those packages there are the higher level functionality commands PDEtools[dpolyform] and PDEtools[casesplit], finally on top of those there is the triad dsolve, pdsolve and PDEtools:-Solve - the later combining into one all of solve, dsolve, pdsolve, fsolve, and also the options series and numeric of these commands.

Surprisingly, it is the year 2019, and still today Maple is the only Computer Algebra System (CAS) that has differential elimination implemented. The piramid of funcionality just mentioned is at the root of Maple's top performance with differential equations when compared with any other CAS.

The two paragraphs above turn on the light on what is going on, that experienced mathematicians and computer algebra users get surprised with this reply by Rouben today. The explanation may look abstract but tell the story; the corresponding help pages the details.

The other two things I wanted to mention: your comment, Rouben, in a previous Mapleprimes reply, about rewriting the help page for PDEtools:-casesplit (regarding usage, this command is the most important one of all those mentioned) is for me, not a clear guideline, I suppose because I wrote these command (but for the original diffalg and rifsimp). It would help if any of you would like to rewrite part or all of that help page according to what you think would be best and post it so that I could see these commands "through your eyes," then see how to integrate your suggestions in these pages.

Finally, Carl, the idea is a simple one: what you call "a scalar" is a function with zero derivative; i.e. replace omega by omega(x) and add the equation diff(omega(x), x) = 0 to the system being solved. All these commands mentioned also handle inequations, not just equations, and accept all the optional arguments you can pass to PDEtools:-casesplit; take a look at its help page you will realize you can compute very tricky things with simple instructions - one of my favorite ones: the option 'independentof' of the PDEtools:-Solve command.

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


... That when assummingusesAssume = true dsolve returns only one branch is a bug. I will give a look at this today.

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

@Carl Love @nm

Another way of rewriting a differential equation in differential polynomial form, that is polynomial in the unknown and its derivatives and with rational coefficients is using PDEtools:-dpolyform with the option no_Fn.

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

Could you please (always) attach the worksheet with the contents you show? For that purpose you can use the Green arrow. Besides that, to give values to the _Fn see ?eval, and if you open the help page for Infinitesimals (always check the help pages before using commands) right at the top you see the option specialize_Fn. You may also want to take a look at the help page of CharacteristicQInvariants.

I can't help you with that. Anyway, for the example you indcate, in my opinion the degree is 1, not 2, since it is linear in the highest derivative, the rest of the terms are not relevant regarding the degree (see the examples in the literature referenced in the wiki page you posted).

Provided the DE you pass is of integer degree in the highest order derivatives (could be an ODE or a PDE), the following procedure gives you the degree:

de_degree := proc(de, f::expects(unknown) := NULL)
local de_in_diff_notation, derivatives, diff_ord;

de_in_diff_notation := convert(`if`(de::`=`, (lhs-rhs)(de), de), diff);
derivatives := indets(de_in_diff_notation, specfunc(diff));
if f <> NULL then
    derivatives := select(has, derivatives, f);
if derivatives = {} then return(0) fi;
diff_ord := PDEtools:-difforder(derivatives);
derivatives := select(u -> PDEtools:-difforder(u) = diff_ord, derivatives);
max(map2(degree, de_in_diff_notation, derivatives));


> DE := (1 + diff(y(x, t), x)^2)^(3/2) = diff(y(x, t), x $ 2)*diff(y(x, t), t, t)^3*diff(f(x, t), x $ 3, t $ 3)^2;
> de_degree(DE);


> de_degree(DE, y(x,t));


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


In the previous reply by Rouben he gave a beautiful/skilled presentation showing that: given the two PDEs for Q(x, z, t) and W(x, z, t) and the relationships for the auxiliary functions phi(x, z, t) and psi(x, z, t) as the laplacians of Q and W:

"phi(x,z,t)=((&PartialD;)^(2 )W)/((&PartialD;)^( )x^2)+((&PartialD;)^(2 )W)/((&PartialD;)^( )z^2)",
"psi(x,z,t)=((&PartialD;)^(2 )Q)/((&PartialD;)^( )x^2)+((&PartialD;)^(2 )Q)/((&PartialD;)^( )z^2)".

then these two auxiliary functions phi and psi satisfy wave equations with different velocities

"((&PartialD;)^2phi)/((&PartialD;)^( )t^2)=mu/(rho)(((&PartialD;)^2phi)/((&PartialD;)^( )x^2)+((&PartialD;)^2phi)/((&PartialD;)^( )z^2)),"
"((&PartialD;)^2psi)/((&PartialD;)^( )t^2)=(lambda+2 mu)/(rho)(((&PartialD;)^2psi)/((&PartialD;)^( )x^2)+((&PartialD;)^2psi)/((&PartialD;)^( )z^2)),"

Rouben's presentation used two commands, not very well known by most people, that are really powerful - gems of the Maple library: PDEtools:-casesplit  and simplify/siderels  (with respect to side relations)


From those two commands, the least known one is PDEtools:-casesplit. In this presentation I want to show that this command is so powerful and versatile that once you mastered it you can do surprising things. For instance, achieve Rouben's result with a single call to PDEtools:-casesplit, not using any particularly obscure syntax or knowledge, and at the same time decoupling of the systems of equations for Q(x, z, t) and W(x, z, t) requested by madany when he started this thread.

For readability use a compact display for all these functions

PDEtools:-declare((Q, W, phi, psi)(x, z, t))

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


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


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


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


Start entering the equations pde1 and pde2 posted by madany and the relationship, all written on a Maple worksheet by Rouben

pde1 := rho*(diff(Q(x, z, t), t, t, x))+rho*(diff(W(x, z, t), t, t, z))-(lambda+2*mu)*(diff(Q(x, z, t), x, x, x)+diff(Q(x, z, t), x, z, z))-mu*(diff(W(x, z, t), x, x, z)+diff(W(x, z, t), z, z, z)) = 0

pde2 := rho*(diff(Q(x, z, t), t, t, z))-rho*(diff(W(x, z, t), t, t, x))-(lambda+2*mu)*(diff(Q(x, z, t), x, x, z)+diff(Q(x, z, t), z, z, z))+mu*(diff(W(x, z, t), x, x, x)+diff(W(x, z, t), x, z, z)) = 0

relationship := phi(x, z, t) = diff(W(x, z, t), x, x)+diff(W(x, z, t), z, z), psi(x, z, t) = diff(Q(x, z, t), x, x)+diff(Q(x, z, t), z, z)

Construct now a system with the four equations, not just pde1 and pde2

sys := [pde1, pde2, relationship]

Now think: we want Q and W expressed as functions of psi and phi, plus equations for them not involving Q or W. According to the help page of casesplit, that can be achieved with a ranking of the form [{Q, W}, {phi, psi}]. In addition, we want Q and W decoupled, so according to the help page, for that purpose, your ranking should be of elimination type with respect to Q and W, explained in the page it should then be written as "[Q, W, {psi, phi}]."

That is all.

PDEtools:-casesplit(sys, [Q, W, {phi, psi}])

`casesplit/ans`([diff(diff(diff(Q(x, z, t), t), t), x) = (mu*(diff(phi(x, z, t), z))+lambda*(diff(psi(x, z, t), x))+2*mu*(diff(psi(x, z, t), x))-rho*(diff(diff(diff(W(x, z, t), t), t), z)))/rho, diff(diff(diff(Q(x, z, t), t), t), z) = (lambda*(diff(psi(x, z, t), z))+2*mu*(diff(psi(x, z, t), z))-mu*(diff(phi(x, z, t), x))+rho*(diff(diff(diff(W(x, z, t), t), t), x)))/rho, diff(diff(Q(x, z, t), x), x) = psi(x, z, t)-(diff(diff(Q(x, z, t), z), z)), diff(diff(W(x, z, t), x), x) = phi(x, z, t)-(diff(diff(W(x, z, t), z), z)), diff(diff(psi(x, z, t), x), x) = (-lambda*(diff(diff(psi(x, z, t), z), z))-2*mu*(diff(diff(psi(x, z, t), z), z))+rho*(diff(diff(psi(x, z, t), t), t)))/(lambda+2*mu), diff(diff(phi(x, z, t), x), x) = (-mu*(diff(diff(phi(x, z, t), z), z))+rho*(diff(diff(phi(x, z, t), t), t)))/mu], [])


Note the derivatives displayed compactly, indexed, due to using PDEtools:-declare in (1). The last two equations of (2) above are, precisely, the wave equations for phi and psi, shown by Rouben, and the first four equations are the decoupled version of the system presented by madany, where W satisfies an equation not involving Q, that is, W[x, x] = phi-W[z, z], and the other equations express Q as a function of the rest, all of which have been decoupled (that is the meaning of sequential decoupling, at the root of how dsolve and pdsolve work).



Download mw_(reviewed).mw

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

Before using any command, always give a look to the help page. The second argument you are passing is not according to the help page. By the way this example doesn't need second argument, and the decoupling is sequential, as explained in the help page; ie you will get a subsystem that only involves Q, and assuming you can solve it, there is another subsystem involving W and Q: put there the value you obtain for Q by solving the first subsystem. In this example, the subsystem for Q consists of a single 4th order PDE.

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

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