Hi Guys
Please refer to the below code. The system being presented two elastic pendulums connected by a rigid rod. I have represented the system using Lagrangian mechanics with a holonomic constraint to develop the required behavior.
I have used the holonomic constraint to represent a variable as a function of other variables. The time derivative is also required, thus all the equations which are presented in the beginning.
The (b+r) term presented as "B" is made the subject of the formula and then the derivative also taken. When I run the code I am attaining an error with regards to a non-numerical value of a derivative being 0. Can anyone assist with regards to what my problem is? Further to this I would appreciate any advice with regards to issues with the manner in which I have coded? I come from a Matlab and Python mindset.
r = 1;
L = 2;
m1 = 3500;
m2 = 3500;
k = 353160;
A = a'[t] + r;
R1 = 2*A*Cos[theta1[t]]*Cos[theta2[t]];
one = -A^2 + L^2;
two = A^2*Cos[2*theta1[t] - 2*theta2[t]];
three = -L^2*Cos[2*theta2[t]];
four = -2*A*L*Sin[theta1[t]];
five = -2*A*L*Sin[theta1[t] - 2*theta2[t]];
R2 = Sqrt[4]*Sqrt[one + two + three + four + five];
six = 2*L*Sin[theta2[t]];
seven = 2*A*Sin[theta1[t]]*Sin[theta2[t]];
R3 = six + seven;
B = R1 + R2 + R3;
R1dot = 2*a'[t]*Cos[theta1[t]]*Cos[theta2[t]] -
2*A*theta1'[t]*Sin[theta1[t]]*Cos[theta2[t]] -
2*A*theta2'[t]*Cos[theta1[t]]*Sin[theta2[t]];
oned = -2*a'[t]*A;
twod = 2*a'[t]*A*Cos[2*theta1[t]]*Cos[2*theta2[t]] -
A^2*2*theta1'[t]*Sin[2*theta1[t]]*Cos[2*theta2[t]] -
A^2*Cos[2*theta2[t]]*2*theta2'[t]*Sin[2*theta2[t]] +
2*a'[t]*A*Sin[2*theta1[t]]*Sin[2*theta2[t]] +
A^2*2*theta1'[t]*Cos[2*theta1[t]]*Sin[2*theta2[t]] +
A^2*Sin[2*theta1[t]]*2*theta2'[t]*Cos[2*theta2[t]];
threed = L^2*s*theta2'[t]*Sin[2*theta2[t]];
fourd = -2*a'[t]*L*Sin[theta1[t]] - 2*A*L*theta1'[t]*Cos[theta1[t]];
fived = -2*a'[t]*L*Sin[theta1[t]]*Cos[2*theta2[t]] -
2*A*L*theta1'[t]*Cos[theta1[t]]*Cos[2*theta2[t]] -
2*A*L*Sin[theta1[t]]*2*theta1'[t]*Sin[2*theta2[t]] +
2 a'[t]*L*Cos[theta1[t]]*Sin[2*theta2[t]] -
2*A*L*theta1'[t]*Sin[theta1[t]]*Sin[2*theta2[t]] +
2*A*L*Cos[theta1[t]]*2*theta2'[t]*Cos[2*theta2[t]];
sixd = 2*L*theta2'[t]*Cos[theta2[t]];
sevend =
2*(a'[t]*Sin[theta1[t]]*Sin[theta2[t]] +
A*theta1'[t]*Cos[theta1[t]]*Sin[theta2[t]] +
A*Sin[theta1[t]]*theta2'[t]*Cos[theta2[t]]);
R2dot = (oned + twod + threed + fourd + fived)/(2*
Sqrt[one + two + three + four + five]);
R3dot = sixd + sevend;
bdot = (0.5*R1dot + R2dot + R3dot)/(4*Sqrt[R1 + R2 + R3]);
NDSolve[{b''[
t] == (m2*theta2'[t]^2*B - k*(B - r) + 9.81*m2*Cos[theta2[t]])/
m2, a''[t] == (m1*theta1'[t]^2*A - k*a[t] +
9.81*m1*Cos[theta1[t]])/m2,
theta1''[t] == (-2*m1*theta1'[t]*a'[t]*A -
9.81*m1*Sin[theta1[t]]*A)/(m1*A^2),
theta2''[t] == (-2*m2*theta2'[t]*bdot*B -
9.81*m2*Sin[theta2[t]]*B)/(m2*B^2), a[0] == 0, a'[0] == 0,
b[0] == 0, b'[0] == 0, theta1[0] == 0, theta1'[0] == 0,
theta2[0] == 0, theta2'[0] == 0}, {a, b, theta1, theta2}, {t, 0,
10}]