Superficie cilíndrica

Supongamos que una partícula parte en reposo de un punto A (R, φ₀, z₀) situado en la superficie cilíndrica a una altura z₀, formando un ángulo φ₀ con el eje X y llega al punto B (R, 0, 0) situado en el eje X, tal como se muestra en la figura.

Aplicamos el principio de conservación de la energía, para calcular la velocidad de la partícula cuando está a una altura z

${\displaystyle \begin{array}{l}mg{z}_0=\frac{1}{2}m{v}^2+mgz\\ v=\sqrt{2g\left(z_0-z\right)}\end{array}}$

Las coordenadas de un punto sobre la superficie cilíndrica son,
x=Rcosφ, y=Rsinφ, z

El desplazamiento infinitesimal ds a lo largo de una línea contenida en dicha superficie

${\displaystyle d{s}^2=d{x}^2+d{y}^2+d{z}^2=R^{2}d{\phi }^2+d{z}^2}$

El tiempo que emplea la partícula en desplazarse desde A a B es

${\displaystyle T={\displaystyle \underset{A}{\overset{B}{\int }}dt}={\displaystyle \underset{A}{\overset{B}{\int }}\frac{ds}{v}}={\displaystyle \underset{A}{\overset{B}{\int }}\frac{\sqrt{R^2+{\left(\frac{dz}{d\phi }\right)}^2}}{\sqrt{2g\left(z_0-z\right)}}}d\phi }$

Aplicamos la ecuación de Euler-Lagrange a la función f que no depende de φ

${\displaystyle f(\phi ,z,\dot{z})=\frac{\sqrt{R^2+{\dot{z}}^2}}{\sqrt{2g\left(z_0-z\right)}}}$

por lo que utilizamos la segunda versión de dicha ecuación

${\displaystyle \begin{array}{l}\frac{\partial f}{\partial \dot{z}}=\frac{\dot{z}}{\sqrt{R^2+{\dot{z}}^2}}\\ \frac{d}{d\phi }\left(\dot{z}\frac{\partial f}{\partial \dot{z}}-f\right)=0\\ \frac{R^2}{\sqrt{2g\left(z_0-z\right)}\sqrt{R^2+{\dot{z}}^2}}=C\end{array}}$

Donde C es una constante. Elevamos al cuadrado, despejamos $\dot{z}=\frac{dz}{d\phi }$, e integramos respecto de z

${\displaystyle \begin{array}{l}\frac{dz}{d\phi }=\sqrt{\frac{c+R^{2}z}{z_0-z}}\\ \phi -{\phi }_0={\displaystyle \underset{z_0}{\overset{z}{\int }}\sqrt{\frac{z_0-z}{c+R^{2}z}}}dz\end{array}}$

Nota: Las raíces de la ecuación x²=k², son x=±k. Habitualmente, tomamos el signo positivo, pero hay ocasiones en que el resultado es un número negativo

Donde c es una constante. Haciendo el cambio de variable, y=z₀-z, dy=-dz

${\displaystyle {\displaystyle \int \sqrt{\frac{z_0-z}{c+R^{2}z}}dz}=-\frac{1}{R}{\displaystyle \int \sqrt{\frac{y}{k-y}}dy}}$

Donde k es una constante a determinar. Hacemos el cambio de variable

${\displaystyle \begin{array}{l}y=k{\sin }^2\theta \\ dy=2k\sin \theta \cos \theta ·d\theta \\ {\displaystyle \int \sqrt{\frac{y}{k-y}}dy=2k{\displaystyle \int \tan \theta ·\sin \theta ·\cos \theta ·d\theta }=2k{\displaystyle \int {\sin }^2\theta ·d\theta =k{\displaystyle \int \left(1-\cos (2\theta )\right)d\theta =}}}\\ k\left(\theta -\frac{1}{2}\sin (2\theta )\right)\end{array}}$

El resultado es

${\displaystyle R\phi -R{\phi }_0=-k\left(\theta -\frac{1}{2}\sin (2\theta )\right)}$

