fork(1) download
copy 
  1. /******************************************************************************
  2. ** [laplace.c]
  3. **	 ラプラス方程式を差分法により解きます。
  4. **
  5. **	 使い方： laplace.exe > result.txt
  6. ******************************************************************************/
  7. #include <stdio.h>
  8. #include <math.h>
  9.  
  10. /*
  11. ** 計算領域の分割数 NX, NY を決めます。理論上は分割数を細かくすれば
  12. ** するほど計算誤差は減るはずですが、コンピュータが扱う浮動小数点数
  13. ** は有限値なので、あまり細かくしすぎると、桁落ちによって逆に
  14. ** 計算精度が落ちると思います。いろいろな DX, DY で結果を比較して、
  15. ** 計算精度が充分に良いところで、かつ計算速度も妥当な分割数を選ぶと
  16. ** 良いと思います。ここでは仮に 41 x 41 とします。
  17. ** NX と NY は同じ数にした方が、微分方程式の特性に合うらしいです。
  18. */
  19.  
  20. #define NX 41
  21. #define NY 41
  22.  
  23.  
  24. /*
  25. ** for (y のループ) { for (x のループ) { (x 方向の計算); } }
  26. ** とする場合は、Ｃ言語の配列のメモリ配置を考慮して、T(x, y)
  27. ** ではなく、T(y, x) とします。
  28. */
  29. double T[NY][NX];
  30.  
  31.  
  32. /*-----------------------------------------------------------------------------
  33. ** [SetBoundaryCondition]
  34. ** 境界条件を設定します。
  35. ** ２次元配列の扱い方はいろいろあります。書籍や web page 等を
  36. ** あたってみてください。
  37. **---------------------------------------------------------------------------*/
  38. void SetBoundaryCondition(double T[NY][NX])
  39. {
  40. 	int x;
  41. 	int y;
  42. 	double left, top, right, bottom;
  43.  
  44. 	fprintf(stderr, "left top right bottom\n");
  45. 	scanf_s("%lf%lf%lf%lf", &left, &top, &right, &bottom);
  46.  
  47. 	/* 領域の底の境界値 */
  48. 	for (x = 0; x < NX; x++)
  49. 	{
  50. 		T[0][x] = bottom;
  51. 	}
  52.  
  53. 	/* 領域の上の境界値 */
  54. 	for (x = 0; x < NX; x++)
  55. 	{
  56. 		T[NY-1][x] = top;
  57. 	}
  58.  
  59. 	/* 左の壁の境界値 */
  60. 	for (y = 0; y < NY; y++)
  61. 	{
  62. 		T[y][0] = left;
  63. 	}
  64.  
  65. 	/* 右の壁の境界値 */
  66. 	for (y = 0; y < NY; y++)
  67. 	{
  68. 		T[y][NX-1] = right;
  69. 	}
  70. }
  71.  
  72.  
  73. int main(void)
  74. {
  75. 	const double dx = 1.0 / (NX-1);
  76. 	const double dy = 1.0 / (NY-1);
  77.  
  78. 	/* 式(16)から、C1 と C2 の値を先に計算しておきます。 */
  79. 	const double C1 = 0.5 * dx*dx / (dx*dx + dy*dy);
  80. 	const double C2 = 0.5 * dy*dy / (dx*dx + dy*dy);
  81.  
  82. 	int i;
  83. 	int j;
  84.  
  85. 	double errorMax;
  86. 	double previousValue;
  87.  
  88. 	/*
  89. 	** 境界条件を設定します。この境界条件は
  90. 	** ディリクレ条件なので、はじめに１度だけ設定
  91. 	** すればＯＫです。
  92. 	*/
  93. 	SetBoundaryCondition(T);
  94.  
  95. 	/*
  96. 	** 前回の計算値と今回の計算値との最大誤差が
  97. 	** 0.00001 を下回るまで繰り返します。
  98. 	*/
  99. 	do
  100. 	{
  101. 		errorMax = 0.0;
  102.  
  103. 		/* 計算領域は境界の１つ内側です。 */
  104. 		for (i = 1; i < NY-1; i++)
  105. 		{
  106. 			for (j = 1; j < NX-1; j++)
  107. 			{
  108. 				previousValue = T[i][j];
  109.  
  110. 				/* ガウス・ザイデル法 */
  111. 				T[i][j] = C1*(T[i+1][j] + T[i-1][j]) + C2*(T[i][j+1] + T[i][j-1]);
  112.  
  113. 				if (errorMax < fabs(T[i][j] - previousValue))
  114. 				{
  115. 					errorMax = fabs(T[i][j] - previousValue);
  116. 				}
  117. 			}
  118. 		}
  119. 	} while (errorMax > 0.00001);
  120.  
  121. 	/* 全領域の計算結果を出力します。 */
  122. 	for (i = 0; i < NY; i++)
  123. 	{
  124. 		for (j = 0; j < NX; j++)
  125. 		{
  126. 			printf("%d\t%d\t%15.15f\n", j, i, T[i][j]);
  127. 		}
  128. 		printf("\n");
  129. 	}
  130.  
  131. 	return 0;
  132. }
  133.  
  134. // gnuplot>cd 'C:\hoge'
  135. // gnuplot>set pm3d
  136. // gnuplot>splot 'result.txt'
  137.  
Not running #stdin #stdout 0s 0KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/GYFqu
language:
C++ (gcc 8.3)
created:
12 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!