How to handle NDSolve::ndsz problem (singularity problem)

Question

I have 2 second order differential equations (non-linear). The physics behind them is correct. I verified the equations many times. It is a solid pendulum with a mass-spring at the end of it.

Now, assuming the mathematical model is correct and that the 2 differential equations are derived correctly based on the physics, I found that under some specific initial conditions, I get the dreaded NDSolve::ndsz error

NDSolve::ndsz: At t == 0.09571750107721341`, step size is effectively zero; 
singularity or stiff system suspected

Here is an example, where a tiny change in initial conditions causes this error:

Clear[theta,x,t];
eq1 = 98  Sin[theta[t]]+x[t]-10 x[t] theta'[t]^2+10 (x''[t]+0.5  theta''[t])==0;
eq2 =-98 (-0.5  Sin[theta[t]]-Cos[theta[t]] x[t])+5 (4 x[t] x'[t] theta'[t]+ 
        (x''[t]+0.5 theta''[t])+2 x[t]^2 theta''[t])==0;

Now depending on the IC, this problem shows up. This one gives the above ndsz error

ic={x[0]==0.23,x'[0]==-2,theta[0]==5,theta'[0]==3};
NDSolve[Flatten[{eq1,eq2,ic}],{x,theta,x',theta'},{t,0,0.1},
       Method->{"StiffnessSwitching"}]

While if I changed x[0] above to be 0.24 instead of 0.23, now it works

ic={x[0]==0.24,x'[0]==-2,theta[0]==5,theta'[0]==3};
sol=NDSolve[Flatten[{eq1,eq2,ic}],{x,theta,x',theta'},{t,0,0.1},
            Method->{"StiffnessSwitching"}]

I tried many of the tricks for stiffness gives here http://reference.wolfram.com/mathematica/tutorial/NDSolveOverview.html but none has helped so far.

I am actually not sure if it is stiffness issue or singularity issue. Singularity means to me that equations are wrong, right? and that should not really happen if the physics is correct, or may be I do not understand exactly what this means.

The error seems to say it is singularity. Either way, I tried all the tricks I know about relating to stiffness methods.

My question is really this: Assuming the equations are correct (and I am almost sure they are). How does one handle such a problem? If I trap this error, what should I do next?

Could a valid and correct differential equations (i.e. derived from a correct mathematical model) generate a singularity error for specific initial conditions? and if so, how to handle this using NDSolve? Looking to see if there are options for NDSOlve that I overlooked to help handle this. Feel free to try any other options above to see if you can get rid of this error for this specific case of initial conditions.

Here is the CDF itself in case you want to see the actual demo http://12000.org/my_notes/mma_demos/solid_pendulum_with_spring_mass/

update: fyi, I just uploaded to the above link a small report showing the derivation of the differential equations from the model. The links are there for the report.pdf and html version, as well as the Mathematica notebook itself report.nb which shows the derivations using the Lagrangian method. The case described here is in report.nb, the last one at the end. This problem again shows up only under very specific initial conditions. This problem shows up when I model the pendulum arm as mass-less instead of solid pendulum. I think I did the model correctly in this case as well. So, right now, for the demo itself, I force the pendulum to be a solid pendulum only to avoid this problem.

thanks

Answer 1

I haven't checked the physical model, but you are getting a very hard rebound.

ic = {x[0] == 20/100, x'[0] == -18610/10000, theta[0] == Mod[5, 2 Pi], theta'[0] == 3};
s = NDSolve[
   Flatten[{eq1, eq2, ic}], {x, theta, x', theta'}, {t, 0, .2}, 
   Method -> {"ExplicitRungeKutta", "StiffnessTest" -> False}];

{Plot[x'[t] /. s, {t, 0, .2}], Plot[theta'[t] /. s, {t, 0, .2}]}

You may try to verify whether that behavior is compatible with the physics. And if it is, you may try to partition your problem at the stiff point and get new initial conditions by taking approximate limits around it.

Edit

Stiffness progression of speed:

Answer 2

I think something like this is expected for a massless Pendulum (I did not look at the physics), so what you are seeing then is a singularity. If you look at

Clear[theta, x, t];
eq1 = 98 Sin[theta[t]] + x[t] - 10 x[t] theta'[t]^2 + 
    10 (x''[t] + 0.5 theta''[t]) == 0;
eq2 = -98 (-0.5 Sin[theta[t]] - Cos[theta[t]] x[t]) + 
    5 (4 x[t] x'[t] theta'[t] + (x''[t] + 0.5 theta''[t]) + 
       2 x[t]^2 theta''[t]) == 0;

Manipulate[
 ic = {x[0] == a, x'[0] == -2, theta[0] == 5, theta'[0] == 3};
 res = NDSolve[
    Flatten[{eq1, eq2, ic}], {x[t], theta[t], x'[t], theta'[t]}, {t, 
     0, 0.1}] // Flatten;
 Plot[Evaluate[{x[t], theta[t], Derivative[1][x][t], 
     Derivative[1][theta][t]} /. res], {t, 0, 0.1}]
 , {a, 0.20, 0.3}]

You'll see how theta and theta' are going south.