2
$\begingroup$

I'm trying to replicate the graphs and animation found in this page that studies the movement of a mass over a dome solving numerically the differential equation r*θ''[t] == g*Sin[θ[t]] but Mathematica gets a wrong result

z = NDSolve[{r*θ''[t] == g*Sin[θ[t]], 
             θ'[0] == ω0, θ[0] == θ0}, θ, {t, 0, 10}]
Plot[Evaluate[θ[t] /. z], Evaluate[Flatten[{t, θ["Domain"] /. z}]]]

even if I change Method of solution in NDSolve.

As you can see, the MATLAB used there gives an adequate model of the situation. I used the code here trying to replicate MATLAB's code but I get the same wrong solution.

DOPRIamat = {{1/5}, {3/40, 9/40}, {44/45, -56/15, 32/9}, 
             {19372/6561, -25360/2187, 64448/6561, -212/729},
             {9017/3168, -355/33, 46732/5247, 49/176, -5103/18656}, 
             {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84}};
DOPRIbvec = {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0};
DOPRIcvec = {1/5, 3/10, 4/5, 8/9, 1, 1};
DOPRIevec = {71/57600, 0, -71/16695, 71/1920, -17253/339200, 22/525, -1/40};
DOPRICoefficients[5, p_] := N[{DOPRIamat, DOPRIbvec, DOPRIcvec, DOPRIevec}, p];

l := NDSolve[
      {r*θ''[t] == g*Sin[θ[t]], θ'[0] == ω0, θ[0] == θ0, 
       WhenEvent[θ[t] >= Pi/4, "StopIntegration"]}, 
      θ, {t, 0, 10}, 
      Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5, 
                 "Coefficients" -> DOPRICoefficients,
                 "StiffnessTest" -> False}
      ]

Plot[
  Evaluate[θ[t] /. l], Evaluate[Flatten[{t, θ["Domain"] /. l}]], 
  AxesLabel -> {"t en s", "θ en rad"}]
$\endgroup$
1
  • 3
    $\begingroup$ Please edit your question to include the Mathematica code you tried as appropriately formatted text, together with any definitions of symbols you used. $\endgroup$ – MarcoB Jun 4 '20 at 20:46
4
$\begingroup$

With your first attempt, there are a number of warnings issued due to $r, g, \omega 0, \theta 0$ being undefined. NDSolve is a numerical solver, and must be able to resolve the functions to numerical values at all steps.

Looking at the page you showed, they define $r = 1$, $g = 9.8$, and I'm not sure about these next two but it looks like $\theta 0 = 0.01$ and $\omega 0 = \sqrt{9.8/1}* 2* \sin(0.01/2) = 0.0313048$. If we provide these values to NDSolve, it has no issues solving it. I'm pretty confident in Mathematica's capabilities, particularly for such a simple equation, and I would be very surprised to find out that it's wrong.

r = 1;
g = 9.8;
θ0 = 0.01;
ω0 = Sqrt[9.8/r]*2*Sin[0.01/2];
sol = θ /. First@NDSolve[{
     r θ''[t] == g Sin[θ[t]],
     θ[0] == θ0,
     θ'[0] == ω0
     },
    θ,
    {t, 0, 10}
    ]
Plot[
 sol[t]/Degree,
 {t, 0, 1.5},
 AxesLabel -> {"Time (s)", "θ (Degrees)"}
]

Plot of angle versus time.

As far as I can tell, this graph agrees exactly with the one on the page you linked. At $t=1$, I get 13.0981 degrees, which looks pretty much like what the Matlab graph shows.

$\endgroup$
4
  • 1
    $\begingroup$ Thanks so much, I forgot setting values for those variables here. Also, I would like to know what's the difference between using /.First@ instead of replacing later $\endgroup$ – Edgar Castro Jun 4 '20 at 22:03
  • $\begingroup$ @EdgarCastro There isn't really any difference most of the time, I just like to do it when I'm getting the solution to get it out of the way early. That way I can just use sol directly without having to perform the replacement all the time. Sometimes I find it more elegant to do the replacement during the plotting, like I did in this question, though that was mostly because I was using Cartesian coordinates for plotting. $\endgroup$ – MassDefect Jun 4 '20 at 22:15
  • 1
    $\begingroup$ You can also use NSolveValue to get the solution directly. $\endgroup$ – tad Jun 5 '20 at 4:01
  • 1
    $\begingroup$ @tad NDSolveValue, but yes, you're right. In this case, that's probably the easiest way to do it. $\endgroup$ – MassDefect Jun 5 '20 at 4:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.