Parameter fixing

L = 0.001;
d = 0.001;
l = 0.001;
B = 0.1;
M = 0.001;
k = 1;
R = 0.001;
g = 10;
rx0 = L/2;
rz0 = L/2 + d + l/2;
r0 = Sqrt[rx0^2 + rz0^2];
A0 = -rz0 rx0;
Fz0 = B/r0^5 (rx0 - 2 rz0 + 5 A0/r0^2 rz0);

Spherical coordinates by time dependent angles

rx = L/2 Sin[alpha[t]] Cos[psi[t]];
ry = L/2 Sin[alpha[t]] Sin[psi[t]];
rz = L/2 Cos[alpha[t]] + d + l/2;
r = Sqrt[rx^2 + ry^2 + rz^2];

 A spherical coordinate tensor A of rank 2 and a vector F

In[18]:= A = (rx^2 + ry^2) Cos[alpha[t]] - ry rz Sin[alpha[t]] Sin[psi[t]] - 
   rz rx Sin[alpha[t]] Cos[psi[t]];
Fx = B/r^5 (-4 rx Cos[alpha[t]] + rz Sin[alpha[t]] Cos[psi[t]] + 5 A/r^2 rx);
Fy = B/r^5 (-4 ry Cos[alpha[t]] + rz Sin[alpha[t]] Sin[psi[t]] + 5 A/r^2 ry);
Fz = B/r^5 (ry Sin[alpha[t]] Sin[psi[t]] + rx Sin[alpha[t]] Cos[psi[t]] - 
     2 rz Cos[alpha[t]] + 5 A/r^2 rz);

Three equations  of  of second order for alpha , theta of the form    

D{rx, ry, rz},t,t]== Fx + fs + Nf

In[48]:= D[rx, t, t] /. {_?NumericQ :> 1} // FullSimplify

Out[48]= Sin[alpha[t]] (Cos[psi[t]] + Sin[psi[t]]) Derivative[1][psi][t] + 
 Derivative[1][alpha][
   t] (Cos[psi[t]] (Cos[alpha[t]] + Sin[alpha[t]]) + 
    Cos[alpha[t]] Sin[psi[t]] Derivative[1][psi][t])

In[46]:= eq1 = M L/
     2 (alpha''[t] Cos[alpha[t]] Cos[psi[t]] - 
      psi''[t] Sin[alpha[t]] Sin[psi[t]] - (alpha'[t]^2 + psi'[t]^2) Sin[
        alpha[t]] Cos[psi[t]] - 
      2 alpha'[t] psi'[t] Cos[alpha[t]] Cos[psi[t]]) == 
   Fx + (fs[t] + k Nf[t] theta'[t]) Sin[psi[t]];


eq2 = M L/
     2 (alpha''[t] Cos[alpha[t]] Sin[psi[t]] + 
      psi''[t] Sin[alpha[t]] Cos[psi[t]] - (alpha'[t]^2 + psi'[t]^2) Sin[
        alpha[t]] Sin[psi[t]] + 
      2 alpha'[t] psi'[t] Cos[alpha[t]] Cos[psi[t]]) == 
   Fy - (fs[t] + k Nf[t] theta'[t]) Cos[psi[t]];

eq3 = M R (alpha''[t] Cos[alpha[t]] - alpha'[t]^2 Sin[alpha[t]]) - 
    M L/2 (alpha''[t] Sin[alpha[t]] + alpha'[t]^2 Cos[alpha[t]]) == 
   Fz + M g + Nf[t];

One algebraic constraint for two unknown functions  Nf and fs, once theta is known


eq4 = theta''[t] - 2 k Nf[t]/(M R) theta'[t] - 2 fs[t]/(M R) == 0;

One  algebraic equation for alpha, once theta and psi are known


eq5 = L/2 Sin[alpha[t]] psi'[t] - R theta'[t] == 0;

Any differential equation  for Nf  and fs  is missing ,  so the given values at one point are void, the system is indefinite with respect to Nf and fs. On the other hand  eq5 is a constraint to the second order system {alpha, psi, theta} . So the system is over-under-determined, an error that NDSolve seems to be unable to detect.

In[37]:= soln = nDSolve[{eq1, eq2, eq3, eq4, eq5, alpha[0] == Pi/2, alpha'[0] == 1, 
   theta[0] == 0, theta'[0] == 1/2, psi[0] == 0, psi'[0] == 1, 
   Nf[0] == -Fz0 - M g, Nf'[0] == 1, fs[0] == 0, fs'[0] == 1}, {alpha[t], 
   theta[t], psi[t], Nf[t], fs[t]}, {t, 0, 10}]

--

Roland Franzius

Hello Roland. I put initial condition for theta dot. But still not working.

Attachments:

pre_iypt_2 (1).nb