狠狠撸

はじめにタイムステッピングについて数値例まとめ
チェビシェフ級数を用いた
タイムステッピングによる
常微分方程式系の精度保証付き数値解法
舩越康太 (筑波大学 M1), 高安亮紀 (筑波大学)
日本応用数理学会 2019 年度年会
2019/9/4
1 / 29

目次
1 はじめに
問題設定
既存手法の紹介
提案手法の紹介
2 タイムステッピングについて
概要
端点評価
3 数値例
問題設定
計算結果
4 まとめ
本手法の特徴
今後の課題
2 / 29

問題設定
常微分方程式系の初期値問題
dx
dt
= f(t, x), x(0) = x0, t ∈ (0, h). (1)
各文字の説明
x(t), x0 ∈ Rn
, h > 0,
f(t, x) =
?
?
?
f1(t, x)
...
fn(t, x)
?
?
? ,
fi(t, x) : R × Rn
→ R, i = 1, . . . , n.
3 / 29

既存手法の紹介
テイラー展開ベース
代表例:
Lohner 法 [Lohner 1987] (AWA)
柏木らによる方法 [Kashiwagi 1994] (kv)
テイラーモデル [Berz-Makino 1998] (COSY)
C1 Lohner 法 [Zgliczynski 2002] (CAPD)
スペクトル法ベース
代表例:
占部によるガレルキン法 [Urabe 1965]
大石による方法 [Oishi 1995]
Radii polynomial approach [Lessard 2014]
4 / 29

1. 近似解の生成
チェビシェフ多項式による展開の形で近似解を生成する.
5 / 29

↓
2. 解の精度保証
簡易ニュートン法とバナッハの不動点定理を用いて, 1 つの時間ス
テップで解の精度保証を行う.
5 / 29

↓
2. 解の精度保証
簡易ニュートン法とバナッハの不動点定理を用いて, 1 つの時間ス
テップで解の精度保証を行う.
↓
3. タイムステッピング
複数の時間ステップに渡って解を延長する.
5 / 29

近似解の生成
Chebfun 1 による (1) 式の近似解
?xi(t) =
M?1∑
l=0
?xl,iTl(t).
各文字の説明
M : チェビシェフ多項式の数,
Tl(t) : l 次第 1 種チェビシェフ多項式.
1
オックスフォード大学の数値解析グループによって開発されている
MATLAB 上で動く数値計算ソフトウェア (http://www.chebfun.org/)
6 / 29

バナッハ空間 X とその閉部分集合 Wγ
バナッハ空間 X
時間区間を J
def
= (0, h) として
X = C(J; Rn
)
def
= { x ∈ Rn
| xi(t) ∈ C(J) }.
?x の近傍 Wγ(?x)
Wγ(?x)
def
= { x ∈ X | ∥xi ? ?xi∥ ≤ γi, γi > 0, x(0) = x0 },
∥x∥
def
= sup
t∈J
|x(t)|.
7 / 29

作用素 F と簡易ニュートン写像 T
作用素 F
(F(x))(t)
def
=
dx
dt
? f(t, x).
簡易ニュートン写像 T
(T (x))(t)
def
= x(t) ? A(F(x))(t).
8 / 29

作用素 A?
作用素 A?
A? def
= DF(?x) =
d
dt
? Df(?x).
f の ?x におけるヤコビ行列 Df(?x)
Df(?x)
def
=
?
?
?
?f1
?x1
(?x) . . . ?f1
?xn
(?x)
...
...
?fn
?x1
(?x) . . . ?fn
?xn
(?x)
?
?
? .
9 / 29

作用素 A
作用素 A
Ag
def
= Φ(t, s)c0 +
∫ t
s
Φ(t, r)g(r) dr.
ただし, Φ(t, s) ≡ Φ(t) ∈ Rn×n は以下の線形化方程式を満たす.
(1) 式の線形化方程式
dΦ
dt
= Df(?x)Φ, Φ(s) = I.
この Φ(t) は精度保証付きで得ることができる. また, AA? = I.
10 / 29

主定理
Theorem 1.1 (Newton-Kantorovich type argument)
いま ?x ∈ X, i = 1, . . . , n に対して
|[x(0) ? ?x(0)]i| ≤ εi,
∥[AF(?x)]i∥X ≤ Yi,
∥[A(DF(b) ? DF(?x))ζ]i∥L(X) ≤ Zε
i (γ), ? b, ζ ∈ Wγ(?x)
が成り立つとする. このとき
Pi(γ)
def
= Zε
i (γ) ? γi + Yi
に対し, P(?γ) ≤ 0 となる ?γ > 0 が存在すれば, F(x?) = 0 を満た
す解 x? が W?γ(?x) 内に ( 局所 ) 一意存在する.
11 / 29

タイムステッピング
t
x
O
?x(t)
t0 t1
1 つの時間ステップでの近似解を得る.

t
x
O
?x(t)
t0 t1
W?γ(?x)
定理 1.1 (J = (t0, t1)) により, 解の包含を得る.

t
x
O
?x(t)
t0 t1
W?γ(?x)
t1
} ε
} ε
t = t1 における誤差評価を行う.
12 / 29

