1

我正在编写一段代码来解决具有电荷密度边界的拉普拉斯算子，用 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、alpha 和 80.0，以及用于计算电荷和 mat[i][j] 的各行中的符号，并根据我输入的内容得到截然不同的结果：例如，将最后一行的电荷系数 [i] 更改为只有 10.0 会导致程序陷入循环，发散到无穷大或 -infinity，或快速收敛到 0（这是意料之外的）。这表明我创建的代码存在问题。

我还尝试将计算浓缩为一行，或者在不同的行中执行相同的步骤，这些也会极大地改变结果。

对此的任何帮助将不胜感激。谢谢。

nb 所有数据类型都是双精度的

编辑 - 完整的循环看起来像：

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);
4

0 回答 0