Simulácia reálneho kyvadla

.exe súbor (1.2MB) odkaz

Vizualizácia pohybu reálneho kyvadla. Na pozadí programu beží nelineárna diferenciálna rovnica, ktorá ho vykresľuje v OpenGL. Rovnica nemá štandardné riešenie, čiže jediná cesta je ju po drobných časových krokoch počítať a výsledok vložiť späť do rovnice a postup opakovať. (Jedná sa o spustiteľný .exe program)

Tu je kompletný zdrojový kód (keby si niekto chcel pozrieť, ako funguje narábanie s neriešiteľnými nelineárnymi diferenciálnymi rovnicami v numerickej simulácii):

#define FREEGLUT_STATIC
#include <GL/freeglut.h>
#include <math.h>
#include <stdio.h>

// Definícia konštánt
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// Fyzikálne parametre kyvadla
const float g = 9.81f; // Gravitačné zrýchlenie (m/s^2)
const float L = 2.0f; // Dĺžka závesu (m)

// Stavové veličiny kyvadla
float theta = 3.1f; // Počiatočný uhol v radiánoch (cca 177 stupňov - extrémna nelinearita)
float omega = 0.0f; // Počiatočná uhlová rýchlosť (rad/s)

// Časový krok simulácie
const float dt = 0.01f; // 10 milisekúnd na krok

// Funkcia na inicializáciu OpenGL
void init(void) {
glClearColor(0.1f, 0.1f, 0.1f, 1.0f); // Tmavosivé pozadie
glShadeModel(GL_FLAT);
}

// Funkcia, ktorá prepočítava fyziku (beží na pozadí)
void updatePhysics(void) {
// Presná nelineárna rovnica: alpha = -(g / L) * sin(theta)
float alpha = -(g / L) * sinf(theta);

// Numerická integrácia (Metóda Euler-Cromer)
omega += alpha * dt;
theta += omega * dt;

// Normalizácia uhla do intervalu <-PI, PI>
if (theta > M_PI) theta -= 2.0f * M_PI;
if (theta < -M_PI) theta += 2.0f * M_PI;

// Vynútenie prekreslenia obrazovky
glutPostRedisplay();
}

// Funkcia na vykreslenie scény
void display(void) {
glClear(GL_COLOR_BUFFER_BIT);
glLoadIdentity();

// Prepočet polohy konca kyvadla z polárnych súradníc na karteziánske
// V OpenGL smeruje os -Y dole, preto: x = L * sin(theta), y = -L * cos(theta)
float x = L * sinf(theta);
float y = -L * cosf(theta);

// 1. Vykreslenie závesu (čiara)
glColor3f(0.7f, 0.7f, 0.7f); // Sivá farba chvosta
glLineWidth(3.0f);
glBegin(GL_LINES);
glVertex2f(0.0f, 0.0f); // Bod uchytenia (stred)
glVertex2f(x, y); // Pozícia závažia
glEnd();

// 2. Vykreslenie uchytenia (malý štvorec v strede)
glColor3f(1.0f, 1.0f, 1.0f);
glBegin(GL_QUADS);
glVertex2f(-0.05f, -0.05f);
glVertex2f( 0.05f, -0.05f);
glVertex2f( 0.05f, 0.05f);
glVertex2f(-0.05f, 0.05f);
glEnd();

// 3. Vykreslenie závažia (kruh aproximovaný n-uholníkom)
glColor3f(0.2f, 0.6f, 1.0f); // Svetlomodrá farba
float radius = 0.15f;
glBegin(GL_POLYGON);
for (int i = 0; i < 32; i++) {
float angle = 2.0f * M_PI * i / 32.0f;
glVertex2f(x + radius * cosf(angle), y + radius * sinf(angle));
}
glEnd();

glutSwapBuffers();
}

// Funkcia na ošetrenie zmeny veľkosti okna
void reshape(int w, int h) {
glViewport(0, 0, (GLsizei)w, (GLsizei)h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();

// Nastavenie ortografického premietania (udržujeme pomer strán okna)
float aspectRatio = (float)w / (float)h;
if (w >= h) {
// Okno je široké
gluOrtho2D(-2.5f * aspectRatio, 2.5f * aspectRatio, -2.5f, 2.5f);
} else {
// Okno je vysoké
gluOrtho2D(-2.5f, 2.5f, -2.5f / aspectRatio, 2.5f / aspectRatio);
}
glMatrixMode(GL_MODELVIEW);
}

void timer(int value) {
// Spustíme prepočet fyziky
updatePhysics();

// Znovu zaregistrujeme tento časovač, aby sa spustil o 10 milisekýnd (10 ms = 100 FPS)
glutTimerFunc(10, timer, 0);
}

// Hlavná funkcia
int main(int argc, char** argv) {
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
glutInitWindowSize(600, 600);
glutInitWindowPosition(100, 100);
glutCreateWindow("Simulacia nelinearneho kyvadla (OpenGL 1.1)");

init();

// Registrácia callback funkcií
glutDisplayFunc(display);
glutReshapeFunc(reshape);

glutTimerFunc(10, timer, 0); //Hlavná počítacia funkcia timer v 10ms intervaloch

glutMainLoop();
return 0;
}

Pri výukách kyvadla sa používa tzv. "matematické kyvadlo", kde sa odstraňuje nelineárny člen dif. rovnice obsahujúci sin(theta). Berie do úvahy len drobný kúsok trajektórie kyvadla. Využíva sa pri tom skutočnosť, že sínus malého uhla (v radiánoch), sa rovná približne hodnote toho uhla. Môžete si to vyskúšať na kalkulačke. Napr. sin 0.1 = 0.0998... Tým odpadne z rovnice sínus a taká rovnica už má riešenie. Lenže taká rovnica nedokáže popísať kyvadlo, ktoré pustím z veľkého uhla, čiže od reality má ďaleko.

Hlavný výpočet prebieha v týchto riadkoch:

// Fyzikálne parametre kyvadla
const float g = 9.81f; // Gravitačné zrýchlenie (m/s^2)
const float L = 2.0f; // Dĺžka závesu (m)

// Časový krok simulácie
const float dt = 0.01f; // 10 milisekúnd na krok

float alpha = -(g / L) * sinf(theta); // 1. Spočítaj zrýchlenie v danom momente
omega += alpha * dt; // 2. Aktualizuj rýchlosť (prirátaj zrýchlenie za čas dt)
theta += omega * dt; // 3. Aktualizuj polohu (posuň uhol o rýchlosť za čas dt)

Samozrejme - v realite ešte dochádza k stratám energie (napr. trenie o vzduch) a kmity spomaľujú a zmenšujú sa. Ale to som už nechcel zozložiťovať kód...