SOR法のC言語実装

参考文献

数値計算の基礎と応用[新訂版]
数値解析学への入門
杉浦　洋(南山大学教授)　著
発行日 2009/12/10

参考ページ

準備

オンラインコンパイラを使用します。

ソースコード

sample.c

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

#define Max(a,b) ((a>b)? a:b)
#define n2 100

double v(double x,double y){
    return exp(x)*cos(y);
}

double g(double x,double y){
    return v(x,y);
}

int main() {
    double u[n2][n2],Pi,d,eps,w,mu,r,rmax,ermax;
    int n,i,j,k,kmax;
    Pi=atan(1)*4;

    n=20;
    d=1.0/(n+1);
    w=2./(1+sin(Pi/(n+1)));
    mu=0;
    for(i=0;i<=n+1;i++){
        u[i][0]=g(i*d,0);
        u[i][n+1]=g(i*d,1);
        u[0][i]=g(0,i*d);
        u[n+1][i]=g(1,i*d);
        mu+=u[i][0]+u[i][n+1]+u[0][i]+u[n+1][i];
    }
    mu/=4*n;
    for(i=1;i<=n;i++){
        for(j=i;j<=n;j++){
            u[i][j]=mu;
        }
    }
    kmax=300;
    eps=1.0e-8;
    for(k=1;k<=kmax;k++){
        rmax=0;
        for(i=1;i<=n;i++){
            for(j=1;j<=n;j++){
                r=(u[i-1][j]+u[i+1][j]+u[i][j-1]+u[i][j+1])/4-u[i][j];
                rmax=Max(rmax,fabs(r));
                u[i][j]+=w*r;
            }
        }
        if(rmax<=eps)break;
    }

    printf("n=%d, eps=%9.2e, rmax=%9.2e\n",n,eps,rmax);
    printf("number of iterarions k=%d, (kmax=%d)\n",k,kmax);
    ermax=0;
    for(i=1;i<=n;i++){
        for(j=1;j<=n;j++){
            ermax=Max(ermax,fabs(u[i][j]-v(i*d,j*d)));
        }
    }
    printf("Maximum Error=%9.2e\n",ermax);
    return 0;
}

実行結果

console

n=20, eps= 1.00e-08, rmax= 7.82e-09
number of iterarions k=77, (kmax=300)
Maximum Error= 4.12e-05