이중 진자(double pendulum)

23.02.17 추가

이중진자 시뮬레이션의 새 코드가 나왔으므로 식 유도가 아니라 코드에 관심있는 사람들은 아래 글 참고 바람.

https://ellipsoid.tistory.com/entry/%EC%A1%B0%EA%B8%88-%EB%8D%94-%EA%B9%94%EB%81%94%ED%95%9C-%EC%9D%B4%EC%A4%91%EC%A7%84%EC%9E%90-%EC%8B%9C%EB%AE%AC%EB%A0%88%EC%9D%B4%EC%85%98

조금 더 깔끔한 이중진자 시뮬레이션

 

 

이중 진자 시뮬레이션을 C++과 SFML로 만들어 보았다.

SFML이 생각보다 불편해서 중간에 좀 헤매긴 했다.

 

#include <string>
#include <exception>
#include <iostream>
#include <SFML/Graphics.hpp>
#include <list>
#include <functional>

struct movedata
{
	double theta;
	double pi;

	double p1;
	double p2;
	
	movedata operator+(const movedata& ref)const
	{
		movedata temp = {theta + ref.theta , pi + ref.pi, p1 + ref.p1, p2 + ref.p2};
		return temp;
	}
};

movedata operator*(double frac, movedata t)
{
	movedata temp = {frac * t.theta, frac * t.pi, frac * t.p1, frac * t.p2};
	return temp;
}

movedata operator/(movedata t, double frac)
{
	movedata temp = {t.theta / frac, t.pi / frac, t.p1 / frac, t.p2 / frac};
	return temp;
}

class System : private movedata
{
public:
	movedata operator()(movedata st)
	{
		movedata temp;

		double u = m2 / m1;

		double A1 = (st.p1 * st.p2 * std::sin(st.theta - st.pi)) / (m1 * length * length * (1 + u * std::sin(st.theta - st.pi) * std::sin(st.theta - st.pi)));
		double A2 = (st.p1 * st.p1 * u - 2.0 * st.p1 * st.p2 * u * std::cos(st.theta -st.pi) + st.p2 * st.p2 * (1 + u)) * (std::sin(2 * (st.theta - st.pi))) / (2.0 * m1 * length * length * (1 + u * pow(std::sin(st.theta - st.pi), 2.0)));

		temp.theta = (st.p1 - (st.p2 * std::cos(st.theta-st.pi))) / (m1 * pow(length, 2.0) * (1 + u * pow(std::sin(st.theta - st.pi) ,2.0)));
		temp.pi = st.p2 * (1 + u) - (st.p1 * u * std::cos(st.theta - st.pi)) / (m1 * pow(length, 2) * (1 + u * pow(std::sin(st.theta - st.pi), 2.0)));
		temp.p1 = -m1 * (1 + u) * 9.806 * length * std::sin(st.theta) - A1 - A2;
		temp.p2 = -m1 * u * 9.806 * length * std::sin(st.pi) + A1 - A2;

		return temp;
	}

private:
	double length = 1.0;

	double m1 = 1.0;
	double m2 = 1.0;
};

template <typename func, typename state>
class movement
{
private:
	std::list<state> move;	//시간, 각속도, 각도

	void RungeKutta(func& sys, int timestep, double stepsize)
	{
		state k1 = {0.0, }, k2 = { 0.0, }, k3 = { 0.0, }, k4 = { 0.0, };

		for (size_t i = 0; i < timestep; i++)
		{
			state temp = move.back();

			k1 = stepsize * sys(temp);

			k2 = stepsize * sys(temp + (0.5 * k1));

			k3 = stepsize * sys(temp + (0.5 * k2));

			k4 = stepsize * sys(temp + k3);

			temp = temp + ((k1 + (2.0 * k2) + (2.0 * k3) + k4) / 6.0);

			move.push_back(temp);
		}
	}

public:
	movement(double theta, double pi, int timestep, double stepsize)
	{
		state m = { theta, pi, 0.0, 0.0 };
		move.push_back(m);

		func sys;
		RungeKutta(sys, timestep, stepsize);
	};

	void calcsteps(int timestep, double stepsize)
	{
		func sys;
		RungeKutta(sys, timestep, stepsize);
	}

	state getcurrentstep()
	{
		state temp = move.front();
		return temp;
	}

	int stepnum()
	{
		return move.size();
	}

	void deletesteps(int number)
	{
		for (size_t i = 0; i < number; i++)
			move.pop_front();
	}
};

