In[]:=
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=-rz0rx0​​Fz0=B/r0^5(rx0-2rz0+5A0/r0^2rz0)​​rx=L/2Sin[alpha[t]]Cos[psi[t]]​​ry=L/2Sin[alpha[t]]Sin[psi[t]]​​rz=L/2Cos[alpha[t]]+d+l/2​​r=Sqrt[rx^2+ry^2+rz^2]​​A=(rx^2+ry^2)Cos[alpha[t]]-ryrzSin[alpha[t]]Sin[psi[t]]-rzrxSin[alpha[t]]Cos[psi[t]]​​Fx=B/r^5(-4rxCos[alpha[t]]+rzSin[alpha[t]]Cos[psi[t]]+5A/r^2rx)​​Fy=B/r^5(-4ryCos[alpha[t]]+rzSin[alpha[t]]Sin[psi[t]]+5A/r^2ry)​​Fz=B/r^5(rySin[alpha[t]]Sin[psi[t]]+rxSin[alpha[t]]Cos[psi[t]]-2rzCos[alpha[t]]+5A/r^2rz)​​eq1=ML/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]]-2alpha'[t]psi'[t]Cos[alpha[t]]Cos[psi[t]])==Fx+(fs[t]+kNf[t]theta'[t])Sin[psi[t]]​​eq2=ML/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]]+2alpha'[t]psi'[t]Cos[alpha[t]]Cos[psi[t]])==Fy-(fs[t]+kNf[t]theta'[t])Cos[psi[t]]​​eq3=MR(alpha''[t]Cos[alpha[t]]-alpha'[t]^2Sin[alpha[t]])-ML/2(alpha''[t]Sin[alpha[t]]+alpha'[t]^2Cos[alpha[t]])==Fz+Mg+Nf[t]​​eq4=theta''[t]-2kNf[t]/(MR)theta'[t]-2fs[t]/(MR)==0​​eq5=L/2Sin[alpha[t]]psi'[t]-Rtheta'[t]==0​​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-Mg,Nf'[0]==1,fs[0]==0,fs'[0]==1},{alpha[t],theta[t],psi[t],Nf[t],fs[t]},{t,0,10}]​​Plot[psi[t]/.soln,{t,0,10}]​​​​
Out[]=
0.001
Out[]=
0.001
Out[]=
0.001
Out[]=
0.1
Out[]=
0.001
Out[]=
1
Out[]=
0.001
Out[]=
10
Out[]=
0.0005
Out[]=
0.002
Out[]=
0.00206155
Out[]=
-1.×
-6
10
Out[]=
-1.57182×
10
10
Out[]=
0.0005Cos[psi[t]]Sin[alpha[t]]
Out[]=
0.0005Sin[alpha[t]]Sin[psi[t]]
Out[]=
0.0015+0.0005Cos[alpha[t]]
Out[]=
2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
Out[]=
-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Cos[psi[t]]
2
Sin[alpha[t]]
-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Sin[alpha[t]]
2
Sin[psi[t]]
+Cos[alpha[t]](2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
)
Out[]=
0.1(0.0015+0.0005Cos[alpha[t]])Cos[psi[t]]Sin[alpha[t]]-0.002Cos[alpha[t]]Cos[psi[t]]Sin[alpha[t]]+
0.0025Cos[psi[t]]Sin[alpha[t]]-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Cos[psi[t]]
2
Sin[alpha[t]]
-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Sin[alpha[t]]
2
Sin[psi[t]]
+Cos[alpha[t]]2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]

2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
5/2
(
2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
)
Out[]=
0.1(0.0015+0.0005Cos[alpha[t]])Sin[alpha[t]]Sin[psi[t]]-0.002Cos[alpha[t]]Sin[alpha[t]]Sin[psi[t]]+
0.0025Sin[alpha[t]]Sin[psi[t]]-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Cos[psi[t]]
2
Sin[alpha[t]]
-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Sin[alpha[t]]
2
Sin[psi[t]]
+Cos[alpha[t]]2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]

