修士論文 要旨

移流現象を示す偏微分方程式の空間シフト付モデル縮減
数理情報学専攻 48186232
指導教員
トA =

1 はじめに(u0u1···)
unt

松尾 洵
松尾 宇泰 教授が得られていることを前提とし, A の特異値分解 A = U SV ⊤ の U の第モデル縮減とは, 高次元の力学系の問題を低次元で近1 列から第 r 列を切り出して Φ とする. A の特異値似して高速に解く手法であり, 特に偏微分方程式 (PDE)

σ1 , σ2 , ... が後ろにいくにつれ急速に減衰する場合にに対しては POD (Proper Orthogonal Decomposition)ui ≈ ΦΦ⊤ ui の近似精度が高くなる. しかし, 移流現象を用いるのが標準的である. POD は, 解をよく表現すを示す PDE のスナップショットの特異値は減衰が比較る低次元空間の基底を選ぶ手法であるが, 少数基底で表的遅く, 精度が悪くなりがちである.

現しやすい拡散現象などではうまくいくものの, 移流2.2 シフト付 POD成分が支配的な解ではうまく解を近似できる空間が選シフト付 POD は, 移流型の PDE に対応した PODべないという問題がある. そこで Reiss et al. [4] では,である [4]. 概要のみ述べるが, 解が複数（NS 個）の異POD の拡張アルゴリズムとして, 複数の速度の移流成なる速度成分を持つ場合に, 解を分を持つような解を少数の基底の空間シフトの重ね合ui ≈わせで近似する方法を提案した. しかしながら, Reissらの論文では Snapshot の近似をひとまずの目標としており, 実際にこの手法を用いたモデル縮減を行う方法については述べられていない.
本研究では, Reiss らの論文をもとに移流現象を示すPDE に対するモデル縮減を考え, また既存の構造保存モデル縮減と組み合わせる方法も提案し, KdV 方程式()
Ticℓ Φℓ ûℓi(4)ℓ=1と表現して, この Φℓ を求めるアルゴリズムがシフト付 POD である. ただし Tic (u) は波の速度 c と刻み
幅 ∆x, ∆t が所与で s := c∆t∆x ∈ Z のとき, u を s 成分だけ循環シフトさせる作用素と定義し, 各速度 cℓ(ℓ = 1, ..., NS ) は今回は既知とする.
2.3に対する数値実験で検証を行った.

NS∑

そのほかの既存研究移流の構造を用いるのではなく, 離散化し縮減した後2 既存研究も保存量が一定となるようにモデル縮減する手法が提案されている [2, 3]. 今回の数値実験にて比較手法とし

2.1 POD を用いたモデル縮減
時間変化する n 次元実ベクトル u(t) の常微分方程式:

