差分法用反復解法プログラム

　　　　　　　　2003/08/24 日立製作所 & 早稲田大学　後 保範 ( Ushiro Yasunori )

1. 概要

　差分法による離散化で発生する規則的疎行列を係数とする連立一次方程式Ax=bを反復解法
で数値解xを計算するFORTRAN及びCのソースプログラムを掲載する。
本プログラムは内容を理解してもらうことを目的とし、計算時間及び使用メモリ量の効率化
は目的としていない。本格的なプログラムは金田研究室から他のプログラムと合わせて
掲載している反復解法は2次元及び3次元差分法用のSOR法及びCG法である。
公開する予定。配列の格納順はCの場合も＃define文を使用しFORTRANに合わせてある。

2. SOR法プログラム

(1) 2次元FDM用プログラム

　(a) 計算方法　： x及びy方向の分割数NX,NY及び加速係数(ω)を入力し、2次元差分法で離散化
　　した、連立一次方程式Ax=bの解xをSOR法(逐次過剰緩和法)で求める。
　　-div(ｋ.grad(u))=f, k=1, f=100を2次元差分法で離散化して疎行列を作成している。
　　行列は5行の疎行列で左からA1,A2,A3(対角),A4,A5の順である。
　　計算結果は、収束状況とy軸の中心の各x-格子点での解xの値を出力する。
　(b) プログラム ： FORTRAN とCのソースプログラム
　(c) テストプログラム : FORTRAN とCのソースプログラム

(2) 3次元FDM用プログラム

　(a) 計算方法　： x,y,z方向の分割数NX,NY,NZ及び加速係数(ω)を入力し、3次元差分法で離散化
　　した、連立一次方程式Ax=bの解xをSOR法(逐次過剰緩和法)で求める。
　　-div(ｋ.grad(u))=f, k=1, f=100を3次元差分法で離散化して疎行列を作成している。
　　行列は7行の疎行列で左からA1,A2,A3,A4(対角),A5,A6,A7の順である。
　　計算結果は、収束状況とy-z軸の中心の各x-格子点での解xの値を出力する。
　(b) プログラム ： FORTRAN とCのソースプログラム
　(c) テストプログラム : FORTRAN とCのソースプログラム

3. CG法プログラム

(1) 2次元FDM用プログラム

　(a) 計算方法　： x及びy方向の分割数NX,NY及び加速係数(ω)を入力し、2次元差分法で離散化
　　した、連立一次方程式Ax=bの解xをCG法(共役勾配法)で求める。
　　-div(ｋ.grad(u))=f, k=1, f=100を2次元差分法で離散化して疎行列を作成している。
　　行列は5行の疎行列で、対称行列のため下三角部分を持ち左からA1,A2,A3(対角)順である。
　　計算結果は、収束状況とy軸の中心の各x-格子点での解xの値を出力する。
　　FORTRANプログラムは結果のy軸対称性を調べるため、y軸対称とならない結果を全て出力できる。
　(b) プログラム-1 ： α=(r,p)/(p,Ap),β=-(r,Ap)/(p,Ap)の計算で FORTRAN とCのソースプログラム
　(c) プログラム-2 ： α=(r,r)/(p,Ap),β=new(r,r)/old(r,r)ので FORTRAN とCのソースプログラム
　(d) テストプログラム : FORTRAN とCのソースプログラム

(2) 3次元FDM用プログラム

　(a) 計算方法　： x,y,z方向の分割数NX,NY,NZ及び加速係数(ω)を入力し、3次元差分法で離散化
　　した、連立一次方程式Ax=bの解xをCG法(共役勾配法)で求める。
　　-div(ｋ.grad(u))=f, k=1, f=100を3次元差分法で離散化して疎行列を作成している。
　　行列は7行の疎行列で、対称行列のため下三角部分を持ち左からA1,A2,A3,A4(対角)順である。
　　計算結果は、収束状況とy-z軸の中心の各x-格子点での解xの値を出力する。
　(b) プログラム-1 ： α=(r,p)/(p,Ap),β=-(r,Ap)/(p,Ap)の計算で FORTRAN とCのソースプログラム
　(c) プログラム-2 ： α=(r,r)/(p,Ap),β=new(r,r)/old(r,r)の計算で FORTRAN とCのソースプログラム
　(d) テストプログラム : FORTRAN とCのソースプログラム