syms x1 x2 x3 x4 g r
syms sys Gamma Gammasecond
sys(1)=x1
sys(2)=x2
sys(3)=x3
sys(4)=x4
%sets g to I
n=3 %dimensions
%metric initialization
for i=1:n
for j=1:n
if(i==j)
g(i,j)=1;
end
end
end
%extra elements for 2 -sphere
%g(1,1)=r^2
%g(2,2)=r^2*(sin(x1))^2
ginv=inv(g)
%Generator of Christoffel symbols of first kind
for p=1:n
for q=1:n
for r=1:n
%Gamma contains Christoffel symbols of the first kind read
%[pq,r]
Gamma(p,q,r)=0.5*(diff(g(q,r),sys(p))+diff(g(r,p),sys(q))-diff(g(p,q),sys(r)));
end
end
end
%Generator of Christoffel symbols of the second kind {i pq}
for m=1:n
for n=1:n
for o=1:n
Gammasecond(m,n,o)=0;
end
end
end
for i=1:n
for p=1:n
for q=1:n
for r=1:n
Gammasecond(i,p,q)=Gammasecond(i,p,q)+ginv(i,r)*Gamma(p,q,r);
end
end
end
end
%generating geodesic equations
%conversion from syms array to string array
for i=1:n
for p=1:n
for q=1:n
Gamma2string(i,p,q)={char(Gammasecond(i,p,q))};
%notice that Gamma2string is a character cell array the{ }
%make it so here on these will be extensively used in this
%program
end
end
end
%creating the equations cell array Eqarr(i) contains the geodesic
%equation for the ith co ordinate
for i=1:n
Eq=strcat('D2x',num2str(i));
for p=1:n
for q=1:n
Eq=strcat(Eq,'+(',Gamma2string(i,p,q),')*Dx',num2str(p),'*Dx',num2str(q))
end
end
Eqarr(i)={Eq}
end
%Extracting the equations from the cells
if n==2
eqn1=cell2mat(Eqarr{1})
eqn2=cell2mat(Eqarr{2})
Paths=dsolve(eqn1,eqn2,'t')
end
if n==3
eqn1=cell2mat(Eqarr{1})
eqn2=cell2mat(Eqarr{2})
eqn3=cell2mat(Eqarr{3})
Paths=dsolve(eqn1,eqn2,eqn3,'t')
end
if n==4
eqn1=cell2mat(Eqarr{1})
eqn2=cell2mat(Eqarr{2})
eqn3=cell2mat(Eqarr{3})
eqn4=cell2mat(Eqarr{4})
Paths=dsolve(eqn1,eqn2,eqn3,eqn4,'t')
end
