Message Boards Message Boards

Error solving differential equations: Infinite expression encountered?

Posted 2 years ago

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}]
POSTED BY: Mishal Mohanlal

It seems to me that your equation for theta2''[t] is singular when t=0 and theta1[0] == theta2[0] == 0.

POSTED BY: Gianluca Gorni
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract