I'm calculating the geodesics to a parametrized system in **R**^{3}. When trying to solve the geodesic equations for a surface of revolution, I'm getting a strange error. The goal is to write the code for any parametrized surface, not just revolution (hence I didn't use the short cut for surfaces of revolution).

For the code, I find the first fundamental for, then the geodesics, then turn the christoffel symbols into a time-variant array. (Gamma for the position of (u,v) in the uv-plane, and C as the array so I can take derivatives with respect to time.)

The code for creating the Christoffel Symbols and the parametrization of the paraboloid:

restart; with(LinearAlgebra):

V := (u, v) -> <v*cos(u), v*sin(u), v^2+1>;

Christoff := proc (X)

local x1, x2, M, N, i, j, k, s, E, F, G, g, Q, Delta, Prelim, cyclicPrelim;

global Gamma, C; #GAMMA

x1 := (u, v) -> <diff(X[1], u), diff(X[2], u), diff(X[3], u)>;

x2 := (u, v)-> <(diff(X[1], v), diff(X[2], v), diff(X[3], v))>;

E := (u, v) -> DotProduct(x1(u, v), x1(u, v), conjugate = false);

F := (u, v) -> DotProduct(x1(u, v), x2(u, v), conjugate = false);

G := (u, v) -> DotProduct(x2(u, v), x2(u, v), conjugate = false);

simplify([E(u, v), F(u, v), G(u, v)]);

M := (u, v) -> <E(u, v), F(u, v); F(u, v), G(u, v)>;

simplify(M(u,v));

printlevel := 3;

Delta := simplify(Determinant(M(u, v)));

N := (1/Delta)*<G(u, v), -F(u, v); -F(u, v), E(u, v)>;

Q[1] := simplify(map(diff, M(u, v), u));

Q[2] := simplify(map(diff, M(u, v), v));

for i to 2 do for j to 2 do for k to 2 do

Prelim[i, j, k] := Q[k][i, j];

simplify(Prelim[i, j, k] ); end do end do end do;

#(OPTIONAL PRINTOUT) print(Prelim);

for i to 2 do for j to 2 do for k to 2 do

cyclicPrelim[i, j, k] := Prelim[i, j, k]+Prelim[j, k, i]-Prelim[k, i, j] ;

end do end do end do;

#(OPTIONAL PRINTOUT) print(cyclicPrelim);

for i to 2 do for j to 2 do for k to 2 do

Gamma[i, j, k] := simplify((1/2)*add(N[i, s]*cyclicPrelim[j, s, k], s = 1 .. 2));

end do end do end do;

# for k from 1 to 2 do

#print(`GAMMA[i,j,k] =` Matrix([[ `%a` , `%a` ],[ `%a` , `%a` ]])` \n`, Gamma[1,1,k], [Gamma[1,2,k], Gamma[2,1,k], Gamma[2,2,k]);

# end do;

#printf('GAMMA[%a,%a,%a] = %a \n', i,j,k, Gamma[i,j,k]);

print([Gamma[1,1,1], Gamma[1,2,1], Gamma[2,1,1], Gamma[2,2,1], Gamma[1,1,2], Gamma[1,2,2], Gamma[2,1,2], Gamma[2,2,2]]);

for i from 1 to 2 do

for j from 1 to 2 do

for k from 1 to 2 do

C[i,j,k]:= apply(Gamma[i,j,k],t);

end do end do end do;

C:=Array(1..2,1..2,1..2,[ [ [ apply(Gamma[1,1,1],t), apply(Gamma[1,2,1],t) ], [ apply(Gamma[1,1,2],t), apply(Gamma[1,2,2],t) ] ], [ [ apply(Gamma[2,1,1],t), apply(Gamma[2,2,1],t) ], [ apply(Gamma[2,1,2],t), apply(Gamma[2,2,2],t) ] ] ]);

print(C);

end proc;

The differential equation solver:

inits:=[u(0)=1, D(u(t))(0)=1,v(0) = 1, D(v(t))(0)=1];

sys1:= [D[1$2](u(t))+C[1,1,1]*(D(u(t)))^(2)+2*C[1,1,2]*(D(u(t)))*(D(v(t)))+C[1,2,2]*(D(v(t)))^2=0, D[1$2](v(t))+C[2,1,1]*(D(u(t)))^(2)+2*C[2,1,2]*(D(u(t))*D(v(t)))^(2)+C[2,2,2]*(D(v(t)))^2=0];

L:=dsolve({sys1} union {inits});

The error message that comes up is:

"Error, (in unknown) invalid input: op expects 1 or 2 arguments, but received 0"

Any help would be greatly appreciated.