連立方程式をSOR法で解くPythonコードと緩和係数のパラスタ

  • このエントリーをはてなブックマークに追加

連立方程式を反復法で解くSOR法はガウス・ザイデル法よりも高速になる可能性を持ちますが、緩和係数(加速パラメータ)の調整が必要です。ここではPythonを使ってSOR法を実装し、緩和係数のパラメータスタディと最適な緩和係数との比較を行います。

こんにちは。wat(@watlablog)です。ここではSOR法を使って連立方程式を解く方法を紹介します

SOR法の概要

SOR法の漸化式

SORSuccessive Over-Relaxation Method)とは、連立方程式を反復法を使って解く手法の一種です。この手法はガウス・ザイデル法に緩和係数\(\omega\)(Relaxation coefficient)を追加した漸化式を持ちます。

Advertisements

緩和係数\(\omega\)は加速パラメータAcceleration parameter)とも呼ばれますが、その名の通り反復計算の収束速度を変化させるパラメータです。

SOR法を理解するためには「Pythonで連立方程式をヤコビ法(反復法)で解く方法」と「Pythonで連立方程式をガウス・ザイデル法(反復法)で解く方法」で紹介したヤコビ法ガウス・ザイデル法の理解が必須となります。

基礎知識は上記2つの記事に記載していますので、ここではいきなり漸化式の説明から始めます。

式(1)にSOR法の漸化式を示します\(^{[1]}\)。

\[ x^{n+1}_{i} = (1-\omega)x^{n}_{i} + \omega \frac{1}{a_{ii}} \left ( b_{i} - \sum_{j=1}^{i-1} a_{ij}x^{n+1}_{j} - \sum_{j=i+1}^{N} a_{ij}x^{n}_{j} \right ) \tag{1} \]

この式をガウス・ザイデル法の時と同様に行列形式で表すと、対角行列\(\mathbf{D}\)と剰余行列\(\mathbf{R}\)を使って式(2)と変形できます。
上添字nはステップ数です。また、剰余行列\(\mathbf{R}\)一つを使って\(x^{n+1}_{i}\)と\(x^{n}_{i}\)の係数を表現できる理由もガウス・ザイデル法の時に考察しました。

\[ x^{n+1}_{i}=(1-\omega) x^{n}_{i} + \omega\mathbf{D}^{-1}_{i}(b_{i} - \mathbf{R}_{i}\mathbf{x}) \tag{2} \]

最適な緩和係数

SOR法の緩和係数\(\omega\)は自分で決める必要があり、値によってはガウス・ザイデル法よりも遅くなるため扱いが難しい手法です。

しかし、ヤコビ法の反復行列が既知であればスペクトル半径\(\rho\)から最適な緩和係数\(\omega_{opt}\)を求めることが可能です。

ヤコビ法の反復行列\(\mathbf{M}_{j}\)を式(3)に示します。

\[ \begin{eqnarray} \mathbf{M}_{j}&=&-\mathbf{D}^{-1}(\mathbf{L}+\mathbf{U}) \\ &=&-\mathbf{D}^{-1}\mathbf{R} \end{eqnarray} \tag{3} \]

スペクトル半径\(\rho\)は正方行列の最大固有値の絶対値\(^{[2]}\)を意味しており、式(3)の行列から得られる\(\rho\)を使った最適な緩和係数\(\omega_{opt}\)は式(4)で計算されます\(^{[3][4]}\)。

\[ \omega_{opt}= \frac{2}{1+\sqrt{1-\rho^{2}}} \tag{4} \]

ちなみに固有値はPythonを使えば簡単に求めることができます。是非「Pythonで線形代数!固有値と固有ベクトルを求める」の記事を参考にしてください。

SOR法で連立方程式を解くPythonコード

動作環境

このページでは以下の環境で動作確認を行ったコードを掲載しています。

Python環境

Python Python 3.9.6
PyCharm (IDE) PyCharm CE 2020.1
Numpy 1.21.1

PC環境

Mac OS macOS Catalina 10.15.7
CPU 1.4[GHz]
メモリ 8[GB]

まずは\(\omega\)を手入力

まずはシンプルさを優先して、緩和係数\(\omega\)を手入力して解くコードを以下に示します。ヤコビ法やガウス・ザイデル法とほぼ同じなので特に詳細な説明はいらないと思います。

以下が結果です。
なんとたった10回の反復回数で正解を近似しました。検証として定数ベクトル\(\mathbf{b}\)の確かめ算やnp.linalg.solveの結果も実施していますが問題はありませんでした。

スペクトル半径から\(\omega_{opt}\)を求めてパラスタ結果と比較する

最適な緩和係数を求める関数

式(3)、式(4)から最適な緩和係数を求めることができるとのことなので、本当なのかどうか試してみます。また、緩和係数というパラメータを色々振ることでパラスタ(パラメータスタディ)を行い、最適値と比較します。

まずは以下に示すように、最適な緩和係数を求める関数を書いてみました。

検証用の全コード

こちらが検証用の全コードです。せっかくなのでヤコビ法ガウス・ザイデル法にも登場していただき、matplotlibによるプロットで視覚的に確かめてみます。

下図が結果です。緩和係数\(\omega\)が1.0の時はガウス・ザイデル法と一致します。今回\(\omega\)は1.0から2.0まで0.1刻みで振りましたが、この連立方程式の場合は\(\omega=\)1.9, 2.0で発散(解がでたらめになる)するといった特徴もありました。

また最適な緩和係数の場合とパラスタで得られた\(\omega\)の値が近く、式(4)はしっかり実装できているのではと考えられます(文献があまりなかったので合っているか不安でした)。

反復法違いと緩和係数違いの比較結果

まとめ

この記事ではSOR法の概要として漸化式と最適な緩和係数加速パラメータ)を説明しました。最適な緩和係数は行列のスペクトル半径を使って求めるということで、固有値計算を使うことがわかりました。

筆者にとって線形代数は計算手法を覚えるということが主な学習方法でしたが、今回のスペクトル半径や固有値の関係性、行列操作の意味に重点を置いて学習するのも有用と感じました。

最後にSOR法をPythonで実装するコードを紹介しました。緩和係数にもよりますがこれまでのヤコビ法やガウス・ザイデル法よりも反復回数が削減できることがわかりました。

しかし、ヤコビ法は並列計算ができるという強みがあるので一概にどの手法が良いかは決まらないと思います。

参考文献

[1]:幸谷智紀, 連立一次方程式の解法2 -- 反復法 (2007), pp109 [2]:[線型代数] 特異値分解と低ランク近似 [3]:Varga, R. S. (2009). Matrix iterative analysis (Vol. 27). Springer Science & Business Media. [4]:山崎郭滋, 偏微分方程式の数値解法入門, 森北出版,(2015), pp28

ヤコビ法、ガウス・ザイデル法、SOR法の3つを使えるようになりました!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメントを残す

*