Double Pendulum returning nan results. Why?

Question

When I run this code the double pendulum shows for a split second and then disappears. It seems most of the variables in the formula, num1, num2, etc., start returning nan after that split second when the pendulum disappears. Can anyone see why? I thought it had to do with converting degrees to radians at 180 degrees so I added the if statement in my convert function.

Here is the link to the Double Pendulum motion formula. I broke down the formula into pieces for simplicity's sake. https://www.myphysicslab.com/pendulum/double-pendulum-en.html

This idea came from Coding Train's Double Pendulum video: https://www.youtube.com/watch?v=uWzPe_S-RVE

#undef __STRICT_ANSI__
#include "window.h"
#include <SFML/Graphics.hpp>
#include <GL/glew.h>
#include <cmath>
#include <iostream>

namespace Window
{
    sf::ContextSettings settings;

    sf::RenderWindow window(sf::VideoMode(600, 600), "Window", sf::Style::Close | sf::Style::Titlebar, settings);

    void init()
    {
        settings.depthBits = 24;
        settings.majorVersion = 4;
        settings.minorVersion = 6;  //OpenGL 4.6

        glewInit();
        glViewport(0,0, 600, 600);
    }

    void close()
    {
        window.close();
    }

    void update()
    {
        window.display();
    }

    void checkForClose()
    {
        sf::Event windowEvent;
        while (window.pollEvent(windowEvent))
        {
            if (windowEvent.type == sf::Event::Closed)
            {
                close();
            }
        }
    }

    bool isOpen()
    {
        return window.isOpen();
    }

    double convDeg(double deg)
    {
        if (deg == 180.0)
        {
            return 0.0;
        }

        else
        {
            return deg * M_PI / 180.0;
        }
    }

    void runLoop()
    {
        double r1 = 100;
        double m1 = 40;
        double a1 = 90;
        double a1_v = 0;

        double r2 = 100;
        double m2 = 40;
        double a2 = 30;
        double a2_v = 0;

        double x1 = 0;
        double x2 = 0;
        double y1 = 0;
        double y2 = 0;

        double g = 0.005;

        double num1 = 0;
        double num2 = 0;
        double num3 = 0;
        double num4 = 0;
        double den = 0;

        double a1_a = 0;
        double a2_a = 0;

        while (window.isOpen())
        {
            num1 = -g * (2 * m1 + m2) * std::sin(convDeg(a1));
            num2 = -m2 * g * std::sin(convDeg(a1-2*a2));
            num3 = -2*std::sin(convDeg(a1-a2))*m2;
            num4 = a2_v*a2_v*r2+a1_v*a1_v*r1*std::cos(convDeg(a1-a2));
            den = r1 * (2*m1+m2-m2*cos(convDeg(2*a1-2*a2)));
            a1_a = (num1 + num2 + num3*num4) / den;

            num1 = 2 * std::sin(convDeg(a1-a2));
            num2 = (a1_v*a1_v*r1*(m1+m2));
            num3 = g * (m1 + m2) * std::cos(convDeg(a1));
            num4 = a2_v*a2_v*r2*m2*std::cos(convDeg(a1-a2));
            den = r2 * (2*m1+m2-m2*std::cos(convDeg(2*a1-2*a2)));
            a2_a = (num1*(num2+num3+num4)) / den;

            x1 = 300 + (r1 * std::sin((convDeg(a1))));
            y1 = 300 + (r1 * std::cos((convDeg(a1))));
            x2 = x1 + (r2 * std::sin((convDeg(a2))));
            y2 = y1 + (r2 * std::cos((convDeg(a2))));

            a1_v += a1_a;
            a2_v += a2_a;
            a1 += a1_v;
            a2 += a2_v;

            sf::Vertex pOne[] =
            {
                sf::Vertex(sf::Vector2f(300,300)),
                sf::Vertex(sf::Vector2f(x1,y1)),
            };

            sf::Vertex pTwo[] =
            {
                sf::Vertex(sf::Vector2f(x1,y1)),
                sf::Vertex(sf::Vector2f(x2,y2)),
            };

            window.clear();
            window.draw(pOne, 2, sf::Lines);
            window.draw(pTwo, 2, sf::Lines);

            update();
            checkForClose();

            std::cout << std::sin(convDeg(a1)) << "\n";

        }
    }
}

What exactly are you asking? Please consider https://stackoverflow.com/help/how-to-ask — Adrian W, May 26 '18 at 16:43
and _what_ is returning `nan`? Please be precise in what you ask. — Adrian W, May 26 '18 at 16:50
Tip: you can enable floating point exceptions to pin down the exact point where a nan is generated. — Matteo Italia, May 26 '18 at 17:23

meowgoesthedog · Answer 1 · 2018-05-26T17:58:40.310

0

List of problems:

convDeg: 180 and 0 degrees are not equivalent.
Angular variables: acceleration / velocity a1_a, a1_v etc are in radians, whereas your angle variables a1, a2 etc are in degrees. Every time the loop runs, your code wrongly converts degrees into radians. Do the conversion before the main loop.
(Speculation based on a previous answer - I could be wrong) The formulas on MyPhysicsLab seem to contradict those from other sources, e.g. here and here; a quick calculation using Lagrangian mechanics agrees with these two sources and not MPL.

edited May 26 '18 at 17:58

answered May 26 '18 at 17:19

meowgoesthedog

14,670
4
27
40

Where are 180 and 0 being equivalent? I have it returning 0 radians for 180 degrees. – Don May 26 '18 at 17:30
@Don why is 0 radians equivalent to 180 degrees? Did you mean 360? – meowgoesthedog May 26 '18 at 17:31
I guess I'm just really confused. When converting 180 degrees to radians and asking for the sin, it is returning a really really small number instead of just 0. That was my goal with my if statement. I assume I'm doing it wrong. Haha – Don May 26 '18 at 17:34
@Don quick tip - you can use `fmod` to "normalize" the angle to [0, 360) degrees or [0, 2pi) radians. I gather that's (the gist of) what you were trying to implement with the if statement in `convDeg`. – meowgoesthedog May 26 '18 at 17:36
I apologize. I'm new to C++. How would I do that? – Don May 26 '18 at 17:38
@Don something like `fmod(deg, 360.0) * (M_PI / 180.0)`. – meowgoesthedog May 26 '18 at 17:40
@Don the equations themselves are in radians because that is the natural unit of angle; you should convert your input angles `a1, a2` etc into radians *before* running the rest of the code. – meowgoesthedog May 26 '18 at 17:51
Alright. I'll try that. Just confused why it worked for the guy in the youtube video but not for me? – Don May 26 '18 at 17:53
@Don skip to [8:42](https://youtu.be/uWzPe_S-RVE?t=522) in that video - his input angles `a1, a2` were **already** in radians. – meowgoesthedog May 26 '18 at 17:54
Wow. Didn't notice. So I made some changes and it looks closer to the expected result. Should I edit my original post with my new code? – Don May 26 '18 at 18:09
@Don yes but don't overwrite the original incorrect code; create a new subsection – meowgoesthedog May 26 '18 at 18:10

Double Pendulum returning nan results. Why?

1 Answers1