🌰
Cyclic Reductionでポアソン方程式を解く【コードあり】

2025/09/26に公開
Cyclic Reductionの"くりっくり"の部分

 アブストラクトポアソン方程式を直交格子上で離散化すると、三重対角行列を係数行列とする線形連立方程式になる。
三重対角行列を係数行列とする線形連立方程式の解法としてCyclic Reduction法がある。
今回は一次元のポアソン方程式を対象とする。この場合Cyclic Reduction法は不等間隔格子にも対応できる。

pythonコードの全体は最後にあります。

 ポアソン方程式を不等間隔格子上で離散化する今回解きたいのは一次元のポアソン方程式

\frac{d^2p}{d x^2}=S
です。p=p(x)、S=S(x)です。pは未知、Sは既知。境界条件は後述。


これを図のような一次元の格子上で離散化します。変数は格子中央で定義されます。格子幅は等間隔でなくともよいです。格子幅を\Delta x_1(i)、格子中央と格子中央の距離を\Delta x_2(i)とします。また、\Delta x_1(0)=\Delta x_1(1)、\Delta x_1(im+1)=\Delta x_1(im)としておきます。

\frac{(p_{i+1}-p_i)/\Delta x_2(i)-(p_i-p_{i-1})/\Delta x_2(i-1)}{\Delta x_1(i)}=S_i
i(=0からim+1)は格子番号。p_0とp_{im+1}は領域外です。等間隔格子であればよくみる形(p_{i+1}-2p_i+p_{i-1})/\Delta x^2=S_iに一致します。

左辺を式変形すると、

\frac{1}{\Delta x_2(i-1)\Delta x_1(i)}p_{i-1}
-\left(\frac{1}{\Delta x_2(i-1)\Delta x_1(i)}+\frac{1}{\Delta x_2(i)\Delta x_1(i)}\right)p_i
+\frac{1}{\Delta x_2(i)\Delta x_1(i)}p_{i+1}
=S_i
となります。

 境界条件今回は境界でp=0とします。

するとp_0=-p_1、p_{im+1}=-p_{im}なので、i=1番目の式は、

-\left(\frac{2}{\Delta x_2(0)\Delta x_1(1)}+\frac{1}{\Delta x_2(1)\Delta x_1(1)}\right)p_1
+\frac{1}{\Delta x_2(1)\Delta x_1(1)}p_2
=S_1
となります。i=im番目も同様。

 （補足）たとえば境界上でp=1のときはp_0=2-p_1となり、p_1の係数に吸収させられません。このような場合、p_1を含まない部分は右辺にまわしてしまって、合わせてあらたにS_1とすればよいと思います。

 Cyclic Reduction法の準備　～三重対角行列～書き下すと、こうなります↓↓

\begin{align}
-\left(\frac{2}{\Delta x_2(0)\Delta x_1(1)}+\frac{1}{\Delta x_2(1)\Delta x_1(1)}\right)p_1
+\frac{1}{\Delta x_2(1)\Delta x_1(1)}p_2
&=S_1\\
\frac{1}{\Delta x_2(1)\Delta x_1(2)}p_{1}
-\left(\frac{1}{\Delta x_2(1)\Delta x_1(2)}+\frac{1}{\Delta x_2(2)\Delta x_1(2)}\right)p_2
+\frac{1}{\Delta x_2(2)\Delta x_1(2)}p_{3}
&=S_2\\
\frac{1}{\Delta x_2(2)\Delta x_1(3)}p_{2}
-\left(\frac{1}{\Delta x_2(2)\Delta x_1(3)}+\frac{1}{\Delta x_2(3)\Delta x_1(3)}\right)p_3
+\frac{1}{\Delta x_2(3)\Delta x_1(3)}p_{4}
&=S_3\\
\vdots \\
&=S_{im}
\end{align}
行列形式で書くと係数行列が三重対角行列になります。
後のために式変形をしておきます。

i番目の式に対して、p_iの係数で割ります。表記が煩雑になるのでp_{i-1}の係数はL_i、p_{i+1}の係数はR_iとおいておきます。右辺S_iは割ったあとの値を改めてS_iと書くとします。

そうすると上式は、行列の形Ap=Sでかくと係数行列Aは次のようになります。

A=
\begin{pmatrix} 
  1   & R_1     &     &        & \\
  L_2 & 1       & R_2 &        &  \huge{0} \\
      & L_3     & 1   & R_3    &  \\
      &         & L_4 & 1      & R_4   &  \\
      & \huge{0}&     & \ddots &\ddots &\ddots\\
      &         &     &        & L_{im} & 1
\end{pmatrix} 
R= np.zeros(im+2)
L= np.zeros(im+2)