t
x
O
?x(t)
t0 t1
W?γ(?x)
t2
次の時間ステップでの近似解を得る.

t
x
O
?x(t)
t0 t1
W?γ(?x)
t2
定理 1.1 (J = (t1, t2)) により, 再度解の包含を得て, これを繰り返す.
13 / 29

端点評価
端点評価の方法
誤差を ?(t)
def
= x(t) ? ?x(t) (t ∈ (0, h)) とすると
?(t) = Φ(t, 0)?(0) +
∫ t
0
Φ(t, s)g(s) ds
と表され, 端点における絶対誤差を考えると
|?(h)| ≤ |Φ(h, 0)|ε + h|Φ(h, 0)| ∥Φ(0, s)∥ ∥g(s)∥ = ε
(∥Φ(0, s)∥ = sup
s∈(0,h)
|Φ(0, s)|)
と評価される.
14 / 29

端点評価 (g(s) について)
各文字の説明
g(s) =
∫ 1
0
D2
f((1 ? θ)?x + θb) dθ(b ? ?x)ζ ?
(
d?x
ds
? f(s, ?x)
)
(θ ∈ (0, 1), b ∈ Wγ(?x), ζ ∈ Sε
γ),
ここで
Sε
γ
def
= { s ∈ X | ∥si∥ ≤ γi, |si(0)| ≤ εi }.
15 / 29

?(t)
dt
=
dx
dt
?
d?x
dt
= f(t, x) ? f(t, ?x) ?
(
d?x
dt
? f(t, ?x)
)
= Df(?x)?(t) + Df(b)?(t) ? Df(?x)?(t) ?
(
d?x
dt
? f(t, ?x)
)
= Df(?x)?(t) +
∫ 1
0
D2
f((1 ? θ)?x + θb) dθ(b ? ?x)ζ ?
(
d?x
dt
? f(t, ?x)
)
= Df(?x)?(t) + g(s)
↓
?(t) = Φ(t, 0)?(0) +
∫ t
0
Φ(t, s)g(s) ds
16 / 29

端点評価 (Φ(t, 0), Φ(0, s) について)
Φ(t, 0), Φ(0, s) の計算方法
Φ(t, 0) ≡ Φ(t), Φ(0, s) ≡ Ψ(s) として
dΦ
dt
= Df(?x)Φ(t), Φ(0) = I,
dΨ
ds
= ?Ψ(s)Df(?x), Ψ(0) = I
のぞれぞれを数値計算により解いた後, Radii polynomial approach
[Lessard 2014] によって精度保証を行う.
17 / 29

Ψ(s)Φ(s) = I
??
dΨ(s)
ds
Φ(s) + Ψ(s)
dΦ(s)
ds
= 0
??
dΨ(s)
ds
Φ(s) + Ψ(s)Df(?x)Φ(s) = 0
??
dΨ(s)
ds
Φ(s) = ?Ψ(s)Df(?x)Φ(s)
??
dΨ(s)
ds
= ?Ψ(s)Df(?x)
18 / 29

問題設定
Lorenz 方程式
?
??
??
d
dt x1(t) = 10(x2(t) ? x1(t)),
d
dt x2(t) = (28 ? x3(t))x1(t) ? x2(t),
d
dt x3(t) = x1(t)x2(t) ? 8
3x3(t).
19 / 29

