新潟工科大学紀要 第 8 号 2003 年 12 月

単独保存則方程式の周期解の数値解析 II
竹野 茂治 ∗ , 小松 幸恵 ∗∗(平成 15 年 10 月 31 日受理)

Numerical results of time–periodic solutions
for a scalar conservation law IIShigeharu TAKENO∗ and Yukie KOMATSU∗∗
In the last report, we mentioned some results for periodic solutions of hyperbolic scalar conservation law equations with a periodic outer force term.
The structure of periodic solutions is complicated and includes period doubling bifurcation. In this paper, we introduce the relation between periodicsolutions and the effect of viscosity, and show some properties for periodic
solutions by using simple examples and numerical computations.

Keywords: scalar conservation law, Burgers equation, time periodic solution, period doubling bifurcation, numerical computation

1

はじめに以前報告 4) したように，時間周期的な外力を含む単独の保存則方程式の境界値問題(ut + f (u)x = g(t, x) (0 < x < 1, t > 0),u(t, 0) = u(t, 1)(t > 0)00(u = u(t, x) ∈ R, f (u) > 0, g(t + T, ·) = g(t, ·))(1)の周期解の様相は単純ではなく，解明されていない多くの現象を引き起こす．例えば周期解の一意性が保証されず，外力の周期の 2 倍，3 倍の周期を持つ周期解が存在することが，それら多くの問題を引き起こしていると考えられる．
初期平均値をū =***

Z 10u(0, x)dx

情報電子工学科 助教授
舞鶴工業高等専門学校 助教授

竹野 茂治 · 小松 幸恵とすると，外力項がある種の条件を満たす場合に，ū の任意の値に対し外力と同じ周期(T ) を持つ周期解が少なくとも一つ存在することが知られている 1) が，前回の報告 4) で見たように，適当な初期値からスタートした解の漸近的な収束先として得られる，ある意味で安定な周期解は必ずしも外力とは同周期ではなく，外力の 2 倍，3 倍などの周期を持ち得る．その場合でも外力と同周期の解の存在は保証されるので，周期解は外力と ū に関して一意的には決定しないことになる．
本報告では，Burgers 方程式 (f (u) = u2 /2) に対して，前回の報告 4) で得られた事実や予想などに関して，その後に行なった解析や数値計算によってわかってきた周期解の性質，特に倍周期解の存在や非存在と粘性との関係，周期解の集合の構造などについて述べる．

2

粘性方程式の周期解

本稿や前回の報告 4) では，保存則方程式の数値計算では段階的な Lax-Friedrichs 法と呼ばれる以下の差分法を用いている．
n+1/2uj=

unj+1/2 + unj−1/2
2
n+1/2
n+1/2
uj+1 + uj−∆t{f (unj+1/2 ) − f (unj−1/2 )},
2∆x

∆t
n+1/2
n+1/2
n+1/2
{f (uj+1 ) − f (uj)} + gj+1/2 ∆t22∆x(n ≥ 0, 0 ≤ j ≤ L − 1, L∆x = 1, unj = u(n∆t, j∆x))

un+1j+1/2 =−この差分法は，(∆x2 /2∆t)uxx の形の誤差を持ち，よって以下の形の粘性保存則方程式の解に似た性質を持つ．
(ut + f (u)x = g(t, x) + εuxx
(0 < x < 1,u(t, 0) = u(t, 1), ux (t, 0) = ux (t, 1) (t > 0)(ε = ∆x2 /2∆t)t > 0),(2)ここで，比 ∆x/∆t は CFL 条件によりほぼ一定に固定し，よって ε = ∆x2 /2∆t = c∆xと考えることができるので，∆x がある程度小さくないと数値計算結果に粘性効果による影響が現われる．
前回の報告 4) においても 4,5 節においてその粘性効果の影響について実験，考察し，∆x が小さい場合，すなわち粘性効果が小さい方が倍周期性がはっきり現われるという数値計算結果を紹介し，逆に粘性効果は倍周期性を消す方向に働くのではないか，と予想しているが，その後，倍周期性と粘性効果の関係を進行波解 u(t, x) = w(x − st) (s は定数)に関して数学的に考察することで，次の結果を得た．
定理 1 f (u) = u2 /2 の非粘性 Burgers 方程式 (1) では，任意の整数 N > 1 に対し，外力の周期の N 倍の周期を持つ進行波型 u(t, x) = w(x − st) の周期解と g(t, x) の例を構成できる．
定理 2 f (u) = u2 /2 の粘性 Burgers 方程式 (2) では，N > 1 に対して，外力の周期のN 倍の周期を持つ進行波型 u(t, x) = w(x − st) の周期解は存在しない．

