この記事は、Yılmaz Yörü氏のhttps://blogs.embarcadero.com/solve-millions-of-unknowns-in-equations-with-c-we-show-you-how/ の抄訳です
C++は、エンジニアリングの問題を計算するための優れたプログラミング言語であり、これらの操作に最適な言語の1つです。 Delphiもまた、エンジニアリングの問題に使用できる、より高速なコンパイルされたプログラミング言語です。このブログでは、SOR反復法を使用して、数百万の未知数を持つ方程式を解く方法を説明します。
まず、方程式の中の未知数の科学について 数学には「陰伏方程式」というものがありますが、私たちエンジニアは、現実の問題に関する数値方程式を解くことに直面しています。熱力学、流体力学、機構論、構造解析など、多くの工学分野で未知の方程式を解くことが求められます。一般的にこれらの陰伏方程式は、1次元、2次元、3次元のグリッドの例(例：非圧縮流体の流れ、熱分布の圧縮)で公開されており、これらのグリッドにはパラメータ（例：速度、圧力、温度など)を持つノードがあります。
数学において陰伏方程式は、R(x1,…, xn)=0の形に表される関係によって、Rは複数の変数の関数です。例えば、単位円の陰伏方程式は、x2 + y2 – 1 = 0です。より複雑な方程式は、熱力学の第一法則と第二法則、運動量の方程式、ナビエ–ストークス方程式などを解く方程式として与えられます。
陰関数とガウス＝ザイデル法 以下の内容は、陰関数 に関するwikipedia からの引用の一部です。
数学 の特に解析学 における陰函数 （いんかんすう、英 : implicit function ; 陰伏函数 ）は、陰伏方程式 すなわち適当な多変数函数（ しばしば多変数多項式 ）R によって R (x 1, …, xn ) = 0 の形に表される関係 によって（その函数の引数 のうちの一つの変数 の値（英語版） を残りの変数に関係付けることによって）陰伏的 (implicitly) に定義される函数 を言う[1] :204–206 。
これらの陰伏方程式は、反復法を用いて解くことができます。最もポピュラーな方程式の一つがSOR(Successive Over Relaxation Iteration )反復法です。この方法は、1次元、2次元、3次元の問題を解決するために使用できます。数値線形代数では、SOR（Successive Over-Relaxation）反復法 (以下、SOR反復法)は、線形連立方程式を解くためのガウス＝ザイデル法 の一種で、収束が速くなります。同様の方法は、ゆっくりと収束する反復プロセスにも用いることができます。
1次元のグリッド問題：3対角線の行列を解くには? 1次元グリッド上の陰伏方程式は、三角行列形式を生成します。問題が1次元のグリッドにある場合（例えば、ワイヤー上の熱伝導が1次元の問題である場合）、これは3つの対角行列形式が生成されます。なぜなら、このワイヤーの各1次元グリッドのノードは、最初と最後のノードを除いて、左ノードと右ノードに結ばれているからです。つまり、2次元行列の各行には3つの変数（左、中、右）が存在し、解決したい3つの未知数に影響を与えます。
行列形式の簡単な方程式は、以下のように記述できます。
この方程式では、
A は2次元の行列 (2本の横棒は2次元の行列であることを示す)U は各ノードからの未知のパラメータの1次元行列形式q は式の右辺の1次元行列形式三角柱のA行列、U行列、q行列は、ここでは以下のように示されます。
2次元のグリッド問題：5対角線の行列を解くには？ 問題が2次元グリッド上にある場合（例えば、プレート上の熱伝導は2次元問題である場合）、これらの方程式は5つの対角行列形式を生成します。なぜなら、この平面上の各2次元グリッドノードは、コーナーとエッジを除いて、左、右、上、下のノード（または、東、西、北、南と呼ばれます）に結ばれているからです。つまり、2次元行列の各行には最大5つの変数が存在することになります。
この方程式では、
A は2次元の行列 (2本の横棒は2次元の行列であることを示す)U は各ノードからの未知のパラメータの1次元行列形式q は式の右辺の1次元行列形式5つの対角線上にあるA行列、U行列、q行列は、以下のように示されます。
3次元のグリッドの例: 7対角線の行列の解き方 上記の例と同様に、問題が3次元グリッド上にある場合（例えば、立方体上の熱伝導が3次元問題である場合）、7つの対角行列形式が生成されます。なぜなら、この立方体の各次元ググリッドのノードは、面、コーナー、エッジを除いて、左、右、前、後ろ、上、下のノード（または東、西、北、南、上、下と呼ばれます）に結ばれているからです。 つまり、2次元行列の各行には最大7つの変数が存在することになります。
SOR反復法による3対角行列の解法例 この方法を理解するために、簡単に3つの対角線上の行列配列を使用してみましょう。新しいC++Builder VCLプロジェクトを作成し、「Create」と「Solve」というキャプションのついた2つのボタン（TButton）を追加し、メモ（TMemo）を追加して出力を簡単に確認します。
次にノード数を定義しましょう。NNは10で、行列配列を作成しましょう。m[] []は2次元行列（上記の方程式では、Aに該当）、u []は未知の値、x []は右側の行列です（ 上記の方程式では、qに該当）、そして最後に、反復で解を保持するためにsol []を定義するので、VCLプロジェクトでそれらを定義しましょう。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop
#include "SOR_Iteration_Unit1b.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
#define max(__a,__b) (((__a) > (__b)) ? (__a) : (__b))
#define min(__a,__b) (((__a) < (__b)) ? (__a) : (__b))
#define NN 10 // Number of Nodes
TForm1 * Form1 ;
double m [ NN + 1 ] [ NN + 1 ] , u [ NN + 1 ] , x [ NN + 1 ] , sol [ NN + 1 ] ;
//---------------------------------------------------------------------------
__fastcall TForm1 :: TForm1 ( TComponent * Owner ) : TForm ( Owner )
{
}
例えば、未知数が判明していて、解が1,2,3,4…と仮定しましょう。これらの解を使用して、以下のような対角形式でランダムな方程式を作ることができます。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
void generate_random_3diag_equations ( )
{
int i , j ;
srand ( time ( 0 ) ) ;
for ( i = 0 ; i < NN ; i ++ ) u [ i ] = i ;
// First Line, Top Left Corner
m [ 0 ] [ 0 ] = ( 10 + rand ( ) % 20 ) / 100.0 ; // central node
m [ 1 ] [ 0 ] = ( 10 + rand ( ) % 20 ) / 100.0 ; // right node
x [ 0 ] = u [ 0 ] * m [ 0 ] [ 0 ] + u [ 1 ] * m [ 1 ] [ 0 ] ; // u0*E + u1*F
for ( j = 1 ; j < NN - 1 ; j ++ )
{
m [ j - 1 ] [ j ] = ( 10.0 + rand ( ) % 20 ) / 100.0 ; // left node
m [ j ] [ j ] = ( 50.0 + rand ( ) % 10 ) / 100.0 ; // central node
m [ j + 1 ] [ j ] = ( 10.0 + rand ( ) % 20 ) / 100.0 ; // right node
x [ j ] = u [ j - 1 ] * m [ j - 1 ] [ j ] + u [ j ] * m [ j ] [ j ] + u [ j + 1 ] * m [ j + 1 ] [ j ] ;
}
// Last Line, Bottom Right Corner
//j = NN-1;
m [ j - 1 ] [ j ] = ( 10 + rand ( ) % 50 ) / 100.0 ; // left node
m [ j ] [ j ] = ( 10 + rand ( ) % 50 ) / 100.0 ; // central node
x [ j ] = u [ j - 1 ] * m [ j - 1 ] [ j ] + u [ j ] * m [ j ] [ j ] ;
}
このプロシージャでは、
FloatToStrF()のフォーマットを使用することで、すべての行列を表示することができます。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
void display_matrixes ( )
{
String str ;
for ( int j = 0 ; j < NN ; j ++ )
{
str = "" ;
for ( int i = 0 ; i < NN ; i ++ )
{
str = str + FloatToStrF ( m [ i ] [ j ] , TFloatFormat :: ffFixed , 9 , 3 ) + " " ;
}
Form1 -> Memo1 -> Lines -> Add ( "|" + str + "| |"
+ FloatToStrF ( u [ j ] , TFloatFormat :: ffFixed , 9 , 3 ) + "| |"
+ FloatToStrF ( x [ j ] , TFloatFormat :: ffFixed , 9 , 3 ) + "| |"
+ FloatToStrF ( sol [ j ] , TFloatFormat :: ffFixed , 9 , 3 ) + "|"
) ;
}
}
SOR反復法とは？ SOA反復法については、こちら(英語) で詳しく説明されています。この方法にはいくつかの定数があります。
max iteration : 計算を終了するまでの最大反復回数は1000回に設定omega : relaxation factorはここでは、0.5に設定min_error : 最小限の誤差で済むように、0.0001に設定C++Builderでは、以下の手順で行列形式を作成することができます。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
unsigned int calc_SOR ( const unsigned int max _ iteration )
{
const double omega = 0.5 ;
const double min_error = 0.0001 ;
double err ;
double prev , next , psol ;
unsigned int iteration = 0 ;
do
{
err = 0 ;
for ( int j = 0 ; j < NN ; j ++ )
{
prev = 0 ; next = 0 ;
for ( int i = 0 ; i < j ; i ++ ) prev += m [ i ] [ j ] * sol [ i ] ; // these lines costly
for ( int i = j + 1 ; i < NN ; i ++ ) next += m [ i ] [ j ] * sol [ i ] ; // these lines costly
psol = sol [ j ] ;
sol [ j ] = ( 1.0 - omega ) * u [ j ] + omega * ( x [ j ] - prev - next ) / m [ j ] [ j ] ;
err += pow ( sol [ j ] - psol , 2 ) ;
}
err = err / NN ;
for ( int i = 0 ; i < NN ; i ++ ) u [ i ] = sol [ i ] ;
iteration ++ ;
Form1 -> Memo1 -> Lines -> Add ( UIntToStr ( iteration ) + ":" + FloatToStrF ( err , TFloatFormat :: ffFixed , 12 , 6 ) ) ;
} while ( min_error < err && iteration < max_iteration ) ;
return iteration ;
}
ここまでで紹介してきた上記の例では、完全な行列を使用していますが、これはより大きなメモリ割り当てる目的では適していません。 以下のこの行列形式を使用して、1次元、2次元、3次元の問題のすべての対角線を計算できます。
struct matrix_line
{
float A , B , D , E , F , H , I ; // values in a line of 7 diagonal matrix form to solve 1D, 2D and 3D problems
short int a , b , d , f , h , i ; // X positions (index) of A,B,D,E,F,H,I avlues
} ;
struct matrix_line mx [ NN + 1 ] ;
しかし、この方法を理解するために、まず完全な行列形式を使用してみましょう。構造体ベースの行列形式については、別の記事で説明します。
最後に、「Create」と「Solve」ボタンの
OnClick()イベントでこれらのメソッドを使用できます。 それでは、3つの対角行列を作成し、以下のようにメモに表示してみましょう。
C++を使用した3対角行列の例
void __fastcall TForm1 :: Button1Click ( TObject * Sender )
{
generate_random_3diag_equations ( ) ;
Memo1 -> Lines -> Add ( "Generated Equations in Matrix Form:" ) ;
display_matrixes ( ) ;
}
未知数をリセットして（解が1,2,3,4…であることはわかっている）、SOR反復法を使用してこれらの方程式を解きましょう。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
void __fastcall TForm1 :: Button2Click ( TObject * Sender )
{
Memo1 -> Lines -> Add ( "Reset Unknowns:" ) ;
for ( int i = 0 ; i < NN ; i ++ ) u [ i ] = 0 ;
display_matrixes ( ) ;
Memo1 -> Lines -> Add ( "Calculating..." ) ;
try
{
calc_SOR ( 1000 ) ;
}
catch ( . . . )
{
} ;
Memo1 -> Lines -> Add ( "Solution:" ) ;
display_matrixes ( ) ;
}
上記のサンプルプログラムを実行し、その出力結果は以下のようになります。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
Generated Equations in Matrix Form :
| 0.210 0.250 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | | 0.000 | | 0.250 | | 0.000 |
| 0.160 0.500 0.100 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | | 1.000 | | 0.700 | | 0.000 |
| 0.000 0.280 0.570 0.260 0.000 0.000 0.000 0.000 0.000 0.000 | | 2.000 | | 2.200 | | 0.000 |
| 0.000 0.000 0.240 0.560 0.120 0.000 0.000 0.000 0.000 0.000 | | 3.000 | | 2.640 | | 0.000 |
| 0.000 0.000 0.000 0.220 0.500 0.230 0.000 0.000 0.000 0.000 | | 4.000 | | 3.810 | | 0.000 |
| 0.000 0.000 0.000 0.000 0.140 0.540 0.160 0.000 0.000 0.000 | | 5.000 | | 4.220 | | 0.000 |
| 0.000 0.000 0.000 0.000 0.000 0.240 0.570 0.200 0.000 0.000 | | 6.000 | | 6.020 | | 0.000 |
| 0.000 0.000 0.000 0.000 0.000 0.000 0.230 0.520 0.270 0.000 | | 7.000 | | 7.180 | | 0.000 |
| 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.100 0.540 0.150 | | 8.000 | | 6.370 | | 0.000 |
| 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.220 0.220 | | 9.000 | | 3.740 | | 0.000 |
Reset Unknowns :
| 0.210 0.250 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | | 0.000 | | 0.250 | | 0.000 |
| 0.160 0.500 0.100 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | | 0.000 | | 0.700 | | 0.000 |
| 0.000 0.280 0.570 0.260 0.000 0.000 0.000 0.000 0.000 0.000 | | 0.000 | | 2.200 | | 0.000 |
| 0.000 0.000 0.240 0.560 0.120 0.000 0.000 0.000 0.000 0.000 | | 0.000 | | 2.640 | | 0.000 |
| 0.000 0.000 0.000 0.220 0.500 0.230 0.000 0.000 0.000 0.000 | | 0.000 | | 3.810 | | 0.000 |
| 0.000 0.000 0.000 0.000 0.140 0.540 0.160 0.000 0.000 0.000 | | 0.000 | | 4.220 | | 0.000 |
| 0.000 0.000 0.000 0.000 0.000 0.240 0.570 0.200 0.000 0.000 | | 0.000 | | 6.020 | | 0.000 |
| 0.000 0.000 0.000 0.000 0.000 0.000 0.230 0.520 0.270 0.000 | | 0.000 | | 7.180 | | 0.000 |
| 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.100 0.540 0.150 | | 0.000 | | 6.370 | | 0.000 |
| 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.220 0.220 | | 0.000 | | 3.740 | | 0.000 |
Calculating . . .
1 : 14.927737
2 : 1.204483
3 : 0.110283
4 : 0.014905
5 : 0.003814
6 : 0.001554
7 : 0.000785
8 : 0.000444
9 : 0.000270
10 : 0.000172
11 : 0.000114
12 : 0.000077
Solution :
| 0.210 0.250 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | | 0.094 | | 0.250 | | 0.094 |
| 0.160 0.500 0.100 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | | 0.944 | | 0.700 | | 0.944 |
| 0.000 0.280 0.570 0.260 0.000 0.000 0.000 0.000 0.000 0.000 | | 2.067 | | 2.200 | | 2.067 |
| 0.000 0.000 0.240 0.560 0.120 0.000 0.000 0.000 0.000 0.000 | | 2.953 | | 2.640 | | 2.953 |
| 0.000 0.000 0.000 0.220 0.500 0.230 0.000 0.000 0.000 0.000 | | 4.031 | | 3.810 | | 4.031 |
| 0.000 0.000 0.000 0.000 0.140 0.540 0.160 0.000 0.000 0.000 | | 4.996 | | 4.220 | | 4.996 |
| 0.000 0.000 0.000 0.000 0.000 0.240 0.570 0.200 0.000 0.000 | | 5.988 | | 6.020 | | 5.988 |
| 0.000 0.000 0.000 0.000 0.000 0.000 0.230 0.520 0.270 0.000 | | 7.018 | | 7.180 | | 7.018 |
| 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.100 0.540 0.150 | | 7.990 | | 6.370 | | 7.990 |
| 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.220 0.220 | | 9.013 | | 3.740 | | 9.013 |
この例では、12回の繰り返しで、0, 1, 2, 3, 4…の結果を得ることができ、誤差は0.000077でした。SOR反復法を使用して、グリッドシステム上の何百万もの未知数、つまり流体力学、熱力学、応力解析、さらに他の多くの線形方程式の問題を解決したりできます。
反復法について 反復法は、いくつかの問題を解決するのに適していますが、動作が収束反復から発散反復に移行するため、理解が難しく、問題を解決するのが難しい場合もあります。よってそれが原因で間違った解決策を引き起こす可能性があります。そのため、SOR反復法を知っているだけでは十分ではなく、緩和因子、エラー係数、仮定、行列の値など、その挙動についてよく知っておく必要があります。結果を出すため、あるいは反復回数を少なくするために、無次元パラメータや再構成された行列を使わなければならないこともあります。
いくつかのランダム生成テストでは、反復回数の誤差が収束ではなく発散することがわかっています。これらの問題を解決するためのSOR改訂版やその他の新しいメソッドがあります。より良い解決策については、学術論文や最新のリリースを確認してください。このような問題を解決するためには、数学やプログラミングの知識だけでは不十分です。そのため数学やプログラミングに精通するだけではなく、これらのメソッドの動作も知っておく必要があります。
