Theoretical and computational results of ODE differ a lot! Why?

Question

Hello. I have a problem in that my theoretical and practical(computational) calculations differ a lot. in 1 order of magnitude actually.

---

Impulse I needed to jump over the saddle point in one case is approx. equals 0.16 in the other approx. 0.06 I do not understand what is the problem.

---

Theoretical part:

---

I have a potential energy that is : $$U(\theta,\phi) = \frac{1}{2}k_2(2M_0)^2cos(\theta)^2 -\frac{1}{4}k_4(2M_0)^4(sin(\theta)^4cos(\phi)^4+sin(\theta)^4sin(\phi)^4+cos(\theta)^4)$$ $k_2, k_4 >0$,$k_2>>k_4$, let $M_0 = 1/2$

Local maximum(Saddle point) is when $(\theta,\phi) = (\frac{\pi}{2},\frac{3\pi}{4})$. $$U = -\frac{1}{4}k_4 \times \frac{1}{2} $$

Local and global minimum is when $(\theta,\phi) = (\pi/2,\pi/2)$, then pot. energy is: $$U = - \frac{1}{4} k_4$$

Logicaly, if I have a particle on the bottom of the potential well and know it`s energy and know the energy of the saddle point I need to go over, the kinetic energy is the difference: $$\Delta U = \frac{1}{8} k_4$$

As my kinetic energy is the following: $$ T = \alpha' 4M_0^2 (\dot\theta^2 + sin^2(\theta)\dot{\phi})^2 $$ Hence: Impulse needed to move in the direction of coordinate $\phi$ is : $$ \dot\phi = \sqrt { \frac{k_4}{8\alpha'} }$$ While I get: $\omega_{af} = \frac{k_4}{2\alpha'}, \gamma_{af} = \frac{\alpha}{2\alpha'} $

$$ \dot\phi = \sqrt { \frac{2\alpha'\omega_{af}}{8\alpha'} }$$

$$ \dot\phi = 0.16$$

---

Computational part

---

$$ [ \ddot{\theta} ]= [ sin(\theta)cos(\theta)\dot{\phi}^2 + \frac{k_2}{2 \alpha'}cos(\theta)sin(\theta)+ \omega_{af}^2(sin(\theta)^3cos(\theta)cos(\phi)^4+sin(\theta)^3cos(\theta)sin(\phi)^4 - cos(\theta)^3sin(\theta)) ] - [ 2 \gamma_{af} \dot{\theta} ];$$

$$[\ddot{\phi}] = - [cot(\theta)\dot{\phi}\dot{\theta}] + [\omega_{af}^2(-sin(\theta)^2cos(\phi)^3sin(\phi) +sin(\theta)^2sin(\phi)^3cos(\phi)) ] - [2\gamma_{af} \dot{\phi} - \frac{1}{2} j ]. $$

here is the matlab code:

function xDot = parallelDegreesW42_mistake_v2(t,x,~,w_af,g_af,gcurr)
k_2 = 125;
k_4 = 12.5;
alpha_stroke = 62.5;
xDot  = [x(2); ... %//= theta'

    sin(x(1)) * cos(x(1)) * ((x(4))^2) + ...
    + ( k_2 * sin( 2 * x(1) ) / ( 4 * alpha_stroke)   ) + ...
    (w_af^2) * (   ( sin(x(1)) )^3  *  cos(x(1)) * ( cos(x(3)) )^4 + ( sin(x(1)) )^3 * cos(x(1)) * ( sin(x(3)) )^4 - ( cos(x(1)) )^3 * sin(x(1))  ) + ...
    - 2 * g_af * x(2);

    x(4); ... %//= phi'

    - x(2) * x(4) * cot(x(1)) + ...
    (w_af^2) * (sin(x(1)))^2 * cos(x(3)) * sin(x(3)) *(-(cos(x(3)))^2 + (sin(x(3)))^2) + ...
    - 2 * g_af * (sin(x(1)))^2 * x(4)  + ...
    + 0.5 * gcurr];

and code for plotting:

function f = callDegrees(v_theta, v_phi,)
    %%
    %
    theta = pi/2;
    phi = (3*pi/4+pi/2)/2;
    w_af = 0.1;
    g_af = 0.01;
    gcurr = 0; % this is j in system of ODES;
    x0 = [theta,v_theta,phi,v_phi];
    tspan = 0:0.01:1000;
    options = odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8 1e-8 1e-8]);
    [t,x] = ode45('parallelDegreesW42_mistake_v2',tspan,x0, options, w_af, g_af, gcurr );
    xx1 = x(:,1);%theta
    xx2 = x(:,3);%phi

    xx3 = x(:,2);%V_theta
    xx4 = x(:,4);%V_phi

    plot(t,xx3,'--','LineWidth' ,1);
    xlabel('t');
    ylabel('phi');

extra remarks: $$w_{af} = 0.1, g_{af}=0.01,j=0$$ FINALY, $$ \dot\phi_{byCOMPUTER} = 0.065$$

Answer 1

There is a mistake in my calculations. I did write that $\omega_{af} = k_4/2\alpha',...$. That is a mistake, I should have written $\omega_{af}^2 = k_4/2\alpha'$. After, I will receive the following: $$\phi_{ini} = \sqrt {\omega_{af}^2 / 4} = \omega_{af}/2 = 0.05$$ Which sounds ok. The code was correct.

Thank you Evgeny for helping to solve my problem. I did actually understood a bit, what have you meant by the phrase:

Your mentioned states were steady states in original equations, without dissipation. When you introduce dissipation they are no longer steady states, but new equilibria points appear.