単独保存則方程式の周期解の数値解析 II
これらはあくまで進行波解に制限した形の結果であるが，粘性のついた方程式ではその粘性効果によって倍周期性が起こらず，倍周期性は非粘性方程式に特別な性質であることが分かる．
次に，粘性以外の効果が働いた場合として，さらに高次導関数項がついた方程式µ∂
ut + f (u)x = g(t, x) + a∂x¶ku(a 6= 0,k ≥ 3)

(3)

ではどうなるかを考察してみた．k = 3 の場合の高次付加項は特に分散項と呼ばれ，いわゆる KdV 方程式と同等である．この場合も進行波解に関しては次の結果が得られた．
定理 3 方程式 (3) では，外力の周期の 2 倍の周期を持つ進行波型の周期解と g(t, x) の例を構成できる．
高次の項がついた場合は，粘性の場合とは異なり倍周期性を持ち得ることになり，よって，その倍周期性を消すのは粘性効果のみであることが，進行波解に関しては示されたことになる．
ただし，進行波解とは限らない一般の周期解に対して，粘性方程式が倍周期性を持たないことはまだ証明されてはおらず，現在のところまだ未解決である．

3

周期解の個数次に，非粘性方程式の周期解の集合の構造について考察する．
周期解は，外力と，周期毎に保存される値である初期平均値ū =Z 10u(0, x)dx

に関して一意的には決まらないことが，倍周期性を保証する定理 1 と，外力と同じ周期を持つ周期解の存在を保証する結果 1) に関する簡単な考察からわかる．
例えば v(t, x) が，方程式 (1) の，平均値 ū を持つ 2 倍周期解であるとする．

v(t + 2T, ·) = v(t, ·),v(t + T, ·) 6= v(t, ·),Z 10v(0, x)dx = ū

一方で，同じ平均値 ū を持ち，外力と同周期の周期解 w(t, x) の存在が保証されている．

w(t + T, ·) = w(t, ·),Z 10w(0, x)dx = ūv,w はもちろん等しくはないが，g の周期を 2T と見れば，この 2 つは同じ平均値を持ち，外力の周期 2T と同じ周期を持つ周期解になっている．よって，周期外力，および平均値 ū だけで周期解は一意には決まらないことになる．
ここから，一つの外力と一つの平均値 ū に対して，はたしてどれくらいの数の周期解が存在するのかという疑問が起こるが，前回の報告 4) の 5 節の例がはからずもその答えを与えていることがわかった．その例を，平行移動し構造を簡単にしたもので紹介する．

竹野 茂治 · 小松 幸恵
外力がµtg(t, x) = πa sin 4π x −2T¶2(a > 0)(4)の場合 (周期 T )，方程式 (1) は次のような 2T 周期を持つ滑らかな周期解を持つ (Fig.
1)．
µ1tu(t, x) =+ a sin 2π x −2T2T¶(5)

一方で，以下の関数が任意の実数 0 < α < 1/2 に対して，同じ平均値 ū = 1/2T を持つ
周期 2T (α = 1/4 なら周期 T ) の不連続周期解であることを容易に検証できる (Fig. 2)．

u(t, x) =1+ aH(x − t/2T ; α, 1 − α) sin 2π(x − t/2T )2T(6)

