Agree with the first paragraph, disagree with the second. As noted before, there is a radial component of the acceleration vector, equal to [itex]-r \left(\theta '^2+\sin ^2(\theta ) \phi '^2\right)[/itex] and so there is a radial component to the net force on the particle. Motion on the sphere implies a radial component of the force on the particle.
So what is the answer to the OP? The force on the particle is:[tex]F_r(\theta,\phi)=-mr \left(\theta '^2+\sin ^2(\theta ) \phi '^2\right)[/tex][tex]F_\theta(\theta,\phi)=mr(-\cos(\theta)\sin(\theta)\phi'^2+\theta'')[/tex][tex]F_\phi(\theta,\phi)=mr(2\cos(\theta)\theta'\phi'+\sin(\theta)\phi'')[/tex]This is three equations in three unknowns ([itex]F_r,\theta,\phi[/itex]), so solve the last two for [itex]\theta[/itex] and [itex]\phi[/itex] (given initial conditions) to get the equations of motion, then solve the first for [itex]F_r[/itex], if you feel like it.
That's a big job, but the above is brute force. Maybe there is a more elegant way? I think solving it by the Lagrangian method will be tricky, since [itex]F_\theta[/itex] and [itex]F_\phi[/itex] might not be described by a potential function.
Maybe the bottom line is that it's a much more complicated question than Thornton and Marion thought.