こんばんは
相変わらず前回の続きになります。
前回、前々回の記事はこちら↓
前回、前々回と1階の（非）線形微分方程式を数値解法して
その推定残差と真値残差（比率）を比較してみました。
今回は
本来の目的だった2階の非線形微分方程式を数値積分してみます。
非線形になると一般に厳密解は求められないので、
とりあえず
非線形要素のパラメータをゼロとおいて線形の微分方程式について
真値残差と推定残差を比較してみます。
そのあとで非線形微分方程式の数値積分解法の推定残差や
その他の解の性質をみることにします。
そういえば、映画「ジュラシックパーク」のなかで、
とある数学者が
「非線形微分法形式を知っているか？カオス理論は？」
というというようなセリフを言っていましたね。
監督のスピルバーグも
このような分野に興味があるのでしょうか？
さて、話題を元に戻して、
解きたい方程式は、Van del pol 方程式と呼ばれているものです。
これは
d^2y(t)/dt^2 + y(t) = μ(1-y(t)^2）dy(t)/dt
但し、μ>0
と表せます。
ここで、μ＝０とおけば
理系のひとにはなじみ深い
2階の線形微分方程式
d^2y(t)/dt^2 + y(t) = 0
になります。
一般的には
d^2y(t)/dt^2 + （ωｎ）＾２・y(t) = 0
のほうが
よく見る形ですね。
ここでは ωｎ＝１として話を進めます。
よく知られているように一般解ｙ（ｔ）は
ｙ（ｔ）＝ ｃ１・cos(t) ＋ ｃ２・sin（ｔ）
ただし、ｃ１，ｃ２は任意の定数
と書けます。
初期値：ｔ＝０のとき、 dy(0)/dt=0、y(0=１とすれば
ｃ１＝１、ｃ２＝０ になりますね。
よって この初期値での厳密解は
ｙ（ｔ）＝ cos(t)となります。
なお、y(t)の周期Ｔは、T=2π=6,28・・・ 秒となり
周波数ｆは、f=1/T＝1/2π Hz
ですね。
一方、ルンゲクッタ法による数値積分では
ｄｙ（ｔ）/dt ＝ ｘ（ｔ） とおいて
ｄｘ（ｔ）／ｄｔ ＝ -y（ｔ） と式を変形して
連立方程式として数値解を求めます。
ということで二つの解 ｘ（ｔ）とｙ（ｔ）が得られます。
厳密解ｙ（ｔ）を時間微分すれば
ｄｙ（ｔ）/dt ＝ -sin（ｔ） ＝ ｘ（ｔ） なので、
厳密解をまとめると
ｙ（ｔ） ＝ cos(t)
ｘ（ｔ） ＝ ｄｙ（ｔ）/dt ＝ -sin（ｔ）
になります。
これを図で表す方法は次のようにふたつがあります。
ひとつめは横軸に時間を縦軸に
ｘ（ｔ）とｙ（ｔ）の振幅を取ったものです。
これが図1です。
図1
また、別の表現として
ｘ（ｔ）を横軸にｙ（ｔ）を縦軸にとった
相空間として表すことができます。
これが図2です。リサジュー図形とも言いますね。
半径が１の円になります。
図2
では、ここから数値積分した結果です。
積分する刻み幅は時間軸ｔで t=0.1秒を初回の刻み幅にして
ｔの範囲は0から409.6秒までとします。
計算結果のステップも0.1秒刻みにします。
そのあとは、前回と同様、刻み幅を1/2、1/4、1/8・・・
と小さくして積分を繰り返します。
そして、計算されたｘ（ｔ）、ｙ（ｔ）と厳密解から
真値残差（比率）を求め
刻み幅との関係をグラフ化する手順になります。
まず、
積分の刻み幅と真値残差（比率）の関係は
ｘ（ｔ）とｙ（ｔ）についてそれぞれ図３と座4のようになります。
ます気が付くのは、右下がりのカーブの傾きです。
1階の微分方程式のときは、刻み幅が
1/2になると残差は(1/2)の4乗になっていましたが、
今回の図3、図4は
刻み幅が1/2になると残差も1/2)になっています。
したがって刻み幅を小さくしても、
残差は1階の方程式のように急激には小さくなりません。
さらに、図３と図4を比べると
ほぼ同じになっていることがわかります。
刻み幅は初回（=0.1秒)のものから、1/2、１／4、1/8、・・・
とし最後は1/（2のお23乗)までにしました。
エクセルのマクロの処理ループさえ増やせば
いくらでも小さくできますが・・・
とりあえず10数時間程連続で動かしてここで止めました。
図3
図４
ついでに今回は、
3種類の刻み幅に対する相空間のグラフ（図５，図6、図7）
も追加してみました。
図３のグラフの赤い楕円がその刻み幅に該当します。
この３つは、
積分区間0から409.5秒の範囲の400から409.5秒まで、
つまり
残差が最も大きいであろう積分区間の最後の部分としてみました。
図５から図７の青いカーブは厳密解＝半径１の円です。
まず、図5は刻み幅＝0.1秒としたときの数値積分した結果です。
原点の青い点が厳密解で・・・
この刻み幅でこの区間まで積分すると
厳密解とは全くかけ離れた解になります。
真値残差（比率）の最大値が1e5なので当然か・・・
図５
図6は刻み幅を小さくして0.1秒の1/（２の5乗）、
つまり、（0.1/32）秒としたものです。
真値残差(比率）がまだ大きいですね。
図6
さらに刻み幅を小さくして図６の刻み幅の1/32つまり(01./1024)秒
としたものです。
このときの真値残差（比率)の最大値は約０．５です。
見た目は厳密解にかなり近くなりましたが、
グラフでもまだ違いが見えます。
図7
ちょっとブレイクして、図5の動きが面白いので
時間ｔについて全積分区間の相空間を見てみます。
図８は、解としては意味を持ちませんが、
模様あるいはアート（？）
として使えそうかもですね。
図８
では、話をもとに戻して
刻み幅が一番小さいときの解と真値残差（比率）の
相関を図にしてみました。（図９から図１１）
図３と図４の横軸の一番右側の刻み幅に対応します。
真値残差（比率）の最大値はｘ（ｔ）、ｙ（ｔ）とも約1e-5です。
図9はx(t)とその真値残差（比率）です。
横軸はx(t)=-sin(t)で t=0のときｘ(0)=０から始まり（図の赤い矢印）
-1→0→１→0→-1・・・とサインカーブで振れます。
グラフの中央ｘ＝０の周辺が開いているのは、
真値残差（比率）の分がゼロになるので、
定義できず除外しているためです。
また、x=0周辺では真値残差（比率）の分母がゼロに近づくので、
黄色い矢印「のように増加します。
さらに数値積分自体が累積加算のため、
緑の矢印のように積分が進行するにつれて
残差も増加していきます。
図９
図１０は積分の開始位置（赤い矢印）から
約2サイクル弱程度を取り出したものです。
図１０
図１１はｙ(t)とその真値残差(比率）です。
図１１
図１２はその部分拡大です。
図1２
ｘ（ｔ）とy(t)の真値残差の傾向は、
初期値以外ほぼ同じになることがわかりました。
図１３はx(t)とy(t)の相関をとったものです。
図の2か所のオレンジの丸印のところにピークがあり、
ｘ（ｔ）とｙ（ｔ）が直交していて、
他方が極大若しくは極小値をもつとき、
つまり時間に対する傾きがゼロのとき
他方の時間に対する傾きが最大になり、
シーソーのような動きになります。
図１３
図１４は対数表示にしたものです。
ピークになっているところ以外は特に特徴がなく
まんべんなく広がっているように見えます。
図１４
真値残差(比率）のようすが概略わかったので、
次は推定残差（比率）です。
ｘ（ｔ）とｙ（ｔ）の推定残差（比率）はぞれぞれ図１５と図１６です。
最大値の横軸が水平になっている部分があります。
これは判定条件を満たさないデータは除外したため
結果としてこの形になりました。
いづれにしてもこの部分のデータは使用しないので、
あまり気にしないことにします。
図１５
図１６も同じ傾向です。
図１６
真値残差（比率）と相対残差（残差）を比較したものが
図１７から図２０です。
図１７を見ると、両者は一見一致しているようにみえますが・・・
縦軸を拡大して違いがわかるようにすると、図１８になります。
図１７
図１８を見ると、
真値残差と推定残差には
刻み幅の一目盛り分の差異がありますね。
これを一致させるのは、
推定残差を刻み幅ひとつ分シフトすれば良さそうです。
図１８
ｙ（ｔ）についても同じ傾向で、
図１９、図２０のようになります。
図１９
図２０
ということで、
ｘ（ｔ）について、図１８の丸印の刻み幅のふたつの残差について
その相関を調べてみました。
すると驚くべきというか拍子抜けというか・・・
図２１のように傾きが１の直線で相関係数１になってしまいました。
図２１
念のため、
エクセルで最小2乗法を使って直線の傾きと切片を求てみると
ｙ＝0.999998 ｘ - 3.5e-13
≒ ｘ
となり、一致していると言っていい結果になりました。
ｙ（ｔ）についても調べてみると、図２２のように同様な結果です。
図２２
同じ計算をやってみると
ｙ＝ 0.999998084 ｘ - 4.25004e-13
≒ ｘ
となりました。
つまり、推定残差(比率）は真値残差(比率）と等しいので、
真値残差(比率）について調べたことは
すべて推定残差(比率）に適応できることがわかります。
ということで、
この2階の線形微分方程式でも
推定残差(比率）と真値残差（比率）は
ほぼ一致していることがわかりました。
せっかくエクセルを使って解析をやっているので、
ついでに
この解の周波数をエクセルの機能を使って求めてみます。
実は時間データ数４０９６個はこのために設定したものでした。
メニューの「データ分析」の中の「フーリエ解析」機能を
マクロプログラムに取り込んで実行します。
FFTを実行するので
入力データ数は２のべき乗であることが前提で、
最大の入力データ数は4096個になっています。
エクセルの入出力セルを指定すれば
複素スペクトルが出力されます。
が・・・
残念ながら、マクロには複素数を扱う機能は標準ではありません。
対応方法は調べればいろいろありますが、
今回は複素数の絶対値を求めるだけなので、
これだけはマクロプログラムから文字列として
セルに直接数式を書き込んで、
セルの処理として実行させる方法をとりました。
文字列の操作はわかってしまえば意外と簡単でしたよ。
FFTの実行前にはハニング窓を掛けてあるので、
周波数軸上ではメインローブが広がります。
また、スペクトルの計算結果の値はFFTの原理上、
時系列データ数だけ積算され（結果として今回は４０９６倍され）
またスペクトルの折り返しで1/2になり、
さらに窓関数のため1/2になります。
したがって、
今回は4096/2/2＝1024の値で
得られた結果を割り算しています。
ということで、ｘ（ｔ）についての結果の縦軸を
デシベルで表現したグラフが図２３です。
ｘ（ｔ）、ｙ（ｔ）のどちらの振幅も１なので、
デシベルで表せば、20*log10(1)＝０ （dB）となります。
必要なデータは周波数領域でこの左半分ですが、
とりあえず全データを載せてみました。
横軸は周波数で数字は
基本周波数ΔFの整数倍の値にしています。
ここでΔF＝サンプリング周波数/時系列データ数
＝(1/0.1)/4096
＝0.002441・・・ Hz です。
本来の理論周波数は 1/２π なので、これをΔFで割ると
Ñ＝ 65..189・・・ になります。が
FFTでは横軸の周波数は整数となるのでピークは６５になります。
厳密解と数値積分した時系列データに対するスペクトルの両方を
載せましたが、両者は完全に重なってしまいました。
当然と言えば当然ですが高調波は見当たらず、基本波のみですね。
赤い丸印部分を拡大したものが図２４です。
図２３
図２４
さらに、窓関数の特性を考慮して
理論上のカーブを入れてみました
図２５の青い矢印が、
左から順に-6dＢ、０DB、-6dＢとなる周波数です。
グラフ上での確認なので定量的なことまではわかりませんが
傾向は同じですね。
図２５
なお、
ｙ（ｔ）はｘ（ｔ）と全く同じ結果になったので、ここでは省略します。
今回はｘ（ｔ）とｙ（ｔ）の2種類の解があるので、両方を調べました。
データ量も2倍になって、ちと大変でした。 (^^;