I'm stuck in a problem where I must simulate the trajectory of a spherical pendulum given the initial conditions. I will start with the things i have already done, then comes the question.
To formalize the problem we have $x(t) = \left [x_1(t), x_2(t), x_3(t)\right ]^T$ that describe the movement of a pendulum of mass $m = 1$ kg in a spherical surface $\phi(x) = x{1}^{2} + x{2}^{2} + x{3}^{2} -1 = 0$. Take $g = 9.8 m/s^2$ and $F = \left [ 0, 0, -g*m\right ]^T$, now we have the system of ODEs:
$x'' = \frac{1}{m}*(F - \frac{mx'^THx'+\bigtriangledown\phi^TF}{\left| \bigtriangledown\phi\right|^2})\bigtriangledown\phi$
Where $H$ is the Hessian matrix of $\phi$ and $\bigtriangledown\phi$ is the gradient of $\phi$.
This part is OK to me. the problem asks movement in terms of it's aceleration, in this case we must do a change of variables to turn this second-order problem into a first-order problem.
$y_1= {x_1}$
$y_2= {x_2}$
$y_3= {x_3}$
$y_4= {x_1}'$
$y_5= {x_2}'$
$y_6= {x_3}'$
Take now $\lambda = \frac{mx'^THx'+\bigtriangledown\phi^TF}{\left| \bigtriangledown\phi\right|^2}$ then we have the 6 variable system:
${y_1}'= {x_1}' = y_4$
${y_2}'= {x_2}' = y_5$
${y_3}'= {x_3}' = y_6$
${y_4}' = {x_1}'' = \frac{1}{m}(F_1 - \lambda \frac{\partial \phi}{\partial y_1}(t))$
${y_5}'= {x_2}'' = \frac{1}{m}(F_2 - \lambda \frac{\partial \phi}{\partial y_2}(t))$
${y_6}'= {x_3}'' = \frac{1}{m}(F_3 - \lambda \frac{\partial \phi}{\partial y_3}(t))$
Now is the part that I define the function on Octave and use it for the methods: Explicit Euler Method, RK2, RK3 and RK4.
All methods worked as expected, but I cannot find a way to calculate the maximum step-size for each method because the problem is not linear since the value of $\lambda$ depends on all variables. I checked in my handouts about stability but found nothing that can help me. Can someone aid me in finding the solution?
You could, if everything else fails, implement step-size control via step doubling, and take a relatively lax tolerance value, this will naturally give step sizes at the maximum for the stability of the method.
Note that $|\nabla\phi|^2=4$, so that the right side is quadratic in the state. This should help to find a Lipschitz constant for the right side.
This then gives--slightly unreliable--estimates $Lh\le1$ for the Euler method, $Lh\le2$ for second-order methods and $Lh\le 2.5$ for RK4 for maximal step sizes that still give somewhat useful results.
With a norm that balances the norms for position and velocity one usually finds smaller Lipschitz constants. Cf. the harmonic oscillator, $\ddot x+w^2x=0$, where a simple norm gives a Lipschitz constant $L=\max(1,w^2)$, vs. a balanced norm where $L=|w|$.