2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
5/2
(
2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
)
Out[]=
0.1-2(0.0015+0.0005Cos[alpha[t]])Cos[alpha[t]]+0.0005
2
Cos[psi[t]]
2
Sin[alpha[t]]
+0.0005
2
Sin[alpha[t]]
2
Sin[psi[t]]
+
5(0.0015+0.0005Cos[alpha[t]])-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Cos[psi[t]]
2
Sin[alpha[t]]
-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Sin[alpha[t]]
2
Sin[psi[t]]
+Cos[alpha[t]]2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]

2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
5/2
(
2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
)
Out[]=
5.×
-7
10
(-2Cos[alpha[t]]Cos[psi[t]]
′
alpha
[t]
′
psi
[t]-Cos[psi[t]]Sin[alpha[t]](
2
′
alpha
[t]
+
2
′
psi
[t]
)+Cos[alpha[t]]Cos[psi[t]]
′′
alpha
[t]-Sin[alpha[t]]Sin[psi[t]]
′′
psi
[t])
0.1(0.0015+0.0005Cos[alpha[t]])Cos[psi[t]]Sin[alpha[t]]-0.002Cos[alpha[t]]Cos[psi[t]]Sin[alpha[t]]+
0.0025Cos[psi[t]]Sin[alpha[t]]-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Cos[psi[t]]
2
Sin[alpha[t]]
-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Sin[alpha[t]]
2
Sin[psi[t]]
+Cos[alpha[t]]2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]

2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
5/2
(
2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
)
+Sin[psi[t]](fs[t]+Nf[t]
′
theta
[t])
Out[]=
5.×
-7
10
(2Cos[alpha[t]]Cos[psi[t]]
′
alpha
[t]
′
psi
[t]-Sin[alpha[t]]Sin[psi[t]](
2
′
alpha
[t]
+
2
′
psi
[t]
)+Cos[alpha[t]]Sin[psi[t]]
′′
alpha
[t]+Cos[psi[t]]Sin[alpha[t]]
′′
psi
[t])
0.1(0.0015+0.0005Cos[alpha[t]])Sin[alpha[t]]Sin[psi[t]]-0.002Cos[alpha[t]]Sin[alpha[t]]Sin[psi[t]]+
0.0025Sin[alpha[t]]Sin[psi[t]]-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Cos[psi[t]]
2
Sin[alpha[t]]
-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Sin[alpha[t]]
2
Sin[psi[t]]
+Cos[alpha[t]]2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]

2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
5/2
(
2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
)
-Cos[psi[t]](fs[t]+Nf[t]
′
theta
[t])
Out[]=
1.×
-6
10
(-Sin[alpha[t]]
2
′
alpha
[t]
+Cos[alpha[t]]
′′
alpha
[t])-5.×
-7
10
(Cos[alpha[t]]
2
′
alpha
[t]
+Sin[alpha[t]]
′′
alpha
[t])0.01+Nf[t]+
0.1-2(0.0015+0.0005Cos[alpha[t]])Cos[alpha[t]]+0.0005
2
Cos[psi[t]]
2
Sin[alpha[t]]
+0.0005
2
Sin[alpha[t]]
2
Sin[psi[t]]
+
5(0.0015+0.0005Cos[alpha[t]])-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Cos[psi[t]]
2
Sin[alpha[t]]
-0.0005(0.0015+0.0005Cos[alpha[t]])
2
Sin[alpha[t]]
2
Sin[psi[t]]
+Cos[alpha[t]]2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]

2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
5/2
(
2
(0.0015+0.0005Cos[alpha[t]])
+2.5×
-7
10
2
Cos[psi[t]]
2
Sin[alpha[t]]
+2.5×
-7
10
2
Sin[alpha[t]]
2
Sin[psi[t]]
)
Out[]=
-2.×
6
10
fs[t]-2.×
6
10
Nf[t]
′
theta
[t]+
′′
theta
[t]0
Out[]=
0.0005Sin[alpha[t]]
′
psi
[t]-0.001
′
theta
[t]0
NDSolve
:Some of the functions have zero differential order, so the equations will be solved as a system of differential-algebraic equations.
NDSolve
:Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.
Out[]=
{}
Out[]=
2
4
6
8
10
-1.0
-0.5
0.5
1.0
​