問題設定
Lorenz 方程式
?
??
??
d
dt x1(t) = 10(x2(t) ? x1(t)),
d
dt x2(t) = (28 ? x3(t))x1(t) ? x2(t),
d
dt x3(t) = x1(t)x2(t) ? 8
3x3(t).
↓
Lorenz 方程式
dx
dt
= f(t, x) , x(0) =
?
?
1
1
1
?
? , t ∈ (0, h).
19 / 29

計算結果 (各 h に対する ?γ と実行時間)
計算環境
macOS 10.13.4
CPU: 2.7GHz Intel Core i5
MATLAB 2016b
INTLAB version 11
Chebfun version 5.7.0
20 / 29

計算環境
macOS 10.13.4
MATLAB 2016b
INTLAB version 11
h max(maxi ?γi) 計算できたステップ数到達時刻実行時間
0.0625 3.6384 × 10?1 61 step 3.8125 sec. 136.03 sec.
0.125 9.3574 × 10?2 52 step 6.5 sec. 124.47 sec.
0.25 4.0847 × 10?3 28 step 7 sec. 66.21 sec.
20 / 29

計算環境
macOS 10.13.4
MATLAB 2016b
INTLAB version 11
0.0625 3.6384 × 10?1 61 step 3.8125 sec. 136.03 sec.
0.125 9.3574 × 10?2 52 step 6.5 sec. 124.47 sec.
0.25 4.0847 × 10?3 28 step 7 sec. 66.21 sec.
0.32 1.6273 × 10?3 16 step 5.12 sec. 40.83 sec.
21 / 29

計算結果 (誤差半径の時間変化)
0 1 2 3 4 5 6 7
t
10 -12
10 -10
10 -8
10 -6
10 -4
10 -2
10 0
h=0.0625
h=0.125
h=0.25
h=0.32
Figure: 誤差半径 γ の時間変化
22 / 29

計算結果 (誤差の時間変化)
0 1 2 3 4 5 6 7
t
10 -16
10 -14
10 -12
10 -10
10 -8
10 -6
10 -4
10 -2
10 0
h=0.0625
h=0.125
h=0.25
h=0.32
Figure: 誤差 ? の時間変化
23 / 29

計算環境
macOS 10.13.4
MATLAB 2016b
INTLAB version 11
0.0625 3.6384 × 10?1 61 step 3.8125 sec. 136.03 sec.
0.125 9.3574 × 10?2 52 step 6.5 sec. 124.47 sec.
0.25 4.0847 × 10?3 28 step 7 sec. 66.21 sec.
0.32 1.6273 × 10?3 16 step 5.12 sec. 40.83 sec.
24 / 29

計算環境
macOS 10.13.4
MATLAB 2016b
INTLAB version 11
0.0625 3.6384 × 10?1 61 step 3.8125 sec. 136.03 sec.
0.125 9.3574 × 10?2 52 step 6.5 sec. 124.47 sec.
0.25 4.0847 × 10?3 28 step 7 sec. 66.21 sec.
0.32 1.6273 × 10?3 16 step 5.12 sec. 40.83 sec.
25 / 29

最も良かった計算結果のプロット
0
30
10
20
20
20
x3
30
10 10
x2
40
x1
0 0
50
-10 -10
-20 -20
Figure: ローレンツアトラクタ
26 / 29

本手法の特徴
本手法の特徴
スペクトル法と簡易ニュートン写像を使った精度保証付き数
値解法である.
基本解を精度保証付き数値計算によって得た.
1 つの時間ステップを大きく取れる.
27 / 29

今後の課題
今後の課題
より複数の時間ステップに渡って解を延長するために · · ·
解の挙動に応じて時間ステップの大きさを変化させる
端点評価において, KKT 方程式を用いて最適化
端点評価において, フローの評価を用いて g(s) の項の影響を
軽減
などを検討中.
28 / 29

ご静聴ありがとうございました.
29 / 29

狠狠撸

JSIAM_2019_9_4

Recommended

More Related Content

What's hot (20)

Similar to JSIAM_2019_9_4 (20)

JSIAM_2019_9_4