Cuando el parámetro θ=0, φ=φ₀. La otra ecuación viene del cambio de variable

${\displaystyle \begin{array}{l}{z}_0-z=k{\sin }^2\theta \\ z=z_0-k{\sin }^2\theta \end{array}}$

Cuando el parámetro θ=0, z=z₀

Un resultado similar al obtenido para la braquistócrona en el plano vertical

Las ecuaciones paramétricas de la trayectoria son

${\displaystyle \{\begin{array}{l}\phi -{\phi }_0=-\frac{k}{R}\left(\theta -\frac{1}{2}\sin (2\theta )\right)\\ z=z_0-k{\sin }^2\theta \end{array}}$

Vamos a determinar el valor de la constante k y del parámetro θ cuando la partícula llega a la posición final B (R,0,0), resolviendo el sistema de dos ecuaciones

${\displaystyle \{\begin{array}{l}{\phi }_0=\frac{k}{R}\left(\theta -\frac{1}{2}\sin (2\theta )\right)\\ z_0=k{\sin }^2\theta \end{array}}$

Dividimos la segunda entre la primera

${\displaystyle \frac{z_0}{R{\phi }_0}=\frac{2{\sin }^2\theta }{2\theta -\sin (2\theta )}}$

Calculamos la raíz θ₁ de la ecuación transcendente. De la segunda ecuación, despejamos k.

Tiempo mínimo

Calculamos el tiempo mínimo que emplea la partícula en desplazarse desde A a B a lo largo de dicha trayectoria

${\displaystyle T_m={\displaystyle \underset{A}{\overset{B}{\int }}\frac{\sqrt{R^2+{\left(\frac{dz}{d\phi }\right)}^2}}{\sqrt{2g\left(z_0-z\right)}}}d\phi }$

De las expresiones de z y φ en términos del parámetro θ

${\displaystyle \begin{array}{l}\frac{dz}{d\phi }=\frac{\frac{dz}{d\theta }}{\frac{d\phi }{d\theta }}=\frac{-2k\sin \theta \cos \theta }{-\frac{k}{R}\left(1-\cos (2\theta )\right)}=R\frac{\cos \theta }{\sin \theta }\\ d\phi =-2\frac{k}{R}{\sin }^2\theta ·d\theta \\ T_m=-\sqrt{\frac{2k}{g}}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}d\theta }=-\frac{{\theta }_1}{\sin {\theta }_1}\sqrt{\frac{2z_0}{g}}\end{array}}$

El resultado debería dar un número positivo

Tiempo geodésica

Calculamos el tiempo que emplea la partícula en desplazarse de A a B a lo largo de la geodésica que une ambos puntos.

En la página titulada Líneas geodésicas, demostramos que para una superficie cilíndrica son líneas helicoidales de paso constante, de ecuación φ=c₃z+c₂. Donde las constantes c₂ y c₃ se determinan a partir de las coordenadas de los puntos A y B

φ₀=c₃z₀+c₂
0=c₃·0+c₂

La ecuación de la línea helicoidal sobre la superficie cilíndrica que pasa por A y B es

${\displaystyle \phi ={\phi }_0\frac{z}{z_0},z=z_0\frac{\phi }{{\phi }_0}}$

Expresamos z y su derivada en términos del ángulo φ

${\displaystyle T_g={\displaystyle \underset{A}{\overset{B}{\int }}\frac{\sqrt{R^2+{\left(\frac{dz}{d\phi }\right)}^2}}{\sqrt{2g\left(z_0-z\right)}}}d\phi ={\displaystyle \underset{{\phi }_0}{\overset{0}{\int }}\frac{\sqrt{R^2+{\left(\frac{z_0}{{\phi }_0}\right)}^2}}{\sqrt{2g\left(z_0-\frac{z_0}{{\phi }_0}\phi \right)}}}d\phi =\frac{\sqrt{R^2{\phi }_0^2+z_0^2}}{{\phi }_0\sqrt{2g{z}_0}}{\displaystyle \underset{{\phi }_0}{\overset{0}{\int }}\frac{1}{\sqrt{1-\frac{\phi }{{\phi }_0}}}}d\phi }$