template <typename T>
double deg2rad(T theta)
{
	return theta * 3.15159265358979 / 180.0;
}

template <typename T>
double rad2deg(T rad)
{
	return rad * 180.0 / 3.14159265358979;
}

int main(void)
{
	sf::RenderWindow window(sf::VideoMode(1200, 675), "SFML window");

	window.setVerticalSyncEnabled(true);

	float fps;
	sf::Clock clock = sf::Clock::Clock();
	sf::Time previousTime = clock.getElapsedTime();
	sf::Time currentTime;

	const double diff = 1.0 / 60.0;
	movement<System, movedata> mov(deg2rad(60.0), deg2rad(75.0), 200, diff/2);

	while (window.isOpen())
	{
		currentTime = clock.getElapsedTime();
		sf::Event event;
		

		while (window.pollEvent(event))
		{
			if (event.type == sf::Event::Closed)
				window.close();

			if (event.type == sf::Event::Resized)
			{
				sf::FloatRect visibleArea(0, 0, event.size.width, event.size.height);
				window.setView(sf::View(visibleArea));
			}
		}

		sf::CircleShape shape(30.f);
		shape.setFillColor(sf::Color::White);
		sf::CircleShape shape2(30.f);
		shape2.setFillColor(sf::Color::White);

		shape.setPosition(300.0 * std::sin(mov.getcurrentstep().theta) - 30.0 + 600.0, 300 * std::cos(mov.getcurrentstep().theta) - 30.0);
		shape2.setPosition(300.0 * std::sin(mov.getcurrentstep().theta) + 600.0 + 300.0 * std::sin(mov.getcurrentstep().pi) - 60.0, 300.0 * std::cos(mov.getcurrentstep().theta) + 300.0 * std::cos(mov.getcurrentstep().pi) - 60.0);

		sf::Vertex line[] =
		{
			sf::Vertex(sf::Vector2f(600.0f, 0.0f)),
			sf::Vertex(sf::Vector2f(300.0 * std::sin(mov.getcurrentstep().theta) + 600.0, 300.0 * std::cos(mov.getcurrentstep().theta)))
		};

		sf::Vertex line2[] =
		{
			sf::Vertex(sf::Vector2f(300.0 * std::sin(mov.getcurrentstep().theta) + 600.0, 300 * std::cos(mov.getcurrentstep().theta))),
			sf::Vertex(sf::Vector2f(300.0 * std::sin(mov.getcurrentstep().theta) + 600.0 + 300 * std::sin(mov.getcurrentstep().pi) - 30.0, 300.0 * std::cos(mov.getcurrentstep().theta) + 300.0 * std::cos(mov.getcurrentstep().pi) - 30.0))
		};

		mov.deletesteps(1);
		if (mov.stepnum() < 200)
		{
			mov.calcsteps(5, diff/2);
		}

		window.clear();
		
		window.draw(shape);
		window.draw(line, 2, sf::Lines);

		window.draw(shape2);
		window.draw(line2, 2, sf::Lines);

		window.display();
		
		fps = 1.0f / (currentTime.asSeconds() - previousTime.asSeconds());
		previousTime = currentTime;
	}
	return 0;
}

 

빨리 만드느라 주석은 안넣었지만 아래에 대충 설명하겠다.

 

사전 지식: 변분법, 라그랑주 역학

ellipsoid.tistory.com/entry/%EB%B3%80%EB%B6%84%EB%B2%95

변분법

 

이중 진자의 라그랑지안은

대충 이렇다

이걸 극좌표로 바꾸어 주면

이렇게 된다(a는 각도)

이제 이 라그랑지안을

요렇게 생긴 오일러-라그랑주 방정식에 넣으면 되는데 넣어보면

요런 식 2개가 나온다

근데 이런 형태는 컴퓨터로 풀기 적당하지 않다. 변수가 2개에 a1의 미분같은것도 있으니까

그래서 르장드르 변환해서 해밀토니안 형태로 써야한다

 

 

요렇게, 여기서 p는 일반화 운동량이다.

A1, A2는 상수로

대충 이런 값을 갖는다. (편의상 l1=l2=l이라고 하자)

사실 나도 유도하기 귀찮아서 www.math24.net/double-pendulum/ 여기서 가져왔다 어짜피 시간 좀 들이면 다 할 수 있긴 하니까

 

이 네개의 일계 미분 방정식을 룽게-쿠타법으로 잘 풀면 된다.

변분법 관련해서 블로그에 적도록 하겠다.

 

코드는 나중에 깃허브에 추가하겠다.