ここで，H(ξ; a, b) は，以下のような周期 1 の符号関数とする．
(H(ξ; a, b) =−1 (a < ξ < b),+1 (0 < ξ < a, b < ξ < 1),

H(ξ + 1; a, b) = H(ξ; a, b)u(t,x)u(0,x)u(t,x)

t=0
t=T1/2T

00.5x1

Fig. 1 a periodic solution with period 2T(t = 0, t = T )1/2T0

α0.5x1-α1Fig. 2 a periodic solution having a freeparameter α (t = 0)

なお，不連続な部分も含む関数が保存則方程式の解であるためには，滑らかな部分が方程式を満たすこと以外に，不連続な部分が Rankine-Hugoniot 条件，およびエントロピー条件と呼ばれる条件を満たす必要がある 2) ．
この例は，一つの外力と一つの平均値に対する周期解が，少なくとも 1 パラメータの自由度を持つ程度に存在し得ることを示している．

4周期解の構造に関する実験

同じ平均値を持つことにこだわらなければ，前節の例 (6) を少し変化させた以下の関数も，0 ≤ α ≤ 1/2 ≤ β ≤ 1 に対し同じ外力に対する周期 2T (α + 1/2 = β なら周期 T )の周期解であることが分かる (Fig. 3)．

u(t, x) =1+ aH(x − t/2T ; α, β) sin 2π(x − t/2T )2T(7)単独保存則方程式の周期解の数値解析 II
これらの解はあたかも，周期 T の極小解u(t, x) =1− a| sin 2π(x − t/2T )|2T(8)と，周期 T の極大解u(t, x) =1+ a| sin 2π(x − t/2T )|2T(9)(Fig. 4) を部分的に取ってそれを衝撃波でつないだような形になっている．

u(0,x)

u(0,x)u(t,x)

1/2T0

α0.5
xβ1

Fig. 3 a periodic solution having two parameters (t = 0)

1/2T00.5
x1

Fig. 4 minimal and maximal periodicsolutions (t = 0)

連立の保存則方程式である気体の方程式の場合に得られた倍周期解の T 毎の関数の形がかなり違っている 3) のに対して，Burgers 方程式の数値計算で得られた数値的な倍周期解は，いずれも T 毎の関数の形にあまり違いがなく 4) ，• 違うのは衝撃波の付近だけで，それ以外の滑らかな部分は完全に一致する• 違うのは衝撃波の位置のみという際だった特徴が見られる．これは (7), (8), (9) の構造に良く似ている．
もちろん，我々の手法，すなわち，適当な初期値からスタートして解を先の方まで追いかけることで漸近的な収束先としての周期解を得るという手法では全ての周期解を捕捉することができるわけではなく，少なくとも数値的に安定なものしか得られないので，上記の性質を持たない周期解が存在する可能性はある．実際，定理 1 で構成した解はこのような性質を持たない解である．
しかし，我々の手法で捕捉可能な安定な周期解は上記の性質を持つようであるので，それらが (7), (8), (9) のような構造から来るものなのかを検証するための数値実験を行なった．目標は

1. 安定な周期解は常に 1 パラメータの族をなす位存在するか
2. もしそうであれば，そのような族の構造は (7), (8), (9) の構造と同様であるか竹野 茂治 · 小松 幸恵
を調べること，とした．
我々の手法では，同じ初期値からスタートすれば同じ周期解にしか漸近し得ないので，別な周期解を得るために，同じ平均値を持つ別な初期値を使ってその漸近収束先を追ってみた．その結果が Fig. 5, 6 である (グラフはいずれも u(mT, x) を，30 個分重ねたもの)．
外力は前回の報告 4) と同じで，Fig. 5 は初期値を定数 u(0, x) = 定数 = ū としたもの，

 0.65

u(mT,x)u(mT,x)0.65

0.60.55

0.60.55
00.5x

1Fig. 5 a solution for constant initial data(m = 300 ∼ 330)00.5x1

Fig. 6 a solution for a step value initialfunction (m = 500 ∼ 530)