u : [0, T ] → Rn , f : Rn × (0, T ) → Rn ,du(t)= f (u(t), t), u(0) = u0dt(1)について, 初期値 u0 ∈ R , f に対して, u(t) (t ≥ 0)nを求める問題を考える. モデル縮減では, 微分方程式をr(≪ n) 次元の微分方程式に落とすため,てこの手法を採用する.
また, 今回の数値実験で対象とする周期境界条件でのKdV 方程式{ut + 6uux + uxxx = 0, x ∈ R,u(x + L, t) = u(x, t) (∀ x, t)t ∈ (0, T )(5)に対し, 離散変分法をもとにしてエネルギーを保存するよう離散化して更新式をと近似する. Galerkin 射影を用いる方法では, (1) 式のui+1 − ui+ D1 {ui+1 ⊙ ui+1 + ui+1 ⊙ ui∆t()}
ui+1 + ui+ui ⊙ ui + D2=02u(t) の位置に Φû(t) を代入した式の両辺に左から Φ⊤と得ると, 比較的安定した解が得られる [1]. ただしu(t) ≈ ũ(t) = Φû(t)(Φ ∈ Rn×r, û(t) ∈ R )r(2)をかけて得た û(t) の更新式

Φ⊤ Φdû(t)= Φ⊤ f (Φû(t), t)
dt(6)D1 , D2 はそれぞれ 1 階, 2 階の中心差分を表す行列で(3)を解いて近似解 ũ(t) = Φû(t) を得る.

POD は縮減行列 Φ の有名な構成方法である. 各時刻の解 ui ≈ u(i∆t) を並べた行列であるスナップショッあり, ⊙ は成分ごとの積である.

3 提案手法
シフト付 POD を用いた移流型の PDE に対するモデル縮減手法を 4 つ提案する.

移流現象を示す偏微分方程式の空間シフト付モデル縮減 (松尾 洵)

提案手法では, シフト行列 K を用いて近似解を

0NS∑ cℓ i∆t 1(K ∆x Φℓ )ûℓi , K = ui ≈ℓ=1···..

0.
1

10 .. . 0(7)と表す. Kℓ := KΦ∗i :=((K1 )i Φ1cℓ ∆t∆xS ⊤ ⊤, û∗i := ((û1i )⊤ ,)..., (ûNi ) ) ,として ui ≈· · · (KNS )i ΦNSΦ∗i û∗i とも表せる. Kℓ は波の速度 cℓ に対し +cℓ t だけ空間シフトさせる行列であり, Φ∗i は時間変化する縮減行列とみなせる.

3.1 提案手法 1, 提案手法 4

図 1. 計算時間と相対誤差の比較 (2-ソリトン). 提案

1:赤× , 提案 2:青＋, 提案 4:緑＊, 比較手法:マゼンタ○.

提案手法 1 の縮減後の更新式は, 縮減前の更新式に対し ui , ui+1 の位置にそれぞれ Φ∗i û∗i , Φ∗i+1 û∗i+1 を代入し, 両辺に左から Pi⊤ ∈ Rr×n を掛けることで得る.
ここで, 速度成分が単一（NS = 1）のとき Pi⊤ として(Φ∗i+1 )⊤ を選ぶことが, 非線形項の近似に関する適当な仮定の下で最適であることを本文中で示した. 速度成分複数（NS ≥ 2）のときは Pi⊤ = (Φ∗i + Φ∗i+1 )⊤ とする.
提案手法 4 は, 提案手法 1 の Pi⊤ を掛ける前の式に,とした. Miyatake の手法のノルム保存手法に対応し,NS = 1 のときはノルム保存が成り立つ.

4 数値実験
周期境界条件の KdV 方程式で 1-ソリトン解, 2-ソリトン解に対しモデル縮減で時間発展を追った. 比較手法として Gong et al. [2], Miyatake [3] の手法を採用す

両辺の差の L2 ノルムを最小にする行列を掛ける.

る. 例えば, 2-ソリトンの提案手法 1, 2, 4 と, 比較手法

3.2 提案手法 2で計算時間と相対誤差を比較したものは図 1 である.

提案手法 1 に加え, 非線形項の計算方法の工夫 (ten-sor manipulation [5]) と, 計算時間が n に依存する部分を事前計算をすることで, 特に r が小さいときの高速化

5 結論
シフト付 POD で求まる縮減行列を用いたモデル縮を図ったものを提案手法 2 とする.

減手法を提案した. 比較手法と同程度・それ以上のパ

3.3 提案手法 3

フォーマンスを示したものの, 速度成分が複数ある場合

提案手法 3 は, Miyatake の手法 [3] と組み合わせの更なる精度向上と高速化が今後の主な課題である.

た手法である. 簡単のため NS = 1 のときで記すと,U = diag(u) として空間離散化した KdV 方程式

参考文献[1] 降籏大介: 偏微分方程式に対する差分スキームの離散的変分による統一的導出の研究. 博士論文, 東京大学 (1997).
(8)[2] Y. Gong, Q. Wang and Z. Wang: Structure-preserving
Galerkin POD reduced-order modeling of Hamilto（ た だ し S(u) := (U (−2D1 ) + (−2D1 )U − D1 D2 ),nian systems. Computer Methods in Applied Mechan∇u H(u) = u とおく）に対し ui ≈ Φ∗i û∗i の近似を
ics and Engineering 315 (2017): 780–798.
[3]Y.
Miyatake: Structure-preserving model reduction
念頭に式変形及び時間離散化しfor dynamical systems with a first integral. Japan)( 1Journal of Industrial and Applied Mathematics 36.3
û1i+1 − û1i
ûi+1 + û1i
11(2019): 1021–1037.
= Sr (ûi+1 , ûi )∇û H̃
(9)∆t2[4] J. Reiss, P. Schulze, J. Sesterhenn and V. Mehrmann:
The shifted proper orthogonal decomposition: Aを得る. ただし ∆i > 0 を十分小さいパラメータとしてmode decomposition for multiple transport phenom( () ena. SIAM Journal on Scientific Computing 40.3
11()(K1 ) 2 Φ1 û1i+1 + (K1 )− 2 Φ1 û1i
(2018): A1322–A1344.
Sr û1i+1 , û1i := Φ⊤S
12[5] Z. Wang, I. Akhtar, J. Borggaard and T. Iliescu:)
Proper orthogonal decomposition closure models for1 1−∆i∆i(K1 − K1 ) Φ1 turbulent flows: a numerical comparison. Computer−∆t 2∆i
Methods in Applied Mechanics and Engineering 237(10)(2012): 10–26.

du = S(u)∇u H(u)dt