Es una integral inmediata que se resuelve haciendo el cambio de variable u=1-φ/φ₀, du=-dφ/φ₀. El resultado es

${\displaystyle T_g=\sqrt{2}\frac{\sqrt{R^2{\phi }_0^2+z_0^2}}{\sqrt{g{z}_0}}}$

Resultados

Sea una superficie cilíndrica de radio R=1, la posición inicial A es (z₀=12, φ₀=4π). La posición final B es (z=0, φ=0)

Representamos la función transcendente para localizar las raíz θ₁

${\displaystyle f(\theta )=2R{\phi }_0{\sin }^2\theta -z_0\left(2\theta -\sin (2\theta )\right)}$

R=1; %radio
z0=12; %posición inicial A
phi_0=4*pi;

f=@(x) 2*R*phi_0*sin(x).^2-z0*(2*x-sin(2*x));
fplot(f,[0,pi])
xlabel('\theta')
ylabel('f(\theta)')
grid on
title('Ecuación transcendente')

La raíz buscada θ₁ estará comprendida entre θ=1 y 1.5 rad. Utilizamos la función fzero de MATLAB para encontrarla. Después, calculamos la constante k, el tiempo mínimo T_m y el tiempo que emplea la partícula en desplazarse por la geodésica que une los puntos A y B, T_g en unidades del tiempo t de caída libre $z_0=\frac{1}{2}g{t}^2$, para tres valores del ángulo inicial φ₀

• φ₀=4π
• φ₀=6π
• φ₀=8π

R=1; %radio
z0=12; %posición inicial A