よって前回の報告 4) の 4 節に紹介したものと同じもので，Fig. 6 は，平均値は ū のまま
初期値を x = 0.5 で段差を持つ階段関数に変えた場合の結果である．これによりわかる
ように，確かに初期値によって漸近形が違っていて，収束先が初期値によらないという気体の連立方程式の場合の結果 3) とも異なっている．
なお，この違いが ū に関する量子誤差によるものでないことは確認済みであり，実際量子誤差よりも大きめの誤差を与えて数値計算したところ，漸近形に対する影響は，Fig.
5, 6 の違いよりもごくわずかなものしか見られなかった．
さらに初期値を色々変えて実験を行なうと，衝撃波の位置の違うものが色々得られ，それらが Fig. 5, 6 を含んで連続的に変形するような 1 パラメータの自由度を持つような族を形成していることが確認できた．
次に，その周期解の集合の構造を調べるために，初期値は定数 (= ū) にして，その ūの値を少しずつ変化させて極限の周期解を求めてみた．それが Fig. 7–10 である (グラフはいずれも u(mT, x) (m = 300 ∼ 330) を重ねたもの)．
ū を細かくより連続的に変化させるとわかるが，ū=0.5891, 0.6295 では衝撃波のない1 倍の周期解 (極小解，極大解) が得られていて，その間の ū ではその両者を衝撃波でつないで平均値を調整しているような倍周期解が現われる．
また，実際には ū=0.6089 付近で，極大解と極小解が接触するところで交差するグラフを持つような，衝撃波を持たない滑らかな 1 倍周期解が存在するようにも見えるが，正確には確認できていない．
その後，ū を 0.6295 よりもさらに大きくした場合を調べてみたが，しばらくはあまり漸近安定性の良くない 1 倍の周期解のようなもの (前回の報告 4) 3 節の Fig. 2 参照) が現われ，その後また新たな極小解と極大解を与える ū の値のペアが現われ (例えば 0.728と 0.739 や 0.805 と 0.827)，その間の ū ではその両者を衝撃波でつなぐ倍周期解が現わ

単独保存則方程式の周期解の数値解析 II  

m.v.=0.594
0.65u(mT,x)u(mT,x)0.65

0.60.55

0.60.55
00.5

10x

0.5

Fig. 7 ū =0.5891, 0.6295 (minimal andmaximal solutions)Fig. 8 ū=0.594

m.v.=0.617

m.v.=0.626
0.65u(mT,x)0.65

u(mT,x)1x0.60.550.60.55
00.5x10Fig. 9 ū=0.617

0.5x1Fig. 10 ū=0.626

れるようである．

5考察

前節の実験結果をまとめると，以下のようになる．

• 倍周期解の構造は (7), (8), (9) と似た構造をしていて，滑らかな 1 倍の極小解と極大解が存在し，その間を衝撃波でつなぐような倍周期解が起きている．
• 倍周期解は，極小解と極大解の間の平均値に対しては (6) と同様，衝撃波の位置に関して 1 パラメータの族を作るくらい存在する．
• 極小解と極大解は，N 箇所で接触し，倍周期解は，それによって分離された N 個の小区間にそれぞれ衝撃波を持ち，その衝撃波で極小解と極大解をつないだ形の解になっていて，これが N 倍周期解を形成する (N ≥ 2)．
• 極小解と極大解を与える平均値の組 (区間) は多数存在し，それらの区間内では倍周期解が発生し，それらの区間外では前回の報告 4) の 2 節に見られたように，周期解への漸近はやや不安定になる．

