How to handle NDSolve::ndsz problem (singularity problem) - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-08-18T01:44:41Z https://mathematica.stackexchange.com/feeds/question/9989 http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://mathematica.stackexchange.com/q/9989 6 How to handle NDSolve::ndsz problem (singularity problem) Nasser https://mathematica.stackexchange.com/users/70 2012-08-28T20:49:19Z 2014-11-22T07:33:03Z <p>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. </p> <p>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 <code>NDSolve::ndsz</code> error</p> <pre><code>NDSolve::ndsz: At t == 0.09571750107721341`, step size is effectively zero; singularity or stiff system suspected </code></pre> <p>Here is an example, where a tiny change in initial conditions causes this error:</p> <pre><code>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; </code></pre> <p>Now depending on the IC, this problem shows up. This one gives the above <code>ndsz</code> error </p> <pre><code>ic={x==0.23,x'==-2,theta==5,theta'==3}; NDSolve[Flatten[{eq1,eq2,ic}],{x,theta,x',theta'},{t,0,0.1}, Method-&gt;{"StiffnessSwitching"}] </code></pre> <p>While if I changed <code>x</code> above to be <code>0.24</code> instead of <code>0.23</code>, now it works</p> <pre><code>ic={x==0.24,x'==-2,theta==5,theta'==3}; sol=NDSolve[Flatten[{eq1,eq2,ic}],{x,theta,x',theta'},{t,0,0.1}, Method-&gt;{"StiffnessSwitching"}] </code></pre> <p>I tried many of the tricks for stiffness gives here <a href="http://reference.wolfram.com/mathematica/tutorial/NDSolveOverview.html">http://reference.wolfram.com/mathematica/tutorial/NDSolveOverview.html</a> but none has helped so far. </p> <p>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. </p> <p>The error seems to say it is singularity. Either way, I tried all the tricks I know about relating to stiffness methods. </p> <p>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?</p> <p>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 <code>NDSolve</code>? Looking to see if there are options for <code>NDSOlve</code> 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.</p> <p>Here is the CDF itself in case you want to see the actual demo <a href="http://12000.org/my_notes/mma_demos/solid_pendulum_with_spring_mass/">http://12000.org/my_notes/mma_demos/solid_pendulum_with_spring_mass/</a> </p> <p><strong>update:</strong> 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 <code>report.nb</code> which shows the derivations using the Lagrangian method. The case described here is in <code>report.nb</code>, 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.</p> <p>thanks</p> https://mathematica.stackexchange.com/questions/9989/-/10027#10027 2 Answer by Dr. belisarius for How to handle NDSolve::ndsz problem (singularity problem) Dr. belisarius https://mathematica.stackexchange.com/users/193 2012-08-29T14:36:06Z 2012-08-29T15:06:31Z <p>I haven't checked the physical model, but you are getting a very hard rebound. </p> <pre><code>ic = {x == 20/100, x' == -18610/10000, theta == Mod[5, 2 Pi], theta' == 3}; s = NDSolve[ Flatten[{eq1, eq2, ic}], {x, theta, x', theta'}, {t, 0, .2}, Method -&gt; {"ExplicitRungeKutta", "StiffnessTest" -&gt; False}]; {Plot[x'[t] /. s, {t, 0, .2}], Plot[theta'[t] /. s, {t, 0, .2}]} </code></pre> <p><img src="https://i.stack.imgur.com/1Y15Z.png" alt="Mathematica graphics"></p> <p>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.</p> <p><strong>Edit</strong></p> <p>Stiffness progression of speed:</p> <p><img src="https://i.stack.imgur.com/JmF1T.png" alt="Mathematica graphics"></p> https://mathematica.stackexchange.com/questions/9989/-/10028#10028 1 Answer by user21 for How to handle NDSolve::ndsz problem (singularity problem) user21 https://mathematica.stackexchange.com/users/0 2012-08-29T14:39:14Z 2012-08-29T14:39:14Z <p>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 </p> <pre><code>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 == a, x' == -2, theta == 5, theta' == 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[x][t], Derivative[theta][t]} /. res], {t, 0, 0.1}] , {a, 0.20, 0.3}] </code></pre> <p>You'll see how <code>theta</code> and <code>theta'</code> are going south.</p>