R[1:im  ]=(1./dx2[1:im]/dx1[1:im  ])/(-1./dx2[0:im-1]/dx1[1:im  ]-1./dx2[1:im  ]/dx1[1:im  ])
L[2:im+1]=(1./dx2[1:im]/dx1[2:im+1])/(-1./dx2[1:im  ]/dx1[2:im+1]-1./dx2[2:im+1]/dx1[2:im+1])
S[2:im  ]=S[2:im  ]/(-1./dx2[1:im-1]/dx1[2:im  ]-1./dx2[2:im  ]/dx1[2:im  ])

R[1] =1./dx2[1   ]/dx1[1] /(-2./dx2[0   ]/dx1[1 ]-1./dx2[1 ]/dx1[1 ])
L[im]=1./dx2[im-1]/dx1[im]/(-1./dx2[im-1]/dx1[im]-2./dx2[im]/dx1[im])
S[1] =S[1] /(-2./dx2[0   ]/dx1[1 ]-1./dx2[1 ]/dx1[1 ])
S[im]=S[im]/(-1./dx2[im-1]/dx1[im]-2./dx2[im]/dx1[im])

 Cyclic Reduction法（前半）Cyclic Reduction法（以下CR法）は次のようなアルゴリズムです。（この記事↓↓を参考にしています）
https://www.bing.com/ck/a?!&&p=e4eb087fa22e845113e95fd6df3194d8a8c25f18f65e3a5897666b2d6448e7b9JmltdHM9MTc1ODg0NDgwMA&ptn=3&ver=2&hsh=4&fclid=09904b1e-77ef-66f7-0f82-4562767e670d&psq=cyclic+reduction&u=a1aHR0cHM6Ly9jYXRhbG9nLmxpYi5reXVzaHUtdS5hYy5qcC9vcGFjX2Rvd25sb2FkX21kLzE3MjgyL3AzNTMucGRm
まず！CR法は格子数がim=2^N-1のときに使える手法です。（im=2^Nのときも少しの変更で適応できるが今回はなし。）そのうえで、

\begin{align}
L_{i-1}p_{i-2}+p_{i-1}+R_{i-1}&p_{i}&=S_{i-1}\\
L_{i}p_{i-1}+&p_{i}+R_{i}p_{i+1}&=S_{i}\\
L_{i+1}&p_{i}+p_{i+1}+R_{i+1}p_{i+2}&=&S_{i+1}
\end{align}
i番目の式から「i-1番目の式をL_i倍したもの」と「i+1番目の式をR_i倍したもの」を引きますと、

\begin{equation}
-L_{i}L_{i-1}p_{i-2}+(1-L_{i}R_{i-1}-R_{i}L_{i+1})p_i-R_iR_{i+1}p_{i+2}=S_i-L_iS_{i-1}-R_iS_{i+1}
\end{equation}
となり、p_{i-2}、p_i、p_{i+2}の式となります。前節と同じようにp_iの係数で両辺を割れば、

\begin{equation}
L_{i}^{(1)}p_{i-2}+p_i+R_i^{(1)}p_{i+2}=S_i^{(1)}
\end{equation}
となり、元の形にできます。係数と右辺は新たに置きなおしています。

すべての偶数番目のiに対してこの操作をすれば式の数は半分になります（正確には2^N-1個が2^{N-1}-1個になる）。したがって、それをN-1回くり返せば、最終的にp_{2^{N-1}}に関する一本の式

\begin{equation}
p_{2^{N-1}}=S_{2^{N-1}}^{(N-1)}
\tag{2}
\end{equation}
が得られます。

前半おわり。
N=int(np.log(float(im+1))/np.log(2.)) 

for n in range(1,N-1+1):
    dm=2**(n-1)
    for m in range(2**n , im+1-2**n+1 , 2**n):
        dc=L[m]*R[m-dm]-1.+R[m]*L[m+dm]
        S[m]=(L[m]*S[m-dm]-S[m]+R[m]*S[m+dm])/dc
        R[m]=(                  R[m]*R[m+dm])/dc
        L[m]=(L[m]*L[m-dm]                  )/dc

 Cyclic Reduction法（後半）得られたp_{2^{N-1}}から、上の手順を逆にたどってすべてのp_iを求めていきます。
具体的には、(2)のひとつ前の段階でp_{2^{N-1}-2^{N-2}}（とp_{2^{N-1}}）の式、p_{2^{N-1}+2^{N-2}}（とp_{2^{N-1}}）の式があるはずなのでここからp_{2^{N-1}-2^{N-2}}とp_{2^{N-1}+2^{N-2}}が得られます。すなわち、（添え字だけ見ると）2^{N-1}の値から\pm 2^{N-2}番目の値が求められます。

