2次元FDM用SOR法(C)プログラム
2003/08/24 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してSOR法で
反復解xを求める。
加速係数(ω)を与える必要がある。収束速度はωにより大きく左右される。
2. プログラム
#include <stdio.h>
#include <math.h>
// Global Define
#define NDX 301
#define NDY 301
#define A(i,j,l) A[l-1][j][i]
#define B(i,j) B[j][i]
#define X(i,j) X[j][i]
extern double A[5][NDY][NDX], B[NDY][NDX], X[NDY][NDX] ;
extern FILE *FT1 ;
//=================================================================C
int SOR2D(int NX,int NY,double Omega,double EPS,int *ITER,double *ERR)
//=================================================================C
// Solve Ax=b by SOR with 2 dimensional FDM C
// Given Omega ( Acceleration factor ) C
//-----------------------------------------------------------------C
// A(0:NX,0:NY,5) R*8, In, A Coefficient Matrix C
// B(0:NX,0:NY) R*8, In, A Right-hand Vector(b) C
// NX I*4, In, Grid Numbers on X-axis C
// NY I*4, In, Grid Numbers on Y-axis C
// X(0:NX,0:NY) R*8, I/O, Initial and Solution vector C
// Omega R*8, In, SOR Acceleration factor C
// EPS R*8, In, if ||r||/||b|| <= EPS --> return C
// ITER I*4, I/O, Number of Iteration C
// ERR R*8, Out, ERR=||r||/||b|| C
// IERR I*4, Out, IERR=0, Normal Return C
// =1, No Convergent C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/08/21 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
{ int i, j, k, IERR ;
double R, BN, RN, RA ;
// Initial
IERR = 0 ;
BN = 0.0 ;
for (j=1; j<=NY-1; j++) {
for (i=1; i<=NX-1; i++) {
BN = BN + B(i,j)*B(i,j) ; }
}
// Main Loop
for (k=1; k<=*ITER; k++) {
RN = 0.0 ;
for (j=1; j<=NY-1; j++) {
for (i=1; i<=NX-1; i++) {
R = (B(i,j) - A(i,j,1)*X(i,j-1) - A(i,j,2)*X(i-1,j)
- A(i,j,4)*X(i+1,j) - A(i,j,5)*X(i,j+1) )
/ A(i,j,3) - X(i,j) ;
X(i,j) = X(i,j) + Omega*R ;
RA = R*A(i,j,3) ;
RN = RN + RA*RA ; }
}
// if(ERR <= EPS) Convergent
*ERR = sqrt(RN/BN) ;
if(*ERR <= EPS) goto M100 ;
}
IERR = 1 ;
M100: *ITER = k ;
return (IERR) ; }