Двойной маятник, возвращающий результаты НАН.Зачем? - PullRequest
Новая
       26

Двойной маятник, возвращающий результаты НАН.Зачем?

0 голосов
/ 26 мая 2018

Когда я запускаю этот код, двойной маятник показывает на доли секунды, а затем исчезает.Кажется, что большинство переменных в формуле, num1, num2 и т. Д., Начинают возвращать nan после той доли секунды, когда маятник исчезает.Кто-нибудь может понять почему?Я подумал, что это связано с преобразованием градусов в радианы на 180 градусов, поэтому я добавил оператор if в свою функцию преобразования.

Вот ссылка на формулу движения двойного маятника.Для простоты я разбил формулу на кусочки.https://www.myphysicslab.com/pendulum/double-pendulum-en.html

Эта идея пришла из видео Двойного маятника Coding Train: 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";

        }
    }
}

1 Ответ

0 голосов
/ 26 мая 2018

Список проблем:

  • convDeg: 180 и 0 градусов не эквивалентны.
  • Угловые переменные: ускорение / скорость a1_a, a1_v и т. Д. В радиан, тогда как ваши угловые переменные a1, a2 и т.д. находятся в градусах .Каждый раз, когда запускается цикл, ваш код неправильно конвертирует градусы в радианы.Выполните преобразование до основного цикла.
  • (предположение, основанное на предыдущем ответе - я могу ошибаться) Формулы на MyPhysicsLab, похоже, противоречат формулам из других источников, например, здесь и здесь ;быстрый расчет с использованием лагранжевой механики согласуется с этими двумя источниками, а не MPL.
...