問題

[future]

#include <stdio.h>
#include <math.h>

/*
** 計算領域の分割数 NX, NY を決めます。理論上は分割数を細かくすれば
** するほど計算誤差は減るはずですが、コンピュータが扱う浮動小数点数
** は有限値なので、あまり細かくしすぎると、桁落ちによって逆に
** 計算精度が落ちると思います。いろいろな DX, DY で結果を比較して、
** 計算精度が充分に良いところで、かつ計算速度も妥当な分割数を選ぶと
** 良いと思います。ここでは仮に 41 x 41 とします。
** NX と NY は同じ数にした方が、微分方程式の特性に合うらしいです。
*/

#define NX 41
#define NY 41


/*
** for (y のループ) { for (x のループ) { (x 方向の計算); } }
** とする場合は、Ｃ言語の配列のメモリ配置を考慮して、T(x, y)
** ではなく、T(y, x) とします。
*/
double T[NY][NX];


/*-----------------------------------------------------------------------------
** [SetBoundaryCondition]
** 境界条件を設定します。
** ２次元配列の扱い方はいろいろあります。書籍や web page 等を
** あたってみてください。
**---------------------------------------------------------------------------*/
void SetBoundaryCondition(double T[NY][NX])
{
int x;
int y;


/* 領域の底1の境界値 0.5 */
if(y=NY/2){
for (x = 0; x < NX/3; x++)
{
T[0][x] = 0.0;
}
}

/* 領域の底2の境界値 0.5 */
for (x = NX/3; x < NX; x++)
{
T[0][x] = 50.0;
}

/* 領域の上の境界値は 0.0 */
for (x = 0; x < NX; x++)
{
T[NY-1][x] = 100.0;
}

/* 左1の壁の境界値は 1.0 */
if(x=NX/3){
for (y = 0; y < NY/2; y++)
{
T[y][0] = 0.0;
}
}

/* 左2の壁の境界値は 1.0 */
for (y = NY/2; y < NY; y++)
{
T[y][0] = 0.0;
}

/* 右の壁の境界値は 0.0 */
for (y = 0; y < NY; y++)
{
T[y][NX-1] = 100.0;
}
}


int main(void)
{
const double dx = 2.0 / (NX-1);
const double dy = 1.0 / (NY-1);

/* 式(16)から、C1 と C2 の値を先に計算しておきます。 */
const double C1 = 0.5 * dx*dx / (dx*dx + dy*dy);
const double C2 = 0.5 * dy*dy / (dx*dx + dy*dy);

int i;
int j;

FILE *output;
output=fopen("output.data","w");

double errorMax;
double previousValue;

/*
** 境界条件を設定します。この境界条件は
** ディリクレ条件なので、はじめに１度だけ設定
** すればＯＫです。
*/
SetBoundaryCondition(T);

/*
** 前回の計算値と今回の計算値との最大誤差が
** 0.00001 を下回るまで繰り返します。
*/
do
{
errorMax = 0.0;

/* 計算領域は境界の１つ内側です。 */
for (i = 1; i < NY-1; i++)
{
for (j = 1; j < NX-1; j++)
{
previousValue = T[j];

/* ガウス・ザイデル法 */
T[j] = C1*(T[i+1][j] + T[i-1][j]) + C2*(T[j+1] + T[j-1]);

if (errorMax < fabs(T[j] - previousValue))
{
errorMax = fabs(T[j] - previousValue);
}
}
}
} while (errorMax > 0.00001);

/* 全領域の計算結果を出力します。 */
for (i = 0; i < NY; i++)
{
for (j = 0; j < NX; j++)
{
fprintf(output,"%10.15f\t%10.15f\t%15.15f\n", j*0.1 / (NX-1), i*0.1 / (NY-1), T[j]);
}
fprintf(output,"\n");
}

fclose(output);

return 0;
}


50505050505050505050505050505050
00000000000000000000000000000001
00000000000000000000000000000001
00000000000000000000000000000001
00000000000000000000000000000001
0000000000000000000001
0000000000000000000001
1111111111111111111111

そんな形状にしたいですけど、四角形状の結果が出ました。よろしくお願いします。

[softya（ソフト屋)]

前のトピックを放置して新たな質問はご遠慮ください。
「幾何形状を書く • C言語交流フォーラム　～ mixC++ ～」
http://dixq.net/forum/viewtopic.php?f=3&t=11304

それとフォーラルルールをご覧ください。
http://dixq.net/board/board.html
ソースコードの添付などに使うcodeタグの使い方が書かれています。

こちらは閉鎖しておきます。