次の段階では、ここまで得られた3つの値から、それぞれの\pm 2^{N-3}番目の値が求められます。これを繰り返せば（N-1回ですかね）、すべてのp_iが計算できます。
for n in range(N-1,0,-1):
    dm=2**(n-1)
    for m in range(2**n , im+1-2**n+1 , 2**(n+1)):
        S[m-dm]=S[m-dm]-L[m-dm]*S[m-2*dm]-R[m-dm]*S[m]
        S[m+dm]=S[m+dm]-L[m+dm]*S[m]     -R[m+dm]*S[m+2*dm]
p[1:im+1]=S[1:im+1]
# 境界条件
p[0]=-p[1] 
p[im+1]=-p[im]

 検証テスト問題として、なんでもいいですが、簡単に、

\begin{equation}
\frac{d^2p}{dx^2}=\cos x +\sin x 
\end{equation}
としてやってみます。0 \leq x \leq 2\piで、境界条件をp(0)=p(2\pi)=0とすると解は

\begin{equation}
p(x)=-\cos x -\sin x +1
\end{equation}
となるはず。

 計算条件格子数im=2^5-1=31（領域外含めると33）
領域サイズL=2\pi
一旦、等間隔格子\Delta x_1=\Delta x_2=L/im

 結果解析解と併せて描いています↓↓



うまくできていそうです。

 検証：不等間隔格子の場合CR法のイイトコロは格子が等間隔でなくてもよいことと、処理がそこそこはやい（らしい↓↓）ことです。
https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&ved=2ahUKEwjZs7iMsvaPAxVBh1YBHdZWGg8QFnoECB8QAQ&url=https%3A%2F%2Fdl.ndl.go.jp%2Fview%2FprepareDownload%3FitemId%3Dinfo%253Andljp%252Fpid%252F10604699%26contentNo%3D1&usg=AOvVaw2fjWb5smiy2tp-X-jRdmOx&opi=89978449
高速なポアソン方程式解法として名高い、高速フーリエ変換やコサイン変換を使うものがあります↓↓
https://qiita.com/eigs/items/cb607d647bc20c7db809
が、これは等間隔格子でのみ適用可能です。その点CR法は不等間隔でもいいというのがいい。計算速度については今回は調べませんが、SOR法とかよりはさすがにはやそうな気がします。
上のテスト問題を格子幅を変えてやってみます。

 ケース1等間隔（1:1:1\cdots）だったのを0.5:0.5:2\cdotsにします。

\begin{equation}
  \Delta x_1(i)=
  \begin{cases}
    2 L/im 　　& \text{if  } i \text{は3の倍数,} \\
    0.5 L/im   & \text{if  } i \text{は3の倍数でない.} 
  \end{cases}
\end{equation}
境界近くだけ\Delta x_1=L/imにしておきます。その他の条件は上と同じ。
# 不等間隔格子（ケース1）------------------------------
dx1[:]=Lx/float(im)
dx1[3:im-1:3]=2.*Lx/float(im)
dx1[4:im:3]=0.5*Lx/float(im)
dx1[5:im:3]=0.5*Lx/float(im)

 結果

できていそうです。少し精度悪いか？格子数を増やせば大丈夫...。間違いがあれば教えてください。

 ケース2領域の端で大きめの格子で、中央に近づくにつれて線形に細かくなっていくようにします。
端　：\Delta x_1(1)=\Delta x_1(im)=1.7L/im
中央：\Delta x_1((im+1)/2)\risingdotseq 0.25L/im
# 不等間隔格子（ケース2）------------------------------
dmax=1.7*Lx/float(im) # L/im < dmax < 2L/(im+1) 
dmin=2./float(im-1)*(Lx-float(im+1)/2.*dmax)
dx1[1]=dmax
dx1[int((im+1)/2)]=dmin
for i in range(2,int((im+1)/2)):
    dx1[i]=(dmin-dmax)/float((im+1)/2-1)*float(i-1)+dmax
dx1[int((im+1)/2+1):im+1:1]=dx1[int((im+1)/2-1):0:-1]
dmaxを自分で決めればあとは勝手に決まるようにしています。L/im < dmax < 2L/(im+1)の範囲でないと、dminがdmaxより大きくなったりdminが負の値になったりするので注意。

 結果

できていそうです。

 おしまい（Unfortunately）うまく実装できていそうでよかったです。短いコードで、特別なライブラリも使わず実装できるかっこいいアルゴリズムだと思います。
残念ながら二次元とか三次元のポアソン方程式をCR法で解こうとなると、難しいようです。二次元の場合は調べるといくつか出てきますが、等間隔格子の場合ばかりです。等間隔ならFFTでいいからなあ。本記事で言うところのR_iとかL_iが二次元の場合は行列になるから逆行列とかを計算しないといけなくて難しい、とか？どなたか教えてください。
なので、x方向に等間隔、y方向に不等間隔みたいな格子のときに、x方向はFFT、y方向はCR法で解く、みたいなことをしたい。そのうちやります。

 コード全体import numpy as np
