@loicva: Ah, I see. Then I recommend that you use the second approach of mine, along the following lines:

with(PDETools):
L := ToJet(L,theta2(t));
expr1 := FromJet(diff(L,theta2[t]),theta2(t));   # wrt. angular velocity
expr2 := FromJet(diff(L,theta2   ),theta2(t));   # wrt. an angle

In the Jet calculus, diff(theta2(t),t) becomes theta2[t], and theta2(t) becomes theta2. A mixed derivative diff(theta2(t,x),t,x), say, would become theta2[t,x]; and so on. In general, the coordinates with respect to which theta2(t,...) are differentiated become indices of theta2.

@loicva: You are welcome.

Now I am confused. You write that you take the derivative with respect to theta2(t). First, how are you able to do that?, in view of your own original problem, as posted above. Second, why take the derivative with respect to theta2(t) after the substitution by u has been made?

But perhaps I misunderstand you. Could you provide some code?

@Kitonum: I am not sure that your solution provides an answer to the question posed. As I understand the latter, a general transformation from spherical to cartesian coordinates is desired, such a transformation requirering, of course, inverse trigonometric functions.

Please provide an example of your equations, perhaps an upload (using the green arrow in the editor) of your worksheet, or an appropiate part of it.

@mskalsi: I think all that is needed is changing

ToJet(I,DepVars)

into

I := ToJet(I,DepVars)

PS: Previously I had clearly not understood the purpose or notation of Jet calculus (with which I have had no previous experience). But now I think I understand it. Loosely speaking, it makes functional differentiation possible by 'converting' derivatives of functions to non-functions. Right?

@lemelinm: You are quite welcome.

I agree. It seems to be far too difficult to perform a task that one would somehow think was a relatively simple one. Of course, it might just be down to ignorance on my behalf concerning the full tensor capabilities of the Physics package. Perhaps Edgardo S. Cheb-Terrab, the lead developer on the Physics package, would be able to perform some magic!?

@lemelinm: From the expression you give, I guess you are working in a local inertial system (in which the Christoffel symbols vanish at some specific spacetime point), right? Note, by the way, that you are missing a derivative. The expression may be defined as follows:

R[mu,nu,alpha,beta] := g_[mu,gamma]*(
   +d_[alpha](Christoffel[~gamma,nu,beta ])
   -d_[beta ](Christoffel[~gamma,nu,alpha])
);

But this R[mu,nu,alpha,beta] is not yet a tensor internally to Maple. To define it as such, I think the correct way is as follows (see the help page on Define):

Define(R[mu,nu,alpha,beta] =
   Array((1..4)$4,(mu,nu,alpha,beta) -> g_[mu,gamma]*(
      +d_[alpha](Christoffel[~gamma,nu,beta ])
      -d_[beta ](Christoffel[~gamma,nu,alpha])
   ))
):

The problem, though, is that it takes absurdly long time; in fact, Maple did never finish before I interrupted it, and that even if the above definition was augmented by the appropiate (anti)symmetries of the Riemann tensor, as available through Define(Riemann,query). So I am stuck for the moment.

PS: I think you have misunderstood me concerning the function linearize. I expect it to work as well for second order derivatives. Maple gets bogged down by the calculation of the Riemannian components themselves, prior to any linearization.

@Preben Alsholm: Sorry for the late reply. I first see your comment now. I have noted the update myself. Although not precisely as before, the update is acceptable.

@lemelinm: You are quite welcome. I am learning new stuff too.

The calculation of the components of the Riemann tensor will require the above linearization function to include as well second order derivatives. The following version, which includes some minor improvements as well, does that:

linearize := proc(expr::algebraic)
   local inds,sub1,sub2;
   # The components of h, and their first two derivatives, i.e.,
   # the indeterminates with respect to which to linearize
   inds := indets(rhs(g_[])) minus {X};
   inds := inds union `union`(
      seq(map(diff,inds,X[rho]
   ),rho = 1..4)) union `union`(
      seq(seq(map(diff,inds,X[rho],X[sigma])
   ,rho = 1..4),sigma = 1..4));
   # Substitution sets for going back and forth
   # between functions and non-functions (in order
   # to be able to apply mtaylor)
   sub1 := {seq(inds[i] = F||i,i = 1..nops(inds))};
   sub2 := {seq(F||i = inds[i],i = 1..nops(inds))};
   # The linearization itself
   eval(mtaylor(eval(expand(expr),sub1),map(rhs,sub1),2),sub2);
end proc:

A word of warning: As before, the function substitutes the components of h, and their derivatives, with the non-functions F1, F2, F3, etc. Thus global variables of yours should preferably not be among these; otherwise, change the name F in the function.

I have just tried to calculate Riemann[~1,2,3,4]. That seems to be a tough job for Maple; in fact, the calculation did not terminate in any 'reasonable' amount of time. So perhaps a brute force approach to the calculation of the components of the Riemann tensor is not really the way to go. But what then?, I ask myself.

@student_md: As vv correctly objects, you will have to give explicit expressions for the functions, at the very least.

@student_md: Sorry to say, but I am still confused. I guess F1,F2,F3,F4 are supposed to be functions of the Xn, n = 1..N, for otherwise the derivatives will, of course, all vanish. But they seem to be given as functions of t only. Also, the derivatives you consider seem to be socalled functional derivatives? Is that intended? Finally, for the infinite sum, you are considering, not to terminate for some finite value N, these functions need, as far as I can see, to depend on infinitely many Xn. Is that really what you want, or have in mind?

@vv: Quite right. I have never used JordanForm, though.

@student_md: I am still confused at what you want to achieve. The variable n must be assigned a specific value first. You could do something like

n := 3:
X := [seq(x||i     ,i = 1..n)]:
Y := [seq(y||i(X[]),i = 1..n)]:
VectorCalculus:-Jacobian(Y,X);

defining X and Y as lists, say. Or, if you do not want the explicit arguments (x1,x2,...,xn) for the functions y1,y2,...,yn, then you could use the inert form, Diff, of diff:

n := 3:
X := [seq(x||i,i = 1..n)]:
Y := [seq(y||i,i = 1..n)]:
Matrix(nops(Y),nops(X),(i,j) -> Diff(Y[i],X[j]));

I am not sure, I understand what you mean. If you know how to define a 6 x 1 matrix, then you also know how to define an n x 1 matrix for any given n. You can just write

n := 10:
Y := Matrix(n,1,symbol = y);

say, and then change n accordingly. An n x 1 matrix cannot be defined without n having first been assigned some (positive integer) value.

@student_md: I have been away from my computer. I guess Kitonum has satisfactorily answered your question in the meantime.