Haskellで数値計算常微分方程式の陰的解法編

以前の記事で陽的Runge-Kutta法を書いたことがあったが、今回は常微分方程式の初期値問題を陰的Runge-Kutta法で解いてみたい。

Implicit Runge-Kutta法
色々な方程式を解いてみよう
- 硬い方程式
- アトラクタなど

Implicit Runge-Kutta法

常微分方程式に初期値があたえられた時には、微分係数を使って端から数値的に積分していけば、(数値誤差の範囲では) 元の関数を得ることができる。この時、それまでのステップの計算値のみを使って次のステップの値を得る陽的な解法と、次のステップの値も使って次のステップの値を得る陰的な解法がある。同じ刻み幅に対しては陰的な解法の方が安定で、特に、硬い方程式に対しては陰的な解法が必要とされる場合があるが、陰的な解法ではステップごとに非線形方程式を解く必要が生じる。

陰的Runge-Kutta法の式と係数行列は以下の通り。

$t_{n+1} = t_n + h$
${\bf x}_{n+1} = {\bf x}_n + h \sum_i b_i {\bf k}_i$
${\bf k}_i = {\bf f}(t_n + c_i h, {\bf x}_n + h \sum_j a_{ij} {\bf k}_j)$

$c_i$	$a_{ij}$
	$b_i$

1段2次

$\frac{1}{2}$	$\frac{1}{2}$
	$1$

2段4次

$\frac{3 - \sqrt{3}}{6}$	$\frac{1}{4}$	$\frac{3 - 2\sqrt{3}}{12}$
$\frac{3 + \sqrt{3}}{6}$	$\frac{3 + 2\sqrt{3}}{12}$	$\frac{1}{4}$
	$\frac{1}{2}$	$\frac{1}{2}$

3段6次

$\frac{5 - \sqrt{15}}{10}$	$\frac{5}{36}$	$\frac{10 - 3\sqrt{15}}{45}$	$\frac{25 - 6\sqrt{15}}{180}$
$\frac{1}{2}$	$\frac{10 + 3\sqrt{15}}{72}$	$\frac{2}{9}$	$\frac{10 - 3\sqrt{15}}{72}$
$\frac{5 + \sqrt{15}}{10}$	$\frac{25 + 6\sqrt{15}}{180}$	$\frac{10 + 3\sqrt{15}}{45}$	$\frac{5}{36}$
	$\frac{5}{18}$	$\frac{4}{9}$	$\frac{5}{18}$

ここではhmatrixの逆行列関数 (LU分解による) を使ってNewton-Raphson法で連立非線形方程式を解くことで、陰的Runge-Kutta法で連立常微分方程式を解いてみる。ほとんどの部分で普通のリストを使っているが、数本の連立方程式で10000ステップ程度を解くなら速度的にも許容できる範囲 (3段6次でも数秒程度) に感じた。速度が必要な場合はvectorパッケージのData.Vector.UnboxedとLLVMのバックエンド (GHCで-fllvmオプション) を使ってみても良いかもしれない。

Numeric Haskell: A Vector Tutorial - HaskellWiki

import Data.List (transpose, unfoldr)
import Data.List.Split (chunksOf)
import qualified Numeric.LinearAlgebra as L

-- 引数をデータ型にまとめる
data Ces = Ces { c :: [Double], b :: [Double], as :: [[Double]], s :: Int }
type Lfs = [[Double] -> Double]
data Ps = Ps { dfs :: Lfs, h :: Double, term :: Double, lim :: Int, err :: Double }

-- 連立非線形方程式のソルバー
delta = 1.0e-9 :: Double -- difference for derivative

derivative :: ([Double] -> Double) -> [Double] -> [Double]
derivative f xs = zipWith dfdx xs [0..]
  where dfdx x i       = (f (replace xs i $ x + delta) - f xs) / delta
        replace xs i x = let (l, _:r) = splitAt i xs in l ++ (x : r)