竹野 茂治 · 小松 幸恵
今回実験に使用した外力の周期 T は T = 0.8 であり，実験で見た 2 倍周期解の極小解，極大解に対しては(極小解: ūT = 0.5891 × 0.8 = 0.471
極大解: ūT = 0.6295 × 0.8 = 0.504となっていて，いずれも 1/2 に近い値になっている．Burgers 方程式では，波は速度 uで移動し，よって平均して ū の速度で移動しているので，ūT ≈ 1/2 は，1 周期分の時間T で波が x 方向に [0,1] 区間を約半分だけ移動することを意味する．
つまり，ū = 1/2T の近くで 2 つの衝撃波を持つ 2 倍周期解が起こるのは，1 周期毎にほぼ互いの衝撃波の位置が入れかわり，それらがわずかにずれているためであると予想される．その付近の ū では極小解と極大解の存在により，その両者を適当にまたぐことで平均値を調整できるので，安定な周期解が存在する区間を保持できるのだろうと考えられ
る．ū = 1/2T の他にも，ū = 2/3T , 3/5T などの近くではそれぞれ 3 倍，5 倍の倍周期性が見られるので，倍周期数はその移動速度 (特性速度) との関係が深いようである．
一方，それらの区間外での ū に対しては，1 周期時間での波の移動距離と x の有限区間長が簡単な整数比にならないので，それにより 1 倍にもなりにくい，またはかなり大きな倍数の倍周期解にまとまろうとしている，などの理由で不安定な状態に見えるのだろうと予想される．
安定な倍周期解が，極小解と極大解をつなぐ衝撃波を持つ，1 パラメータ族をなすことは単独の保存則方程式特有の性質ではないかと思われる．単独保存則方程式の場合，衝撃波 x = x(t) が満たすべき条件である Rankine-Hugoniot 条件はx0 (t) =f (u(t, x(t) + 0)) − f (u(t, x(t) − 0))u(t, x(t) + 0) − u(t, x(t) − 0)の 1 本の方程式であり，もう一つの条件のエントロピー条件は右に向かって跳び下がることなので，大小関係のある滑らかな解が 2 つ存在する場合，その 2 つの解を接続するような衝撃波は，どの出発点 x(0) に対しても上の常微分方程式を満たすように必ず 1 つ作れることが常微分方程式の理論よりわかる．
つまり，単独保存則方程式の場合，そのように容易に作れる衝撃波を利用して平均を調節することが可能で，そのことがそのような解が安定でかつ現われやすい，という理由になっているのではないかと思われる．
連立方程式の場合は Rankine-Hugoniot 条件は複数の条件になり，それを衝撃波に対する常微分方程式と見ても過剰決定系になってしまい，解の存在は明らかではなく特定の状況の元でしか解は存在し得ない．よって，そういう解が容易に安定なものとしては現われないのだろうと考えられる．
また、今回得られた極小解と極大解は，丁度その接点の位置に関する空間対称性を持っているように見える (Fig. 7)．これは，2 倍周期解の場合だけでなく 3 倍周期解の極小解と極大解の組に対しても見られる．それが，現在数値計算で用いている周期外力の空間対称性が原因なのか，それ以外の何らかの性質を表しているものなのかは今のところは不明であり，今後それを検証する数値計算も必要であろうと考えられる．

単独保存則方程式の周期解の数値解析 II

6

まとめ

今回の数値実験や考察によって，(1) の周期解の様子はかなり明らかになってきた．特に前回の報告 4) で今後の課題としてあげたことの多くはその概要が見えてきたと言える．
しかしいずれも数値的に得られた具体的な例ばかりで，一般性を持って厳密に証明された事実はまだ少ない．今後，今回明らかになった事柄を含む以下の数学的な解析を進めていきたいと考えている．

• 一般に粘性保存則方程式が倍周期性を持たないこと• 一般の周期外力に対する極小解と極大解の存在，および 1 パラメータの自由度を持
つ倍周期解の存在• 非倍周期区間での漸近極限の解析 (周期解への遅い収束 ?)• 気体の連立方程式の場合の周期解の構造の解明
そのためには，さらなる数値実験が必要であろうし，厳密な解析にも新たな数学的手法が必要ではないかと考えている．

参考文献
1) S. Takeno: “Time-periodic solutions for a scalar conservation law”, Nonlinear Anal.
45, pp1039-1060 (2001)
2) J. A. Smoller: Shock waves and reaction-diffusion equations. 2nd ed. (Springer,1994)
3) 相京 敦: “理想気体の 1 次元流の数値解析”, 新潟大学工学部卒業論文, pp1–20 (1994)4) 竹野茂治: “単独保存則方程式の周期解の数値解析”, 新潟工科大学紀要, 2, pp19–26(1997)