Simulating ASCII Double Pendulums in C++
Finishing off the theme of “chaos” that persists within my earlier articles, I’d like to explore how to simulate the motion of double pendulums using the console and C++. The implementation closely follows this resource, so I highly recommend taking a look at that. To make things simpler, I won’t be using a time library to calculate each frame of motion. You will also need a Windows machine to do this since it relies on Windows.h! to draw to the console.
What is a double pendulum?
We think of a double pendulum as a simple way to demonstrate chaotic motion. You can imagine it as a set of two strings of an arbitrary length, each with a ball of a varying mass attached to the end. Given some set of initial conditions, we can use an equation that takes those conditions and predicts the next frame of motion.
Implementation
Global Variables
In order to write characters to the console, we will need a few global variables alongside some other constants that will be used later in the program. (that means, make sure to include Windows.h!)
// These two variables control the size of the window.
int xSize = 100;
int ySize = 100;// These two variables control the position of the pendulum within the window.
int offsetx = 25;
int offsety = -32;// This is the gravity constant.
float g = 1;// These two variables control the time-step and dampening, respectively.
float dt = 0.3;
float DAMPENING = 0.997;// These variables handle writing characters directly to the console.
HANDLE hOutput = GetStdHandle(STD_OUTPUT_HANDLE);
DWORD dwWritten;
HWND myconsole = GetConsoleWindow();
HDC dc = GetDC(myconsole);
Next, let’s define a struct called a Pendulum. Inside it, we will add some variables that we would want to modify from the introduction above. These include length, angle, position, velocity, and acceleration values for each pendulum. We also store the mass of the circles at the end of each pendulum. For the purposes of this implementation, “p1” stands for the first pendulum, and “p2” stands for the second pendulum.
struct Pendulum
{
float p1Length = 10;
float p2Length = 15; float p1angle = M_PI/2;
float p2angle = M_PI/2; float p1x = p1Length*sin(p1angle);
float p1y = p1Length*cos(p1angle); float p2x = p2Length*sin(p2angle);
float p2y = p2Length*cos(p2angle); float p1v = 0;
float p2v = 0; float p1a = 0;
float p2a = 0; float m1 = 2;
float m2 = 2;
};
Update Function
Now we will need an update function that runs within a while(true) loop, which will calculate the values of the position, velocity, and acceleration quantities for both pendulums. It also controls how the pendulums are displayed within the console. The basic idea behind the display is as follows:
Using Windows’ FillConsoleOutputCharacter function with Bresenham’s line drawing algorithm, we can accomplish this display pretty easily by redrawing each frame with the updated values outputted from the formulae which I’ll talk about later. However, there’s a small problem. Once a frame of motion is drawn on the screen, it stays there forever because there’s nothing to replace it. We need some way to handle the cleanup of the last frame, so it doesn’t look like this:
To do this, we can store the X and Y coordinates of each iteration with std::pair<int,int>, and store each pair within a std::vector. Doing this for both pendulums allows us to iterate through both vectors and remove each character by replacing it with an empty space at the end of each iteration. The rest of the function simply updates each of the values using the formula shown below.
This equation of motion calculates the acceleration of each pendulum, which is added to the velocity. The velocity changes the angle of the pendulum, and all these values are displayed once again using the method discussed above. There’s another problem, though. The while(true) loop iterates really quickly, so the acceleration can easily get out of control and cause the program to crash due to a NaN error. We need some way to scale down the equations, and so we multiply each output by some value that is called dt (we can think of this as the time-step). I’ve found that 0.3 works best.
I’ve also multiplied the velocity by a dampening factor so that it gradually comes to a stop. I’ve found that 0.997 works quite well. The program is still a bit too fast when displaying each iteration, so a quick ad-hoc solution to fix this is to add a Sleep command. I won’t go over how to implement Bresenham, but my recommendation is to add it as a method within the Pendulum struct to use it within the update function like I ended up doing. I think that covers everything within the update function, so here’s the code:
void update(Pendulum *pendulum)
{while (true)
{
vector<pair<int,int>> v;
vector<pair<int,int>> v2; pendulum->p1x = pendulum->p1Length*sin(pendulum->p1angle);
pendulum->p1y = pendulum->p1Length*cos(pendulum->p1angle); pendulum->p2x = pendulum->p2Length*sin(pendulum->p2angle);
pendulum->p2y = pendulum->p2Length*cos(pendulum->p2angle); v = pendulum->bresenham(0,pendulum->p1x,0,pendulum- >p1y);
v2 = pendulum->bresenham2(pendulum->p1x,pendulum->p2x+pendulum->p1x,pendulum->p1y,pendulum->p2y+pendulum->p1y); FillConsoleOutputCharacter(hOutput, char(219), 1, {offsetx+xSize/2,(ySize/2)+offsety}, &dwWritten); pendulum->p1v += pendulum->p1a;
pendulum->p2v += pendulum->p2a; pendulum->p1angle += pendulum->p1v;
pendulum->p2angle += pendulum->p2v; float num1 = -g*(2*pendulum->m1+pendulum->m2)*sin(pendulum- >p1angle)-pendulum->m2*g*sin(pendulum->p1angle-2*pendulum->p2angle)
-2*sin(pendulum->p1angle-pendulum- >p2angle)*pendulum->m2*(pendulum->p2v*pendulum->p2v*pendulum->p2Length+pendulum->p1v*pendulum->p1v*pendulum->p1Length*cos(pendulum->p1angle-pendulum->p2angle)); float den1 = pendulum->p1Length*(2*pendulum->m1+pendulum->m2-pendulum->m2*cos(2*pendulum->p1angle-2*pendulum->p2angle)); float num2 = 2*sin(pendulum->p1angle-pendulum->p2angle)*(pendulum->p1v*pendulum->p1v*pendulum->p1Length*(pendulum->m1+pendulum->m2)+g*(pendulum->m1+pendulum->m2)*cos(pendulum->p1angle)+pendulum->p2v*pendulum->p2v*pendulum->p2Length*pendulum->m2*cos(pendulum->p1angle-pendulum->p2angle)); float den2 = pendulum->p2Length*(2*pendulum->m1+pendulum->m2-pendulum->m2*cos(2*pendulum->p1angle-2*pendulum->p2angle)); pendulum->p1a = (num1/den1)*dt;
pendulum->p2a = (num2/den2)*dt; pendulum->p1v = pendulum->p1v*DAMPENING;
pendulum->p2v = pendulum->p2v*DAMPENING; Sleep(25); for (int i = 0; i < v.size(); i++)
{
FillConsoleOutputCharacter(hOutput, ' ', 1, {(xSize/2)+v[i].first,(ySize/2)+v[i].second}, &dwWritten);
} for (int i = 0; i < v2.size(); i++)
{
FillConsoleOutputCharacter(hOutput, ' ', 1, {(xSize/2)+v2[i].first,(ySize/2)+v2[i].second}, &dwWritten);
}
}
}
Now the program can display one pendulum. While this is quite cool on its own, I was searching to see who else had done an ASCII implementation and found this video by Dino1729, in which he has the functionality to add multiple pendulums. I decided to figure out how I could implement that myself.
Adding More Pendulums
One way to do this is to use threads. Using kbhit() and a while(true) loop, we can continually check if a particular key has been pressed and use that to create a new Pendulum object and spawn a new thread to handle it. We will then add the thread into the vector container and call the update function, passing in the new object. In order to synchronize the threads, we will use another std::vector and an enhanced for loop to join each thread after the while loop. We’ll also be doing this in main().
int main()
{ vector<thread> tv; while(true)
{
if(kbhit())
{
char ch = getch(); if (ch == 'g') // you can use whatever key you want here.
{
Pendulum *temp = new Pendulum();
tv.emplace_back([&](){update(temp);});
}
}
}
for (auto& it : tv)
{
it.join();
}
}
An example of the code running with multiple pendulums can be found at the start of the article. Another run is shown below. We can see how even if new pendulums are spawned very close to each other, over time their motion becomes increasingly different.
Conclusions
One of the best ways to solidify a mathematics or physics concept in your head is to find a way to implement it in code. With that being said, I hope that this post inspires people to implement this and helps people understand the concept of how to implement it a little bit better. As for me, it might be time to move on from coding console simulations in C++, and switch to a more practical library and language for this sort of thing like Processing using Java or JavaScript.