newtonRaphson :: Int -> Double -> Lfs -> [Double] -> Maybe [Double]
newtonRaphson lim err fs xs
  | lim < 1           = Nothing
  | L.norm_2 fx < err = Just xs
  | otherwise         = newtonRaphson (lim + 1) err fs sx
  where sx = L.toList $ L.vector xs - (L.inv j L.#> fx)
        fx = L.vector $ map ($ xs) fs
        j  = L.fromLists $ map (`derivative` xs) fs

-- Implicit Runge-Kuttaのソルバーの生成器
genRungeKutta :: Ces -> Ps -> [Double] -> Maybe ([Double], [Double])
genRungeKutta (Ces cec ceb ceas s) (Ps dfs h term lim err) st@(t:xs)
  | t > term  = Nothing
  | otherwise = do ks <- newtonRaphson lim err kfs ks0
                   return (st, (t + h):xsf ks ceb)
  where ks0        = concat . transpose $ replicate s xs 
        kfs        = concat $ zipWith3 kfi cec ceas [0..]
        kfi c as i = zipWith kf dfs [0..]
          where kf f n ks = let l  = (t + h * c):xsf ks as
                                k0 = ks !! (i * length dfs + n)
                            in f l - k0
        xsf ks cs  = zipWith csx xs $ transpose $ chunksOf (length xs) ks
          where csx x rks = x + h * sum (zipWith (*) cs rks)

-- 係数行列
ce36 = Ces {
  c  = [(5.0 - sqrt 15.0)/10.0, 0.5, (5.0 + sqrt 15.0)/10.0],
  b  = [5.0/18.0, 4.0/9.0, 5.0/18.0],
  as = [[5.0/36.0, (10.0 - 3.0*sqrt 15.0)/45.0, (25.0 - 6.0*sqrt 15.0)/180.0],
        [(10.0 + 3.0*sqrt 15.0)/72.0, 2.0/9.0, (10.0 - 3.0*sqrt 15.0)/72.0],
        [(25.0 + 6.0*sqrt 15.0)/180.0, (10.0 + 3.0*sqrt 15.0)/45.0, 5.0/36.0]],
  s  = 3 }

ce24 = Ces {
  c  = [(3.0 - sqrt 3.0)/6.0, (3.0 + sqrt 3.0)/6.0],
  b  = [1.0/2.0, 1.0/2.0],
  as = [[1.0/4.0, (3.0 - 2.0*sqrt 3)/12.0 ],
        [(3.0 + 2.0*sqrt 3.0)/12.0, 1.0/4.0]],
  s  = 2 }

ce12 = Ces {
  c  = [1.0/2.0],
  b  = [1.0],
  as = [[1.0/2.0]],
  s  = 1 }

-- ソルバーの生成
rk36 = genRungeKutta ce36
rk24 = genRungeKutta ce24
rk12 = genRungeKutta ce12

1段2次、2段4次、3段6次の係数を用意した。Newton-Raphsonが収束しない可能性を考えてMaybeにしたので、Runge-Kuttaもunfoldrでリスト生成するようにした。最近のunfoldrはFusionが効率的に展開してくれるらしい。hmatrixを使う部分でVectorとリストの変換が入っていたり、陰的Runge-Kuttaの中間係数の計算部分で任意個の方程式を連立する必要があったりして、少しごちゃごちゃしてしまった気がする。termを関数にして、もっと細かい終了条件をあたえてもいいかもしれない。

performance - Efficiency of unfoldr versus zipWith - Stack Overflow

色々な方程式を解いてみよう

可視化にはScheme系のLisp方言のRacketを使ってみた。

https://racket-lang.org

-- 出力のための関数
toString :: Show a => [a] -> String
toString = (++ "\n") . concatMap ((++ " ") . show)

writeLists :: Show a => String -> [[a]] -> IO ()
writeLists file = appendFile file . concatMap toString

-- 方程式を解く例
l_dxdt [t,x,y,z] = -10.0*x + 10.0*y
l_dydt [t,x,y,z] = -x*z + 28.0*x - y
l_dzdt [t,x,y,z] = x*y - 8.0/3.0 * z

lorentz solver = do
  let file = "lorentz.txt"
      sps  = solver $ Ps [l_dxdt,l_dydt,l_dzdt] 0.01 100.0 100 1.0e-9
  writeFile file "t x y z \n"
  writeLists file $ unfoldr sps [0.0,1.0,1.0,1.0]

main :: IO ()
main = lorentz rk36

硬い方程式

方程式に急激に変化する成分と緩やかに変化する成分が混じっている時、急激な方の変化部分に合わせると刻み幅を小さくしなければならないが、緩やかな方の変化部分を全て計算するためには長い時間がかかってしまう。陽的な解法でも刻み幅を誤差に応じて変化させることで、ある程度対応することもできるが、場合によっては必要な刻み幅が浮動小数表現の精度を超えてしまうこともある。そうした硬い方程式は、陽的な解法で解けない場合でも、陰的な解法を使うと上手く積分できることもある。

微分方程式系の硬さを示す指標として、係数行列の固有値の最大値と最小値の比 (の絶対値) があげられる。

${\bf y'} = A{\bf y}$

これを元に、割と硬そうな方程式系を作って解いてみよう。

$\frac{dx}{dt} = 2x + 3y$
$\frac{dy}{dt} = -14400x - 11900y$

陰的Runge-Kutta (3段6次) では10.0という非常に大きな刻み幅でさえ、振動しているものの積分はでき、この場合では右端では正しい答えに収束して行っている。一般的には許容誤差に応じて刻み幅を変化させた方が良いが、ある程度小さな固定刻み幅を取れば十分な場合もある。

f:id:lqtmirage:20181127013307p:plain:w400

陽的Runge-Kutta (4段4次) では刻み幅を0.001まで小さくしても積分が発散してしまい、0.0001というとても小さな刻み幅 (固定刻みの場合) を取らなければならなかった。

f:id:lqtmirage:20181127013311p:plain:w400

こんな風に、陰的な解法では刻み幅を比較的大きく取ることができ、硬い方程式も解けることがある。ただ、ステップごとの計算量は陽的な解法の方が少ないし、陽的な解法でも刻み幅を変えればある程度対応できることもあるので、一概にどちらが優れているとは言えない。実際には解きたい方程式系の特性に応じて様々な解法を使い分けた方が良い。科学的な検証に用いる場合などは、本当に目的の方程式が解けているのか、方程式や解法の特性を踏まえて何れにせよよく考えなければならない。

アトラクタなど

この辺りは特にImplicitに解く必要はないかもしれないが、とりあえず微分方程式を解いて綺麗な絵を描いてみたい。

スパイラルサドル (相空間の特徴)

f:id:lqtmirage:20181127014245p:plain:w400

フィッツフュー南雲方程式

矩形の入力を入れてみた様子
f:id:lqtmirage:20181127091521p:plain:w400
閾値を超えた固定の入力でリミットサイクルがまわる様子
f:id:lqtmirage:20181127091517p:plain:w400
リミットサイクルに引き込まれていく様子
f:id:lqtmirage:20181130092646p:plain:w400