for phi_0=[4,6,8]*pi
    f=@(x) 2*R*phi_0*sin(x).^2-z0*(2*x-sin(2*x));
    th_1=fzero(f,[0.2,2]);
    k=z0/sin(th_1)^2;
    Tm=2*th_1/sin(th_1); %en unidades de tiempo de caída libre
    Tg=2*sqrt((R*phi_0/z0)^2+1);
    fprintf('posición, %1.3f, tiempo mínimo, %1.3f, tiempo geodésica, 
%1.3f\n',phi_0, Tm,Tg)
end

posición, 12.566, tiempo mínimo, 2.626, tiempo geodésica, 2.896
posición, 18.850, tiempo mínimo, 3.142, tiempo geodésica, 3.724
posición, 25.133, tiempo mínimo, 3.661, tiempo geodésica, 4.642

Comprobamos que T_g es mayor que T_m

Representamos la trayectoria (color rojo) y la geodésica (color negro) para los tres valores del ángulo inicial φ₀

R=1; %radio
z0=12; %posición inicial A
phi_0=4*pi; %cambiar 6*pi y 8*pi

f=@(x) 2*R*phi_0*sin(x).^2-z0*(2*x-sin(2*x));
th_1=fzero(f,[0.2,2]);
k=z0/sin(th_1)^2;
Tm=2*th_1/sin(th_1);
Tg=2*sqrt((R*phi_0/z0)^2+1);

hold on
r=1*ones(30,1);
phi=linspace(0,2*pi,30);
[r,phi]=meshgrid(r,phi);
x=r.*cos(phi);
y=r.*sin(phi);
z=repmat(linspace(-1,12.5,30),30,1);
surfl(x,y,z);
shading interp
colormap([0.5 0.5 0.5]);

%espiral
th=linspace(0,th_1,400);
phi=-k*(th-sin(2*th)/2)/R+phi_0;
z=z0-k*sin(th).^2;
plot3(cos(phi),sin(phi),z, 'r','lineWidth',1.5) 
phi=linspace(0,phi_0,400);
z=z0*phi/phi_0;
plot3(cos(phi),sin(phi),z, 'k') 

hold off
axis equal
xlabel('x'); ylabel('y'); zlabel('z')
title('Superficie cilíndrica')
view(60,15)

Para ver mejor las trayectorias eliminamos el cilindro de color gris

R=1; %radio
z0=12; %posición inicial A
phi_0=4*pi;

f=@(x) 2*R*phi_0*sin(x).^2-z0*(2*x-sin(2*x));
th_1=fzero(f,[0.2,2]);
k=z0/sin(th_1)^2;
Tm=2*th_1/sin(th_1);
Tg=2*sqrt((R*phi_0/z0)^2+1);

hold on
%espiral
th=linspace(0,th_1,400);
phi=-k*(th-sin(2*th)/2)/R+phi_0;
z=z0-k*sin(th).^2;
plot3(cos(phi),sin(phi),z, 'r','lineWidth',1.5) 
phi=linspace(0,phi_0,400);
z=z0*phi/phi_0;
plot3(cos(phi),sin(phi),z, 'k') 

hold off
axis equal
xlabel('x'); ylabel('y'); zlabel('z')
title('Superficie cilíndrica')
view(60,15)

Observamos que la partícula, en los primeros instantes, cae casi verticalmente, luego, describe la hélice de paso variable. Lo más sorpendente, es que en el tercer caso, φ₀=8π, la partícula se mueve por debajo del origen z<0 y luego, regresa a la posición final z=0. Parece mucho rodeo para un tiempo mínimo. Como podemos comprobar, la posición final B (z=0) se alcanza para el valor del parámetro θ₁=1.7875 y tambien, para π-θ₁=1.3541 (antes)

Estos tres casos se corresponden, con los tres mencionados en el caso del plano vertical

Superficie esférica

Supongamos que una partícula parte en reposo de un punto A (φ₀, θ₀) situado en la superficie esférica de radio R a una altura z₀=Rcosθ₀ y llega al punto B (φ₁, θ₁) tal como se muestra en la figura.

Aplicamos el principio de conservación de la energía, para calcular la velocidad de la partícula cuando está a una altura z

${\displaystyle \begin{array}{l}mg{z}_0=\frac{1}{2}m{v}^2+mgz\\ v=\sqrt{2g\left(z_0-z\right)}\end{array}}$

Las coordenadas de un punto sobre la superficie esférica son
x=Rcosφ·sinθ, y=Rsinφ·sinθ, z=Rcosθ

El desplazamiento infinitesimal ds a lo largo de una línea contenida en dicha superficie

${\displaystyle d{s}^2=d{x}^2+d{y}^2+d{z}^2=R^2{\sin }^2\theta ·d{\phi }^2+R^{2}d{\theta }^2}$

El tiempo que emplea la partícula en desplazarse desde A a B es

${\displaystyle T={\displaystyle \underset{A}{\overset{B}{\int }}dt}={\displaystyle \underset{A}{\overset{B}{\int }}\frac{ds}{v}}=R{\displaystyle \underset{A}{\overset{B}{\int }}\frac{\sqrt{{\sin }^2\theta +{\left(\frac{d\theta }{d\phi }\right)}^2}}{\sqrt{2g\left(z_0-z\right)}}}d\phi =\sqrt{\frac{R}{2g}}{\displaystyle \underset{{\phi }_0}{\overset{{\phi }_1}{\int }}\sqrt{\frac{{\sin }^2\theta +{\left(\frac{d\theta }{d\phi }\right)}^2}{\cos {\theta }_0-\cos \theta }}}d\phi }$

Aplicamos la ecuación de Euler-Lagrange a la función f que no depende de φ

${\displaystyle f(\phi ,\theta ,\dot{\theta })=\frac{\sqrt{{\sin }^2\theta +{\dot{\theta }}^2}}{\sqrt{\cos {\theta }_0-\cos \theta }}}$

por lo que utilizamos la segunda versión de dicha ecuación

${\displaystyle \begin{array}{l}\frac{d}{d\phi }\left(\dot{\theta }\frac{\partial f}{\partial \dot{\theta }}-f\right)=0\\ \frac{{\sin }^2\theta }{\sqrt{{\sin }^2\theta +{\dot{\theta }}^2}\sqrt{\cos {\theta }_0-\cos \theta }}=C\end{array}}$

Donde C es una constante. Elevamos al cuadrado, despejamos $\dot{\theta }=\frac{d\theta }{d\phi }$ e integramos respecto de θ

${\displaystyle \begin{array}{l}\frac{d\phi }{d\theta }=\sqrt{\frac{\cos {\theta }_0-\cos \theta }{c{\sin }^4\theta -{\sin }^2\theta \left(\cos {\theta }_0-\cos \theta \right)}}\\ \phi -{\phi }_0={\displaystyle \underset{{\theta }_0}{\overset{\theta }{\int }}\sqrt{\frac{\cos {\theta }_0-\cos \theta }{c{\sin }^4\theta -{\sin }^2\theta \left(\cos {\theta }_0-\cos \theta \right)}}d\theta }\end{array}}$

Donde c es una constante a determinar de modo que se alcance la posición final B (φ₁, θ₁)

${\displaystyle {\phi }_1-{\phi }_0={\displaystyle \underset{{\theta }_0}{\overset{{\theta }_1}{\int }}\sqrt{\frac{\cos {\theta }_0-\cos \theta }{c{\sin }^4\theta -{\sin }^2\theta \left(\cos {\theta }_0-\cos \theta \right)}}d\theta }}$

Utilizaremos la función fzero de MATLAB para calcular c y el procedimiento integral para calcular la integral entre los límites especificados.

Conocida la constante c, determinamos la trayectoria mediante la ecuación φ-φ₀=f(θ)

Tiempo mínimo

Conocida la constante c calculamos el tiempo mínimo T_m que emplea la partícula en desplazarse a la largo de dicha trayectoria entre A y B

${\displaystyle \begin{array}{l}{T}_m=\sqrt{\frac{R}{2g}}{\displaystyle \underset{{\theta }_0}{\overset{{\theta }_1}{\int }}\sqrt{\frac{{\sin }^2\theta {\left(\frac{d\phi }{d\theta }\right)}^2+1}{\cos {\theta }_0-\cos \theta }}}d\theta =\sqrt{\frac{R}{2g}}{\displaystyle \underset{{\theta }_0}{\overset{{\theta }_1}{\int }}\frac{\sqrt{{\sin }^2\theta \frac{\cos {\theta }_0-\cos \theta }{c{\sin }^4\theta -{\sin }^2\theta \left(\cos {\theta }_0-\cos \theta \right)}+1}}{\sqrt{\cos {\theta }_0-\cos \theta }}}d\theta =\\ \sqrt{c\frac{R}{2g}}{\displaystyle \underset{{\theta }_0}{\overset{{\theta }_1}{\int }}\frac{\sin \theta ·d\theta }{\sqrt{\left(\cos {\theta }_0-\cos \theta \right)\left(c{\sin }^2\theta -\cos {\theta }_0+\cos \theta \right)}}}\end{array}}$

Tiempo geodésica

Calculamos el tiempo que emplea la partícula en desplazarse de A a B a lo largo de la geodésica que une ambos puntos.

En la página titulada Líneas geodésicas, determinamos la ecuación de la línea geodésica sobre una superficie esférica de radio R que pasa por los puntos A y B

${\displaystyle A\{\begin{array}{l}{x}_0=\cos {\phi }_0\sin {\theta }_0\\ y_0=\sin {\phi }_0\sin {\theta }_0\\ z_0=\cos {\theta }_0\end{array}B\{\begin{array}{l}{x}_1=\cos {\phi }_1\sin {\theta }_1\\ y_1=\sin {\phi }_1\sin {\theta }_1\\ z_1=\cos {\theta }_1\end{array}}$

La ecuación es

${\displaystyle \begin{array}{l}a=\arctan \left(\frac{y_0{z}_1-y_1{z}_0}{x_0{z}_1-x_1{z}_0}\right)\\ b=\frac{x_0\sin A-y_0\cos A}{z_0}\\ \tan \theta =\frac{b}{\sin (a-\phi )}\end{array}}$

Damos valores a φ entre 0 y 2π, obtenemos el ángulo θ, dibujamos el punto de coordenadas (φ, θ), unimos los puntos para trazar la circunferencia máxima que pasa por A y B

Despejamos el ángulo φ para calcular la derivada dφ/dθ y obtener el tiempo que emplea la partícula en desplazarse desde A a B a través de la geodésica que pasa por dichos puntos

${\displaystyle \phi =a-\arcsin \left(\frac{b}{\tan \theta }\right),\frac{d\phi }{d\theta }=\frac{b}{{\sin }^2\theta }\frac{1}{\sqrt{1-\frac{b^2}{{\tan }^2\theta }}}}$

El tiempo T_g es

${\displaystyle \begin{array}{l}{T}_g=\sqrt{\frac{R}{2g}}{\displaystyle \underset{{\theta }_0}{\overset{{\theta }_1}{\int }}\sqrt{\frac{{\sin }^2\theta {\left(\frac{d\phi }{d\theta }\right)}^2+1}{\cos {\theta }_0-\cos \theta }}}d\theta =\sqrt{\frac{R}{2g}}{\displaystyle \underset{{\theta }_0}{\overset{{\theta }_1}{\int }}\sqrt{\frac{{\sin }^2\theta \frac{b^2}{{\sin }^4\theta }\frac{1}{1-\frac{b^2}{{\tan }^2\theta }}+1}{\cos {\theta }_0-\cos \theta }}}d\theta =\\ \sqrt{\frac{R}{2g}}\sqrt{b^2+1}{\displaystyle \underset{{\theta }_0}{\overset{{\theta }_1}{\int }}\frac{\sin \theta ·d\theta }{\sqrt{\left({\sin }^2\theta -b^2{\cos }^2\theta \right)\left(\cos {\theta }_0-\cos \theta \right)}}}\end{array}}$

Resultados

Sea una superficie esférica de radio R=1, la posición inicial A es (φ₀=π/3, θ₀=π/6). La posición final B es (φ₁=π/10, θ₀=2π/5)

1. Representamos la función transcendente, φ₁-φ₀±f(c), para localizar la intersección con el eje X (cero(s)). El signo puede ser positivo o negativo

phi_0=pi/3; %posición A
th_0=pi/6;
phi_1=pi/10; %posición B
th_1=2*pi/5;
cc=0:0.01:1.4;
res=zeros(1, length(cc));
i=1;
for c=cc
    f=@(x) sqrt((cos(th_0)-cos(x))./(c*sin(x).^4-sin(x).^2.*
(cos(th_0)-cos(x))));
    res(i)=phi_1-phi_0+integral(f,th_0, th_1);
    i=i+1;
end
plot(cc,res)
grid on
xlabel('c')
ylabel('f(c)')
title('Constante c')



Hay dos posibles valores de c uno próximo a 0.3 y otro próximo a 1. Descartamos el primero ya que obtenemos un valor complejo del tiempo T_m.

2. Utilizamos la función fzero de MATLAB para obtener una mejor aproximación de la constante c

phi_0=pi/3;
th_0=pi/6;
phi_1=pi/10;
th_1=2*pi/5;
f=@(x,c) sqrt((cos(th_0)-cos(x))./(c*sin(x).^4-sin(x).^2.*
(cos(th_0)-cos(x))));
g=@(c) phi_1-phi_0+integral(@(x) f(x,c),th_0, th_1);
c=fzero(g,[0.8,1.2]);
disp(c)

    0.9867

3. Calculamos el tiempo mínimo T_m y el tiempo T_g a lo largo de la geodésica que une los puntos A y B.

R=1; %radio
phi_0=pi/3; %posición A
th_0=pi/6;
phi_1=pi/10; %posición B
th_1=2*pi/5;

%tiempo mínimo
c=0.9867; %constante c
f=@(x) sin(x)./sqrt((cos(th_0)-cos(x)).*(c*sin(x).^2-cos(th_0)+cos(x)));
Tm=sqrt(c*R/(2*9.8))*integral(f,th_0,th_1);

%tiempo geodésica
x0=cos(phi_0)*sin(th_0);
y0=sin(phi_0)*sin(th_0);
z0=cos(th_0);
x1=cos(phi_1)*sin(th_1);
y1=sin(phi_1)*sin(th_1);
z1=cos(th_1);
A=atan2((y0*z1-y1*z0),(x0*z1-x1*z0));
B=(x0*sin(A)-y0*cos(A))/z0;
g=@(x) sin(x)./sqrt((sin(x).^2-B^2*cos(x).^2).*(cos(th_0)-cos(x)));
Tg=sqrt(R/(2*9.8))*sqrt(B^2+1)*integral(g,th_0,th_1);
disp([Tm,Tg])

    0.5997    0.6840

Comprobamos que el tiempo a lo largo de la línea geodésica es mayor

4. Representamos la trayectoria (color rojo) y la geodésica (negro) que pasa por los puntos A y B

%esfera
phi=linspace(0,pi,30);
theta=linspace(0,2*pi,40);
[phi,theta]=meshgrid(phi,theta);
x=sin(phi).*cos(theta);
y=sin(phi).*sin(theta);
z=cos(phi);
hold on
h1=mesh(x,y,z);
set(h1,'EdgeColor',[0.6,0.6,0.6]);

%punto A
phi_0=pi/3;
th_0=pi/6;
x0=cos(phi_0)*sin(th_0);
y0=sin(phi_0)*sin(th_0);
z0=cos(th_0);
plot3(x0,y0,z0,'o','markersize',4,'markeredgecolor','b','markerfacecolor','b')

%punto B
phi_1=pi/10;
th_1=2*pi/5;
x1=cos(phi_1)*sin(th_1);
y1=sin(phi_1)*sin(th_1);
z1=cos(th_1);
plot3(x1,y1,z1,'o','markersize',4,'markeredgecolor','b','markerfacecolor','b')

%trayectoria
c=0.9876; %constante
thh=linspace(th_0,th_1,50);
phi=zeros(1,length(thh));
i=1;
for th=thh
    f=@(x) sqrt((cos(th_0)-cos(x))./(c*sin(x).^4-sin(x).^2.*
(cos(th_0)-cos(x))));
    phi(i)=phi_0-integral(f,th_0,th);
    i=i+1;
end
plot3(cos(phi).*sin(thh),sin(phi).*sin(thh),cos(thh),'r','lineWidth',1.5)

%geodésica
A=atan2((y0*z1-y1*z0),(x0*z1-x1*z0));
B=(x0*sin(A)-y0*cos(A))/z0;
phi=(0:360)*pi/180;
th=atan2(B,sin(A-phi));
plot3(cos(phi).*sin(th),sin(phi).*sin(th),cos(th), 'k','lineWidth',1)

axis equal
axis off
view(110,10)
hold off
xlabel('x'); ylabel('y'); zlabel('z')
title('Superficie esférica')

Conclusiones

El cálculo numérico con MATLAB es muy complejo y produce errores y avisos en muchas ocasiones, se vuelve inestable. Quizás algún lector interesado pueda encontrar la solución a los distintos casos

Referencias

H. A. Yamani, A. A. Mulhem. A cylindrical variation on the brachistochrone problem. Am. J. Phys. 56 (5) May 1988, pp. 467-469

Rafael Xavier Deiga Ferreira, Hector Leny Carrion Salazar. Braquistócrona na esfera. Revista Brasileira de Ensino de Física, vol. 41, n° 4, e20190020 (2019). https://www.scielo.br/j/rbef/i/2019.v41n4/