最速降下線の問題を解きサイクロイドを描く
文責：八重樫 和之 （2011年3月18日）
カテゴリ：ゴールドスタイン(3)
最速降下線の問題
最速降下線の問題は高さの違う二点間を，最初に高い点で静止していた質点が， 重力の作用の下で，低い点の方に向かって移動する場合に，最も短い時間で到達することのできるような曲線を求めるという問題である． 質点が曲線に沿って動く速さをvとすると，長さdsに落ちて行くのに要する時間はds/vである． そうすると問題は積分
の極小値を求めることになる．質点のエネルギー保存則から
と書くことができる．そうすると
よりt_{12}に対する式は
となり,fは
であることがわかる． これを，ラグランジュの方程式を用いて解くと，θを用いて
と表される．ただしAは定数である．これを解くと質点を離した点を先点とするサイクロイドが得られる．
実行結果
上記式を解いた結果を下記に示す．
プログラムソース
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 | #include <math.h>#include <fstream>#include <sstream>#include <iostream>#include <direct.h>#include <time.h>#include <glut.h>#include <gl_screenshot.h> using namespace std;int WindowPositionX = 200; //生成するウィンドウ位置のX座標int WindowPositionY = 200; //生成するウィンドウ位置のY座標int WindowWidth = 512; //生成するウィンドウの幅int WindowHeight = 512; //生成するウィンドウの高さchar WindowTitle[] = "最速降下線の問題"; //ウィンドウのタイトルstatic GLfloat floor_planar[4];static GLfloat floor_s = 50.0f;static GLfloat pM[16];static GLfloat lightpos[4] = { -30, -100, 50, 1 };typedef struct _QUADS_VERTEX{ GLfloat v0[3]; GLfloat v1[3]; GLfloat v2[3]; GLfloat v3[3];}QUADS_VERTEX;static QUADS_VERTEX floor_v = { { floor_s, floor_s, -1.0f }, { -floor_s, floor_s, -1.0f }, { -floor_s, -floor_s, -1.0f }, { floor_s, -floor_s, -1.0f },};#define _BITMAP 1int tn = 0;double dt = 0.05;gl_screenshot gs; //bmpファイルの出力struct { double x, y, z; double vx, vy, vz; double theta;}p[100];int pn = 0;double ax = 0.0 , ay = 0.0 , az = -4.0;double vx = 5.0 , vy = 5.0 , vz = 20.0; double hanpatu = 0.9;//----------------------------------------------------// 物質質感の定義//----------------------------------------------------struct MaterialStruct { GLfloat ambient[4]; GLfloat diffuse[4]; GLfloat specular[4]; GLfloat shininess;};//jade(翡翠)MaterialStruct ms_jade = { {0.135, 0.2225, 0.1575, 1.0}, {0.54, 0.89, 0.63, 1.0}, {0.316228, 0.316228, 0.316228, 1.0}, 12.8};//ruby(ルビー)MaterialStruct ms_ruby = { {0.1745, 0.01175, 0.01175, 1.0}, {0.61424, 0.04136, 0.04136, 1.0}, {0.727811, 0.626959, 0.626959, 1.0}, 76.8};int list;//----------------------------------------------------// 関数プロトタイプ（後に呼び出す関数名と引数の宣言）//----------------------------------------------------void Initialize(void);void Display(void);void Idle();void Keyboard(unsigned char key, int x, int y);void Ground(void); //大地の描画void Axis(void);//軸の描画void findPlane(GLfloat plane[4], GLfloat v0[3], GLfloat v1[3], GLfloat v2[3]); void shadowMatrix(GLfloat *m, GLfloat plane[4], GLfloat light[4]);void DrawFloor(bool bTexture);void DrawShadow(void);void DrawStructure(bool);void DRAW_STRING(int x, int y, char *string, void *font = GLUT_BITMAP_TIMES_ROMAN_24);void DISPLAY_TEXT(int x, int y, char *string);//----------------------------------------------------// メイン関数//----------------------------------------------------int main(int argc, char *argv[]){ srand((unsigned)time(NULL));#if _BITMAP _mkdir("bitmap"); //bmpファイル保存用のフォルダの作成#endif glutInit(&argc, argv);//環境の初期化 glutInitWindowPosition(WindowPositionX, WindowPositionY);//ウィンドウの位置の指定 glutInitWindowSize(WindowWidth, WindowHeight); //ウィンドウサイズの指定 glutInitDisplayMode(GLUT_RGBA | GLUT_DEPTH | GLUT_DOUBLE);//ディスプレイモードの指定 glutCreateWindow(WindowTitle); //ウィンドウの作成 glutDisplayFunc(Display); //描画時に呼び出される関数を指定する（関数名：Display） glutKeyboardFunc(Keyboard);//キーボード入力時に呼び出される関数を指定する（関数名：Keyboard） glutIdleFunc(Idle); //プログラムアイドル状態時に呼び出される関数 Initialize(); //初期設定の関数を呼び出す glutMainLoop(); return 0;}//----------------------------------------------------// 初期設定の関数//----------------------------------------------------void Initialize(void){ glClearColor(1.0, 1.0, 1.0, 1.0); //背景色 glEnable( GL_DEPTH_TEST ); //デプスバッファを使用：glutInitDisplayMode() で GLUT_DEPTH を指定する glDepthFunc( GL_LEQUAL ); glClearDepth( 1.0 ); findPlane( floor_planar, floor_v.v0, floor_v.v1, floor_v.v2 ); //透視変換行列の設定------------------------------ glMatrixMode(GL_PROJECTION);//行列モードの設定（GL_PROJECTION : 透視変換行列の設定、GL_MODELVIEW：モデルビュー変換行列） glLoadIdentity();//行列の初期化 gluPerspective(30.0, (double)WindowWidth/(double)WindowHeight, 0.1, 1000.0); //透視投影法の視体積gluPerspactive(th, w/h, near, far);}//----------------------------------------------------// 描画の関数//----------------------------------------------------void Display(void) { glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT); //視点の設定------------------------------ gluLookAt( 0.0, 0.0, 150.0, // 視点の位置x,y,z; 0.0, 10.0, 0.0, // 視界の中心位置の参照点座標x,y,z 0.0, 0.0, 1.0 ) ; //視界の上方向のベクトルx,y,z //---------------------------------------- //モデルビュー変換行列の設定-------------------------- glMatrixMode(GL_MODELVIEW);//行列モードの設定（GL_PROJECTION : 透視変換行列の設定、GL_MODELVIEW：モデルビュー変換行列） glLoadIdentity();//行列の初期化 glViewport(0, 0, WindowWidth, WindowHeight); //---------------------------------------------- //ステンシルバッファクリア値の設定-------------------------- glClearStencil( 0 ); glCullFace( GL_BACK ); glEnable( GL_CULL_FACE ); glEnable( GL_AUTO_NORMAL ); glEnable( GL_NORMALIZE ); //---------------------------------------- // 平面射影行列の算出-------------------------- shadowMatrix(pM,floor_planar,lightpos); //-------------------------- // 光源ON----------------------------- glEnable( GL_LIGHTING ); glEnable( GL_LIGHT0 ); glLightfv( GL_LIGHT0,GL_POSITION,lightpos ); //----------------------------------- glPushMatrix(); DrawStructure(false); DrawShadow(); glPopMatrix(); glDisable(GL_AUTO_NORMAL); glDisable(GL_NORMALIZE); //陰影OFF----------------------------- glDisable(GL_LIGHTING); //----------------------------------- Ground(); Axis(); //文字の描画 char t_char[20]; char t_char2[20]; strcpy_s(t_char2, "t = "); sprintf_s(t_char, "%d", tn); strcat_s(t_char2, t_char); DISPLAY_TEXT(5, 95, t_char2 );#if _BITMAP ostringstream fname; int tt = tn +10000; fname << "bitmap/" << tt << ".bmp" ;//出力ファイル名 string name = fname.str(); gs.screenshot(name.c_str(), 24);#endif glutSwapBuffers(); //glutInitDisplayMode(GLUT_DOUBLE)でダブルバッファリングを利用可 tn++;}//----------------------------------------------------// 物体の描画//----------------------------------------------------void DrawStructure(bool flag){ int i; p[1].theta += dt ; p[1].x = 2 * (p[1].theta - sin(p[1].theta)); p[1].y = - 2 * (1 - cos(p[1].theta)); for(int i=1; i<=pn; i++){ if(!flag || p[i].z >0){ glPushMatrix(); glMaterialfv(GL_FRONT, GL_AMBIENT, ms_ruby.ambient); glMaterialfv(GL_FRONT, GL_DIFFUSE, ms_ruby.diffuse); glMaterialfv(GL_FRONT, GL_SPECULAR, ms_ruby.specular); glMaterialfv(GL_FRONT, GL_SHININESS, &ms_ruby.shininess); glTranslated(p[i].x , p[i].y , p[i].z );//平行移動値の設定 glutSolidSphere(1.0, 20, 20);//引数：(半径, Z軸まわりの分割数, Z軸に沿った分割数) glPopMatrix(); } }}//----------------------------------------------------// 大地の描画//----------------------------------------------------void Ground(void) { double ground_max_x = 300.0; double ground_max_y = 300.0; glColor3d(0.8, 0.8, 0.8); // 大地の色 glBegin(GL_LINES); for(double ly = -ground_max_y ;ly <= ground_max_y; ly+=10.0){ glVertex3d(-ground_max_x, ly, -1.1); glVertex3d(ground_max_x, ly , -1.1); } for(double lx = -ground_max_x ;lx <= ground_max_x; lx+=10.0){ glVertex3d(lx, ground_max_y , -1.1); glVertex3d(lx, -ground_max_y, -1.1); } glEnd();}//----------------------------------------------------//軸の描画//----------------------------------------------------void Axis(void) { glColor3d(.0, .0, .0); glBegin(GL_LINES); glVertex3d(-100, 0, -1.1); glVertex3d(100, 0 , -1.1); glEnd(); glBegin(GL_LINES); glVertex3d(0, -100, -1.1); glVertex3d(0, 100 , -1.1); glEnd(); glBegin(GL_POINTS); for(int i=0; i < 1000 ; i++){ glVertex3d(2 * (i * 0.1 - sin(i * 0.1)), - 2 * (1 - cos(i * 0.1)), -1.1); } glEnd();}//----------------------------------------------------// アイドル時に呼び出される関数//----------------------------------------------------void Idle(){ glutPostRedisplay(); //glutDisplayFunc()を１回実行する}//----------------------------------------------------// キーボード入力時に呼び出される関数//----------------------------------------------------void Keyboard(unsigned char key, int x, int y){ switch ( key ) { case 'a': pn++; p[1].theta = 0; p[pn].x = -0.0; p[pn].y = -10.0; p[pn].z = 5.0; p[pn].vx = vx * ( (double)rand()/(double)RAND_MAX - (double)rand()/(double)RAND_MAX ); p[pn].vy = vy * ( (double)rand()/(double)RAND_MAX - (double)rand()/(double)RAND_MAX ); p[pn].vz = vz * ( (double)rand()/(double)RAND_MAX ); break; case 'q': exit(0); break; default: break; }}//----------------------------------------------------// 床平面の方程式と行列の計算//----------------------------------------------------void findPlane( GLfloat plane[4], // 作成する平面方程式の係数 GLfloat v0[3], // 頂点１ GLfloat v1[3], // 頂点２ GLfloat v2[3]) // 頂点３{ GLfloat vec0[3], vec1[3]; // Need 2 vectors to find cross product. vec0[0] = v1[0] - v0[0]; vec0[1] = v1[1] - v0[1]; vec0[2] = v1[2] - v0[2]; vec1[0] = v2[0] - v0[0]; vec1[1] = v2[1] - v0[1]; vec1[2] = v2[2] - v0[2]; // find cross product to get A, B, and C of plane equation plane[0] = vec0[1] * vec1[2] - vec0[2] * vec1[1]; plane[1] = -(vec0[0] * vec1[2] - vec0[2] * vec1[0]); plane[2] = vec0[0] * vec1[1] - vec0[1] * vec1[0]; plane[3] = -(plane[0] * v0[0] + plane[1] * v0[1] + plane[2] * v0[2]);}void shadowMatrix( GLfloat *m, // 作成する行列のポインタ GLfloat plane[4], // 射影する表面の平面方程式の係数 GLfloat light[4]) // 光源の同時座標値{ GLfloat dot; // Find dot product between light position vector and ground plane normal. dot = plane[0] * light[0] + plane[1] * light[1] + plane[2] * light[2] + plane[3] * light[3]; m[0] = dot - light[0] * plane[0]; m[4] = 0.f - light[0] * plane[1]; m[8] = 0.f - light[0] * plane[2]; m[12] = 0.f - light[0] * plane[3]; m[1] = 0.f - light[1] * plane[0]; m[5] = dot - light[1] * plane[1]; m[9] = 0.f - light[1] * plane[2]; m[13] = 0.f - light[1] * plane[3]; m[2] = 0.f - light[2] * plane[0]; m[6] = 0.f - light[2] * plane[1]; m[10] = dot - light[2] * plane[2]; m[14] = 0.f - light[2] * plane[3]; m[3] = 0.f - light[3] * plane[0]; m[7] = 0.f - light[3] * plane[1]; m[11] = 0.f - light[3] * plane[2]; m[15] = dot - light[3] * plane[3];}//----------------------------------------------------// 床の描画と影の描画//----------------------------------------------------void DrawFloor(bool bTexture){ if( bTexture ){ // 床にテクスチャを使う時はココで設定する // glBindTexture( GL_TEXTURE_2D, ); glDisable(GL_LIGHTING); glBegin(GL_QUADS); // glTexCoord2f( , ); glVertex3fv( floor_v.v0 ); // glTexCoord2f( , ); glVertex3fv( floor_v.v1 ); // glTexCoord2f( , ); glVertex3fv( floor_v.v2 ); // glTexCoord2f( , ); glVertex3fv( floor_v.v3 ); glEnd(); glEnable(GL_LIGHTING); }else{ glDisable(GL_LIGHTING); glBegin(GL_QUADS); glVertex3fv( floor_v.v0 ); glVertex3fv( floor_v.v1 ); glVertex3fv( floor_v.v2 ); glVertex3fv( floor_v.v3 ); glEnd(); glEnable(GL_LIGHTING); }}void DrawShadow(void){ ///////////////////////////////////////////// //床のステンシルを付ける glEnable(GL_STENCIL_TEST); glStencilFunc( GL_ALWAYS, 1, ~0); //これから描画するもののステンシル値にすべて１タグをつける glStencilOp(GL_KEEP,GL_KEEP ,GL_REPLACE); glColor4f(0.7f, 0.4f, 0.0f, 1.0f); //DrawFloor( false );//床の描画 ///////////////////////////////////////////// //カラー・デプスバッファマスクをセットする //これで以下の内容のピクセルの色の値は、書き込まれない。 glColorMask(0,0,0,0); glDepthMask(0); ///////////////////////////////////////////// //床にオブジェクトの影のステンシルを付ける glEnable(GL_STENCIL_TEST); glStencilFunc( GL_EQUAL, 1, ~0); //これから描画するもののステンシル値にすべて１タグをつける glStencilOp(GL_KEEP,GL_KEEP ,GL_INCR); glDisable(GL_DEPTH_TEST); glPushMatrix(); glMultMatrixf(pM); DrawStructure(true); glPopMatrix(); glEnable(GL_DEPTH_TEST); ///////////////////////////////////////////// //ビットマスクを解除 glColorMask(1,1,1,1); glDepthMask(1); ///////////////////////////////////////////// //影をつける glStencilFunc( GL_EQUAL, 2, ~0 ); glEnable(GL_BLEND); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glColor4f(0.1f, 0.1f, 0.1f, 0.5f); glDisable(GL_DEPTH_TEST); DrawFloor( false );//床の描画 glEnable(GL_DEPTH_TEST); glDisable(GL_BLEND); glDisable(GL_STENCIL_TEST);}//////////////////////////////////////////////////////////////////////////// 文字描画//////////////////////////////////////////////////////////////////////////void DISPLAY_TEXT(int x, int y, char *string){glDisable(GL_LIGHTING);glDisable(GL_LIGHT0); glPushAttrib(GL_ENABLE_BIT); glMatrixMode(GL_PROJECTION); glPushMatrix(); glLoadIdentity(); gluOrtho2D(0, 100, 0, 100); glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadIdentity(); glColor3f(0.0, 0.0, 0.0); glCallList(list); glPopMatrix(); glMatrixMode(GL_PROJECTION); glPopMatrix(); glPopAttrib(); glMatrixMode(GL_MODELVIEW); list=glGenLists(1); glNewList(list,GL_COMPILE); DRAW_STRING(x, y, string , GLUT_BITMAP_TIMES_ROMAN_24); glEndList();glEnable(GL_LIGHTING);glEnable(GL_LIGHT0);}void DRAW_STRING(int x, int y, char *string, void *font){ int len, i; glRasterPos2f(x, y); len = (int) strlen(string); for (i = 0; i < len; i++){ glutBitmapCharacter(font, string[i]); }}</gl_screenshot.h></glut.h></time.h></direct.h></iostream></sstream></fstream></math.h> |
参考
- ゴールドスタイン (著), 瀬川 富士 (翻訳) 『古典力学 (上) (物理学叢書 (11a)) 』吉岡書店,1983,p54-56
- 物理のかぎしっぽ-最速降下曲線
- VisualC++ を使った OpenGL 入門【９日目】 文字