import matplotlib.pyplot as plt

im=32-1
Lx=2.*np.pi
dx1=np.zeros(im+2)
dx2=np.zeros(im+2)
x  =np.zeros(im+2)
p  =np.zeros(im+2)

# 等間隔格子----------------------------
dx1[:]=Lx/float(im)
# -------------------------------------

# # 不等間隔格子（ケース1）-------------------
# dx1[:]=Lx/float(im)
# dx1[3:im-1:3]=2.*Lx/float(im)
# dx1[4:im:3]=0.5*Lx/float(im)
# dx1[5:im:3]=0.5*Lx/float(im)
# # ----------------------------------------

# # 不等間隔格子（ケース2）------------------------------
# dmax=1.7*Lx/float(im) # L/im < dmax < 2L/(im+1) 
# dmin=2./float(im-1)*(Lx-float(im+1)/2.*dmax)
# dx1[1]=dmax
# dx1[int((im+1)/2)]=dmin
# for i in range(2,int((im+1)/2)):
#     dx1[i]=(dmin-dmax)/float((im+1)/2-1)*float(i-1)+dmax
# dx1[int((im+1)/2+1):im+1:1]=dx1[int((im+1)/2-1):0:-1]
# #--------------------------------------------------------

dx1[0]=dx1[1]
dx1[im+1]=dx1[im]

x[0]=-0.5*dx1[0]
for i in range(1,im+2):
    x[i]=x[i-1]+0.5*(dx1[i-1]+dx1[i])
dx2[0:im+1]=x[1:im+2]-x[0:im+1]

S= np.zeros(im+2)
for i in range(0,im+2):
    S[i]=np.cos(x[i])+np.sin(x[i])
S[0]=0.
S[im+1]=0.

R= np.zeros(im+2)
L= np.zeros(im+2)

R[1:im  ]=(1./dx2[1:im]/dx1[1:im  ])/(-1./dx2[0:im-1]/dx1[1:im  ]-1./dx2[1:im  ]/dx1[1:im  ])
L[2:im+1]=(1./dx2[1:im]/dx1[2:im+1])/(-1./dx2[1:im  ]/dx1[2:im+1]-1./dx2[2:im+1]/dx1[2:im+1])
S[2:im  ]=S[2:im  ]/(-1./dx2[1:im-1]/dx1[2:im  ]-1./dx2[2:im  ]/dx1[2:im  ])

R[1] =1./dx2[1   ]/dx1[1] /(-2./dx2[0   ]/dx1[1 ]-1./dx2[1 ]/dx1[1 ])
L[im]=1./dx2[im-1]/dx1[im]/(-1./dx2[im-1]/dx1[im]-2./dx2[im]/dx1[im])
S[1] =S[1] /(-2./dx2[0   ]/dx1[1 ]-1./dx2[1 ]/dx1[1 ])
S[im]=S[im]/(-1./dx2[im-1]/dx1[im]-2./dx2[im]/dx1[im])

N=int(np.log(float(im+1))/np.log(2.)) 

for n in range(1,N-1+1):
    dm=2**(n-1)
    for m in range(2**n , im+1-2**n+1 , 2**n):
        dc=L[m]*R[m-dm]-1.+R[m]*L[m+dm]
        S[m]=(L[m]*S[m-dm]-S[m]+R[m]*S[m+dm])/dc
        R[m]=(                  R[m]*R[m+dm])/dc
        L[m]=(L[m]*L[m-dm]                  )/dc
for n in range(N-1,0,-1):
    dm=2**(n-1)
    for m in range(2**n , im+1-2**n+1 , 2**(n+1)):
        S[m-dm]=S[m-dm]-L[m-dm]*S[m-2*dm]-R[m-dm]*S[m]
        S[m+dm]=S[m+dm]-L[m+dm]*S[m]     -R[m+dm]*S[m+2*dm]
p[1:im+1]=S[1:im+1]
# 境界条件
p[0]=-p[1] 
p[im+1]=-p[im]

fig1=plt.figure(figsize=(12,6))
ax1= fig1.add_subplot()
ax1.plot(x,p,'.',label='CR sol.')
x=np.linspace(0,Lx,10000)
print(x)
ax1.plot(x,-np.sin(x)-np.cos(x)+1.,label='analytical sol.')
ax1.legend()
# plt.savefig('CR.png')
plt.show()
アブストラクト

ポアソン方程式を不等間隔格子上で離散化する

境界条件

（補足）

Cyclic Reduction法の準備 ～三重対角行列～

Cyclic Reduction法（前半）

Cyclic Reduction法（後半）

検証

計算条件

結果

検証：不等間隔格子の場合

ケース1

結果

ケース2

結果

おしまい（Unfortunately）

コード全体

Discussion

Cyclic Reduction法の準備　～三重対角行列～
