C言語で、電荷密度との境界でラプラシアンを解くコードを書いています。次のループを使用しています。
chargeold[i] = charge[i];
charge[i] = -0.05*sgn(mat[i][j])*mat[i][j]*mat[i][j];
charge[i] = (1.0 - alpha)*charge[i] + alpha*chargeold[i];
mat[i][j] = (( mat[i][j-1] + di2*mat[i][j+1] ) / ( di2 + 1.0)) + 80.0*charge[i];
ここで、定数 alpha は緩和不足のパラメータ 0 < alpha < 1 です。 0.05、アルファ、80.0、および電荷とマット[i] [j]を計算するためのさまざまな行の符号と、私が何を入れたかによって大きく異なる結果が得られます。たとえば、最後の行の電荷[i]の係数を次のように変更します10.0 だけでは、プログラムがループ内に閉じ込められたり、無限大または無限大に発散したり、すぐに 0 に収束したりする可能性があります (これは予期されていません)。これは、作成したコードに問題があることを示唆しています。
また、計算を 1 行にまとめたり、別の行で同じ手順を実行したりしてみましたが、これらも結果を大幅に変更します。
これに関する助けをいただければ幸いです。ありがとう。
nb すべてのデータ型は double です
編集 - 完全なループは次のようになります。
do
{
sum_matdiff = 0;
for (i = 1; i < meshno; i++)
{
for (j = 1; j < meshno; j++)
{
if (bound[i][j] == 1) // holds boundary conditions
continue;
else
matold[i][j] = mat[i][j];
if ((i + j + count) % 2 == 1)
{
continue;
}
else if (j == (int)(0.3 * meshno))
{ // if statement to calculate at boundary
chargeold[i] = charge[i];
charge[i] = -0.05 * sgn(mat[i][j]) * mat[i][j] * mat[i][j];
charge[i] = (1.0 - alpha) * charge[i] + alpha * chargeold[i];
mat[i][j] =
((mat[i][j - 1] + di2 * mat[i][j + 1]) / (di2 + 1.0)) +
80.0 * charge[i];
if (i == 50)
{
printf("%f\n", charge[50]);
}
}
else
{ // calculates outside boundary
omega = 1.0 / (1.0 - 0.25 * omega * rho_sq);
mat[i][j] =
0.25 * (mat[i + 1][j] + mat[i - 1][j] + mat[i][j + 1] +
mat[i][j - 1] + (mat[i + 1][j] - mat[i - 1][j]) / (2 * i));
}
mat[i][j] = (1.0 - omega) * matold[i][j] + omega * mat[i][j];
sum_matdiff += fabs(1.0 - matold[i][j] / mat[i][j]);
}
}
count += 1;
av_diff = sum_matdiff / N;
}
while (av_diff > 0.01 || count < meshno * 2);