In[]:=
k=266*1000;m=3755;L=2;r=0.5;A=L^2/4+L*x[t]+x[t]^2+y[t]^2-0.5*L^2*Cos[theta[t]]-L*x[t]*Cos[theta[t]]+0.25*L^2*Cos[theta[t]]^2-L*y[t]*Sin[theta[t]]+0.25*L^2*Sin[theta[t]]^2;B=L^2/4-L*x[t]+x[t]^2+y[t]^2-0.5*L^2*Cos[theta[t]]^2+L*x[t]*Cos[theta[t]]+0.25*L^2*Cos[theta[t]]^2+L*y[t]*Sin[theta[t]]+0.25*L^2*Sin[theta[t]]^2;oneA=0.5*k*L^2*Sin[theta[t]]+0.25*k*L^2*(-2*Cos[theta[t]]*Sin[theta[t]])+0.25*k*L*(2*Sin[theta[t]]*Cos[theta[t]]);oneB=-k*r*(0.5*L^2*Sin[theta[t]]+L*x[t]*Sin[theta[t]]+0.25*L^2*(-2*Cos[theta[t]]*Sin[theta[t]])-L*y[t]*Cos[theta[t]]+0.25*L^2*(2*Sin[theta[t]]*Cos[theta[t]]))/(2*Sqrt[A]);oneC=-k*r*(0.5*L^2*Sin[theta[t]]-L*x[t]*Sin[theta[t]]+0.25*L^2*(2*Cos[theta[t]]*Sin[theta[t]])+L*y[t]*Cos[theta[t]]+0.25*L^2*(2*Sin[theta[t]]*Cos[theta[t]]))/(2*Sqrt[B]);ONE=oneA+oneB+oneC;two=2*k*x[t]-k*r*(L+2*x[t]-L*Cos[theta[t]])/(2*Sqrt[A])-k*r*(-L+2*x[t]+L*Cos[theta[t]])/(2*Sqrt[B]);three=-9.81*m+2*k*y[t]-k*r*(2*y[t]-L*Sin[theta[t]])/(2*Sqrt[A])-k*r*(2*y[t]+L*Sin[theta[t]])/(2*Sqrt[B]);S=NDSolve[{x''[t]==-two/(m),y''[t]==-three/(m),theta''[t]==-ONE*4/(m*L^2),x[0]==0.0522,y[0]==0.6,theta[0]==0.0349,x'[0]==0,y'[0]==0,theta'[0]==0},{x,y,theta},{t,0,20}]
Out[]=
xInterpolatingFunction,yInterpolatingFunction,thetaInterpolatingFunction
In[]:=
Plot[{x[t],y[t],theta[t]}/.S,{t,0,20},PlotLegends->{"x","y","Theta"}]Plot[{x'[t],y'[t],theta'[t]}/.S,{t,0,20},PlotLegends->{"x","y","Theta"}]Plot[{x''[t],y''[t],theta''[t]}/.S,{t,0,20},PlotLegends->{"x","y","Theta"}]
Out[]=
Out[]=
Out[]=