我们正在尝试将 FORTRAN 程序转换为 C。这是一个 SOR 实现，它应该在大约 200 次迭代时收敛，但 C 代码不起作用。它在大约第 100 次迭代时变得“奇怪”(它开始非常缓慢地收敛)。
原始代码如下:
PROGRAM Main
! *** Solution of Laplace's Equation.
! ***
! *** Uxx + Uyy = 0
! *** 0 <= x <= pi, 0 <= y <= pi
! *** U(x,pi) = sin(x), U(x,0) = U(0,y) = U(pi,y) = 0
! ***
! *** then U(x,y) = (sinh(y)*sin(x)) / sinh(pi)
! ***
! *** Should converge with
! *** tol = 0.001 and M = 20 in 42 iterations.
! *** and with tol = 0.001 and M = 100 in 198 iterations.
! ***
INTEGER M
PARAMETER (M = 100)
REAL*8 PI
REAL*8 DATAN
REAL*8 unew(0:M+1,0:M+1), uold(0:M+1,0:M+1)
REAL*8 solution(0:M+1,0:M+1)
REAL*8 omega, tol, h
INTEGER i, j, iters
PI = 4.0D0*DATAN(1.0D0)
h = PI/FLOAT(M+1)
DO i=0,M+1
uold(i,M+1) = SIN(FLOAT(i)*h)
END DO
DO i=0,M+1
DO j=0,M
uold(i,j) = FLOAT(j)*h*uold(i,M+1)
END DO
END DO
DO i=0,M+1
DO j=0,M+1
solution(i,j) = SINH(FLOAT(j)*h)*SIN(FLOAT(i)*h)/SINH(PI)
END DO
END DO
omega = 2.0/(1.0+SIN(PI/FLOAT(M+1)))
tol = 0.001
CALL SOR (unew, uold, solution, omega, tol, m, iters)
PRINT *, " "
PRINT *, " Omega = ", omega
PRINT *, " It took ", iters, " iterations."
STOP
END PROGRAM Main
REAL*8 FUNCTION ComputeError (solution, u, m, iters)
INTEGER m, iters
REAL*8 u(0:m+1,0:m+1)
REAL*8 solution(0:m+1,0:m+1)
! *** Local variables
REAL*8 error
INTEGER i, j
error = 0.0
DO j=1,m
DO i=1,m
error = MAX(error, ABS(solution(i,j)-u(i,j)))
END DO
END DO
PRINT *, "On iteration ", iters, " error = ", error
ComputeError = error
RETURN
END FUNCTION ComputeError
SUBROUTINE SOR (unew, uold, solution, omega, tol, m, iters)
INTEGER m, iters
REAL*8 unew(0:m+1,0:m+1), uold(0:m+1,0:m+1)
REAL*8 solution(0:m+1,0:m+1)
REAL*8 omega, tol
! *** Local variables
REAL*8 error
INTEGER i, j
! *** External function.
REAL*8 ComputeError
EXTERNAL ComputeError
! *** Copy bonudary conditions.
DO i=0,m+1
! print *,i
unew(i,m+1) = uold(i,m+1)
unew(m+1,i) = uold(m+1,i)
unew(i, 0) = uold(i, 0)
unew(0, i) = uold(0, i)
END DO
! *** Do SOR until 'tol' satisfied.
iters = 0
error = ComputeError (solution, uold, m, iters)
DO WHILE (error .GE. tol)
! *** Do one iteration of SOR
DO j=1,m
DO i=1,m
unew(i,j) = uold(i,j) + 0.25*omega* &
(unew(i-1,j) + unew(i,j-1) + &
uold(i+1,j) + uold(i,j+1) - &
4.0*uold(i,j))
END DO
END DO
! *** Copy new to old.
DO j=1,m
DO i=1,m
uold(i,j) = unew(i,j)
END DO
END DO
! *** Check error every 20 iterations.
iters = iters + 1
IF (MOD(iters,2) .EQ. 0) THEN
error = ComputeError (solution, uold, m, iters)
END IF
END DO
RETURN
END SUBROUTINE SOR
这是我们在 C 中得到的:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
#define M 100
double compute_error(double solution[][M + 2], double u[][M + 2], const int m, const int iters);
void sor(double unew[][M + 2], double uold[][M + 2], double solution[][M + 2], double omega, double tol, int m, int iters);
int main(int argc, char ** argv)
{
/*
Solution of Laplace's Equation.
==============================
! ***
! *** Uxx + Uyy = 0
! *** 0 <= x <= pi, 0 <= y <= pi
! *** U(x,pi) = sin(x), U(x,0) = U(0,y) = U(pi,y) = 0
! ***
! *** then U(x,y) = (sinh(y)*sin(x)) / sinh(pi)
! ***
! *** Should converge with
! *** tol = 0.001 and M = 20 in 42 iterations.
! *** and with tol = 0.001 and M = 100 in 198 iterations.
*/
double pi;
double unew[M + 2][M + 2];
double solution[M + 2][M + 2];
double uold[M + 2][M + 2];
double omega, tol, h;
int i, j, iters;
// Initialization
for(i = 0; i <= M + 1; ++i)
{
for(j = 0; j <= M + 1; ++j)
{
unew[i][j] = 0;
uold[i][j] = 0;
solution[i][j] = 0;
}
}
pi = 4.0 * atan(1.0);
h = pi / (M + 2);
for(i = 0; i <= M + 1; ++i)
{
uold[i][M + 1] = sin(i * h);
}
for(i = 0; i <= M + 1; ++i)
{
for(j = 0; j <= M; ++j)
{
uold[i][j] = j * h * uold[i][M + 1];
}
}
for(i = 0; i <= M + 1; ++i)
{
for(j = 0; j <= M + 1; ++j)
{
solution[i][j] = sinh(j * h) * sin(i * h) / sinh(pi);
}
}
omega = 2.0 / ( 1.0 + sin(pi / (M + 1)) );
printf("omega = %.10f\n", omega);
tol = 0.001;
sor(unew, uold, solution, omega, tol, M, iters);
printf(" \n");
printf(" Omega = %f\n", omega);
printf(" It took %d iterations.\n", iters);
return 0;
}
double compute_error(double solution[][M + 2], double u[][M + 2], const int m, const int iters)
{
double error = 0.0;
int i, j;
for(j = 1; j <= m; ++j)
{
for(i = 1; i <= m; ++i)
{
if(error >= fabs(solution[i][j] - u[i][j]))
error = error;
else
error = fabs(solution[i][j] - u[i][j]);
}
}
printf("On iteration %d error = %f\n", iters, error);
return error;
}
void sor(double unew[][M + 2], double uold[][M + 2], double solution[][M + 2], const double omega, double tol, int m, int iters)
{
double error;
int i, j;
for(i = 0; i <= m + 1; ++i)
{
unew[i][m + 1] = uold[i][m + 1];
unew[m + 1][i] = uold[m + 1][i];
unew[i][0] = uold[i][0];
unew[0][i] = uold[0][i];
}
iters = 0;
error = compute_error(solution, uold, m, iters);
while(error >= tol)
{
for(j = 1; j <= m; ++j)
{
for(i = 1; i <= m; ++i)
{
unew[i][j] = uold[i][j] + 0.25 * omega *
(unew[i - 1][j] + unew[i][j - 1] +
uold[i + 1][j] + uold[i][j + 1] -
4.0 * uold[i][j]);
}
}
for(j = 1; j <= m; ++j)
{
for(i = 1; i <= m; ++i)
{
uold[i][j] = unew[i][j];
}
}
iters += 1;
if(iters % 2 == 0)
{
error = compute_error(solution, uold, m, iters);
}
}
}
我们怀疑这可能与精度有关。
最佳答案
iters 有问题，您按值传递但看起来应该按引用传递(即作为 C 中的指针)。 (请注意，gcc 会为此生成警告 - 确保您使用 -Wall 进行编译以捕获此类问题。)
您似乎还误译了这一行:
h = pi / (M + 1);
它最有可能是:
h = pi / (M + 2);
通过此更改，我在 196 次迭代后收敛:
...
On iteration 192 error = 0.001293
On iteration 194 error = 0.001104
On iteration 196 error = 0.000930
Omega = 1.939676
It took 196 iterations.
LIVE DEMO (这还包括对 iters 错误的修复)。
当然可能还有其他bug有待修复...
关于c - SOR 程序不收敛(FORTRAN 到 C 的转换)，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/26548188/