/* * heat1d-e-glut.c -- 1次元熱伝導方程式の初期値境界値問題を陽解法で解く。 * gcc -I/usr/local/include -o heat1d heat1d-e-glut.c -lglut -lGL -lX11 -lm * * オリジナルは fplot ライブラリィを利用した * http://nalab.mind.meiji.ac.jp/%7Emk/program/fdm/heat1d-e.c */ #define WIDTH (600) #define HEIGHT (600) #include <stdio.h> #include <stdlib.h> #include <math.h> #include <GL/glut.h> double xc = 0.0, yc = 0.0, ratio_x = 0.9, ratio_y = 0.9; void screen(double xl, double yb, double xr, double yt) { /* [xl,xr]×[yb,yt] を [-0.9,0.9]×[-0.9,0.9] */ xc = (xl + xr) / 2; yc = (yb + yt) / 2; ratio_x = 1.8 / (xr - xl); ratio_y = 1.8 / (yt - yb); } double X(double x) { return ratio_x*(x-xc); } double Y(double y) { return ratio_y*(y-yc); } void keyboard(unsigned char key, int x, int y) { switch (key) { case 'q': case 'Q': case '\033': exit(0); default: break; } } void init(void) { glClearColor(0.0, 0.0, 1.0, 1.0); } void resize(int w, int h) { /* ウィンドウ全体をヴューポートにする */ glViewport(0, 0, w, h); /* 変換行列の初期化 */ glLoadIdentity(); /* */ glOrtho(- w / WIDTH, w / WIDTH, -h / HEIGHT, h / HEIGHT, -1.0, 1.0); } void drawBitmapString(double x, double y, void *font, char *str) { glRasterPos2d(X(x), Y(y)); while (*str) glutBitmapCharacter(font, *str++); } int N; double lambda, Tmax; void display(void) { int i, n, nMax; double tau, h; double *u, *newu; double f(double); int win; char message[100]; /* h, τ を計算する */ h = 1.0 / N; tau = lambda * h * h; printf("τ=%g\n", tau); /* ベクトル u, newu を用意する */ u = malloc(sizeof(double) * (N+1)); newu = malloc(sizeof(double) * (N+1)); /* 初期値の代入 */ for (i = 0; i <= N; i++) u[i] = f(i * h); screen(-0.2, -0.2, 1.2, 1.2); glClear(GL_COLOR_BUFFER_BIT); /* 座標軸を表示する */ glColor3d(0.0, 1.0, 0.0); glBegin(GL_LINES); glVertex2d(X(-0.1), Y(0.0)); glVertex2d(X( 1.1), Y(0.0)); glVertex2d(X(0.0), Y(-0.1)); glVertex2d(X(0.0), Y( 1.1)); glEnd(); /* t=0 の状態を表示する */ glColor3d(1.0, 0.0, 0.0); glBegin(GL_LINE_LOOP); glVertex2d(X(0.0), Y(u[0])); for (i = 1; i <= N; i++) glVertex2d(X(i * h), Y(u[i])); glEnd(); /* タイトルと入力パラメーターを表示する */ #ifdef NONE eggx_drawstr(win, 0.1, 0.8, 14, 0, "heat equation, homogeneous Dirichlet boundary condition"); sprintf(message, "N=%d, lambda=%g, Tmax=%g", N, lambda, Tmax); eggx_drawstr(win, 0.1, 0.7, 14, 0, message); #endif drawBitmapString(0.1, 0.8, GLUT_BITMAP_TIMES_ROMAN_24, "heat equation, u(0,t)=u(1,t)=0"); /* ループを何回まわるか計算する (四捨五入) */ nMax = rint(Tmax / tau); /* 時間に関するステップを進めるループ */ for (n = 0; n < nMax; n++) { /* 差分方程式 (n -> n+1) */ for (i = 1; i < N; i++) newu[i] = (1.0 - 2 * lambda) * u[i] + lambda * (u[i+1] + u[i-1]); /* 計算値を更新 */ for (i = 1; i < N; i++) u[i] = newu[i]; /* Dirichlet 境界条件 */ u[0] = u[N] = 0.0; /* この時刻 (t=(n+1)τ) の状態を表示する */ glBegin(GL_LINE_LOOP); glVertex2d(X(0.0), Y(u[0])); for (i = 1; i <= N; i++) glVertex2d(X(i * h), Y(u[i])); glEnd(); } glFlush(); } int main(int argc, char **argv) { /* N, λ を入力する */ printf("区間の分割数 N = "); scanf("%d", &N); printf("λ (=τ/h^2) = "); scanf("%lf", &lambda); /* 最終時刻を入力する */ printf("最終時刻 Tmax = "); scanf("%lf", &Tmax); printf("計算終了後、プログラムを終了するには q または ESC\n"); /* ***************** グラフィックスの準備 ***************** */ glutInitWindowPosition(0, 0); glutInitWindowSize(WIDTH, HEIGHT); glutInit(&argc, argv); glutInitDisplayMode(GLUT_RGBA); glutCreateWindow(argv[0]); glutDisplayFunc(display); glutReshapeFunc(resize); glutKeyboardFunc(keyboard); init(); glutMainLoop(); } double f(double x) { if (x <= 0.5) return x; else return 1.0 - x; }