工学系大学院生のブログ

2人の大学院生による雑記ブログ

第1-A回 pythonにおける対角行列の作り方

数値解析において用いる行列は対角行列に近いです。
対角成分とその前後の対角部分でできています。

今回の目標は非定常拡散方程式における境界条件も加味した下記の様なn行×n列の行列を作ることとします。

$$ \left[\begin{array}{cccccc}1&0& & &&\\ c&1-2c&c && &\\ & \ddots &\ddots &\ddots &&\\ && \ddots &\ddots& \ddots &\\ && &c&1-2c&c \\ && & &0&1 \end{array}\right] $$

先にコードを書いておきます。

import numpy
from scipy.sparse import dia_matrix,identity

""" matrix A """
data = numpy.array([numpy.ones(n), -2.0*numpy.ones(n), numpy.ones(n)])

data[2][1] = 0
data[0][n-2] = 0
data[1][0] = 0
data[1][n-1] = 0

offsets = numpy.array([-1, 0, 1])
B = dia_matrix((data, offsets), shape = (n, n))

E = identity(n)
A = E + c*B

対角行列を作るために,scipy.sparse.dia_matrixを使用します.
13行目でdataとoffsetsを与えています。これはdata[0]をoffsets[0]の列に配置します.
このoffsetsは行列の対角成分を位置0として入れます.下の行列において、bの部分が位置0となり、aは位置-1,cは位置1となります.

$$\left[\begin{array}{cccccc}b_0&c_1& & &&\\ a_0&b_1&c_2 && &\\ & \ddots &\ddots &\ddots &&\\ && \ddots &\ddots& \ddots &\\ && &a_{n-3}&b_{n-2}&c_{n-1} \\ && & &a_{n-2}&b_{n-1} \end{array}\right] $$

つまりdata[0],[1],[2]にそれぞれa, b, c の成分を入れればいいわけです.また\(a_{n-1}\) と\(c_0\)は行列に含まれないことに注意してください.


下記5行目では、n 個の要素の入った3つの列を作っています。

data = numpy.array([numpy.ones(n), -2.0*numpy.ones(n), numpy.ones(n)])


また次の4行でそれぞれ下記の通り与えています。

data[2][1] = 0
data[0][n-2] = 0
data[1][0] = 0
data[1][n-1] = 0

$$c_1=0 $$ $$ a_{n-2}=0 $$ $$ b_0=0 $$ $$b_{n-1}=0$$

そして、先ほど説明した dia_matrix を用いてB行列を作成します。この時点でB行列は次のようになっています。

$$B= \left[\begin{array}{cccccc}0&0& & &&\\ 1&-2&1 && &\\ & \ddots &\ddots &\ddots &&\\ && \ddots &\ddots& \ddots &\\ && &1&-2&1 \\ && & &0&0 \end{array}\right] $$

15行目でn×nの単位対角行列Eを作成します。

E = identity(n)

最後の16行目でB行列にcをかけ、Eと足すことで目標の行列が完成しました。

A = E + c*B

by hide

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

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

トラックバックURL: https://teru-hide.com/cfd-python-diagonal-matrix/trackback/
にほんブログ村 英語ブログへ
にほんブログ村