こんにちは
「その３」に続いて今回はいよいよ、
μ≠０ のときの方程式の解を考えてみます。
ここからが、とてもおもしろく（？）なってきます。(私だけか？）
あらためて微分方程式は
Van del pol 方程式と呼ばれているものです。
これは
d^2y(t)/dt^2 + y(t) = μ(1-y(t)^2）dy(t)/dt
但し、μ>0
と表せます。
ちなみに、
この方程式の歴史的な話や全体像の解説が
ネット検索したらPDFで見つかりました。
情報はこちらです。↓
「 電気学会雑誌 J. I. E. E. J. 7/'73 技術総説
電気工学における非線形問題 」
電気屋さん、特にハード屋さんには
なじみがあるかもしれません。
話を戻します。
まず、
μを一定値に固定して初期値を変えた場合の動きを調べてみます、
いくつでもいいのですが、今回はμ＝0.5 としました。
線形の微分方程式では初期値が変われば
その解も変わるのが当然なのですが、
この方程式の場合全く異なる動きになります。
とりあえず、
初期値を色々変えて数値積分して
その解を相平面に表したものが図1です。
積分の刻み幅は
0.1/（2の13乗）=0.1/8192のものを使用しています。
あとでグラフが出てきますが
推定残差（比率）の最大値は
およそ1e-3、最小値はおよそ5e-7です。
あくまで「その３」の結果から類推したものですが、
かなりザックリ言って
解の桁の上から４桁目が怪しいことになります。
（やっとこれを書くことができました。）
これを知ることが私の本来の目的なのでした。
(^.^)
「その１」から「その３」の内容は
たったこれだけを知るためだったのですが・・・
知るのと知らないのは大違いなのであります。
では、話を本筋に戻します。
図１
グラフにマーカ（太い丸印）があると太すぎて見づらいので、
直線だけの表示にしたものが図２です。
＊＊＊＊＊
あ、前回の記事｛その３」はこちらです。↓
サムネイル（？）の画像にどうしても図１を入れたかったので、
ここに入れました。必要のないかたはスルーで。 (^^;
＊＊＊＊＊
図２
図１(図２）では6種類の初期値を適当に選んで、その解を計算し
x(t)とy(t)をプロットしたものです。
（ひとつは別のカーブに完全に重なって見えませんが・・・）
原点を除いた任意の初期値ｘ(0)、ｙ(0)がら初めても、
時間が経過するにつれて同じ周期解に収束しています。
まるで巻き付いているようですね。
これがこの非線形微分方程式のおもしろい特徴のひとつです。
この周期解はリミットサイクルと呼ばれ、
色々な非線形現象にも出てきます。
単純に数値積分をやった
（正しくはPCがエクセルのマクロのプログラムを実行した）
だけですが、とりあえずリミットサイクルは確認できました。
よかった 、よかった (^^;
電気・電子系のかたはご存じかも知れませんが、
有限桁の数値計算でフィルタリングする
ＩＩＲデジタルフィルタにも
リミットサイクルが発生する場合があることがわかっています。
普通は好ましくないのですが、
現在も有効な解決策はないようです。
40年ほど前にデジタルフィルタを設計する機会があり、
そのときに初めてこの現象を知ったのですが、
いまだに未解決のようです。
結局ＩＩＲは止めてFIRデジタルフィルタにした経験があります。
若干話がそれました。
もとに戻します。
相平面では解の波形がピンとこないので、
時間変化をグラフで表してみます。
6種類の解のなかから変化の大きい2種類を選びました。
これが図３と図4です。
ふたつとも赤い丸印のところが初期値で、
ここから急激に変化して周期解になるようすがわかりますね。
図３
図４
さて、
肝心の残差（比率）ですが、
この非線形微分方程式は解析的な厳密解はない
（近似的に求める方法はあります）ので、
推定残差（比率）だけになります。
これを6種類の解について求めてみました。
6種類のｘ（ｔ）、ｙ（ｔ）はほぼ同じ形になったので、
ここではその中の一つを図５に示します。
図５
図５の赤い丸印の刻み幅（最も推定残差（比率）が小さいもの）
の解を図１(図２）に使いました。
ということで、
６種類の解の推定残差（比率）をひとつにまとめたもの
図６です。
最小値に凸凹がありますが、
桁的にはほぼ同じなのがわかりますね。
図６
では、
ここで条件を変えて
初期値を一定にして、
パラメータμを変化させたときのようすをみてみます。
初期値はx(0)＝0、y(0)＝1の一定値で、
μは0.1から0.9まで、0.1刻みで変えてみました、。
その前に各々の推定残差(比率）を確認しておきます。
図７と図８がそのグラフです。
最大値、最小値のいずれも若干の変動はありますが、
桁的にはほぼ同じと言えますね。
図７
図８
それでは
相平面のグラフが図９、図１０、図１１です。
9種類の解をひとつにまとめると見づらいので
3種類をまとめて、3枚のグラフにしました。
（マーカも見づらいので省略し、直線だけにしています。）
図９
図１０
図１１
μが増加するにつれて、
リミットサイクルが次第に円から形がくずれ、
歪んでいくのがわかりますね、
ただ、ｘ（ｔ）つまり横軸方向は歪みが大きいのですが、
ｙ（ｔ）＝縦軸方向の歪みは横軸ほど大きくはないですね。
さらに、縦軸ｙ（ｔ）の最大値(=振幅）はμの値に無関係に
ほぼ２になります。
この方程式の解のおもしろい特徴です。
（実は若干増加しますが、これはこのあとで調べます。）
相平面よりも時間変化のほうがピンとくるので
初期値（＝0秒）から30秒までの時間を横軸とした
ｙ（ｔ）の３枚のグラフを図１２、図１３、図１４に示します。
図１２
図１３
図１４
最初の30秒は、さほと大きな差は見当たりませんが・・・
今度は終わり辺りの時間、
400秒から409.5秒までのｙ（ｔ）の時間変化を見てみます。
これが図１５，図１６、図１７です。
図１５
図１６
図１７
この3枚は明らかに違いますね。
波の形が時間的にずれています。
これを位相がずれていると言い、
一般的に連続した周期波形では周波数が違うとこうなります。
また、
波形の縦軸の大きさ（振幅）はどれも
±２でグラフでは同じように見えます。
（マーカが太いのでやや見づらいですが）
ついでに
ｙ（ｔ）の時間微分ｘ（ｔ）の時間変化も400秒から409.5秒まで
グラフに表してみました。
これが図１８、図１９、図２０です。
図１８
図１９
図２０
これら３つのグラフのカーブは、
個人的には全く違和感がなく、どこかで見たことがある(気がする）
波形ばかりです。
どこで見たかは思い出せませんが・・・
ではでは
ここからはμが変わるとy(t)の周波数が変わっていることと
振幅が２から僅かに増加することの
ふたつについて調べてみます。
ます、μに対する周波数の変化です。
図１２から図１７のような時系列データｙ（ｔ）について、
各々のμことにエクセルでフーリエ解析を行い
周波数分布を求めたものが図２１です。
横軸はスペクトル番号Nで、
周波数ｆとの関係は
分解周波数Δｆを
Δｆ＝ 時系列のサンプリング周波数 / 時系列データ数
＝ (1/0.1) / 4096
＝ 0.00244140625 (Hz)) （＝一定値）
としたとき
ｆ ＝ N ＊ Δｆ
となります。
ここでのフーリエ解析では、
Nは 0,1,2,3,・・・,2045,2046,2047の自然数で、
Nから求めた周波数ｆはΔｆに比例し離散化された値となります。
縦軸は振幅レベル（デシベル）で表しています。
9種類のμについて1枚のグラフにまとめました。
これだけでは概略しかわかりませんが・・・
図２１の赤い丸印が最も低い周波数（基本周波数）で、
その奇数倍の高調波を含んだ周波数分布になっていますね。
また、高調波の減衰度（対数表記）は、
周波数に対し直線的に変化しています。
普段よく見るパターンは周波数（横軸）も対数表示にしたときに
縦軸が直線的に変化するのですが・・・
この形も特徴になりそうです、
図２１
図２１の基本周波数部分を部分拡大したものが図２２です。
どの解の振幅もほぼ２なので、
グラフでは20*log10(2)=６デシベルにピークを持ちます。
μが増加するに従い
ピークの周波数が小さくなっているのがわかります。
図２２
ここで、μに対する周波数変化を示す2次の近似式を
手持ちの本から見つけました。
基本周波数をｆbとすると
ｆb ＝｛ １－ (μ＾2)/16 ｝ / 2π
となります。
これをスペクトル番号Nに対応させると
N = ｆb / Δｆ
になり、プロットしたものが図２３になります。
μが増加するとNが減少します、
図２３
さらに
図２２を縦横の各方向に拡大して、
図２３のNの値を追加したものが図２４です。
比較が簡単なように
図２３のNは縦軸の値が10dＢの直線としました。
図２４の赤い矢印がμ=0のときのNになり、
緑色の矢印が図２３のNと
そのときのμに対応する周波数分布のピークとをつないだものです。
矢印がまっすぐ下に降りていればふたつのNは同じになります。
図２４
グラフでの比較なので定性的になりますが、
μが小さいときは
フーリエ解析したスペクトルのピーク（エクセルでの近似カーブ）と
図２３のNは殆ど違いはなく
μが大きくなるにつれて差が大きくなる傾向です。
図２３は2次の近似式なので、
おそらくμが増加するにつれて確度が下がるためだと思いますが、
詳細は調べていません、
ということでこの話題はここで終了です。
では、次の話題です。
図９，図１０，図１１のところにも書きましたが、
この方程式の解ｙ（ｔ）の周期解の振幅は
約2で微小増加することがわかっています。
たまたま、ネット検索で、
とある修士論文要旨にμと振幅の数値データを見つけました。
μの範囲は0.1から0.9まで刻み幅が0.1だったので、
実は、これに合わせて数値積分をやってみた次第です。
（そのほかにも色々と条件を変えては計算しています。）
では、
入手したデータ（これを理論値とします）の値はいくつかというと
μが0.1のとき 2.00010398
：
：
μが0.9のとき 2.00724521
という殆ど２と変わらない程度のものです。
このグラフが図２５です。
図２５の横軸がμを、縦軸がｙ（ｔ）の周期解の振幅です。
ここに数値積分した結果
（2種類：数値積分結果の最大値と最小値の絶対値））も入れたので
3本のカーブになりますが、
完全に重なってしまい1本に見えます
図２５
では、ここで各μごとに残差を比較します。
今回は真値（＝理論値）がわかっているので
真値残差(比率）が計算でき、推定残差(比率）とも比較できますね。
図２６は、図８のｙ（ｔ）だけを表記したグラフに
振幅の真値残差(比率）を追加したもので赤い丸印で囲みました、
とりあえず最大値と最小値の間にあるので、
おかしくはありませんね。
図２６
次に、真値残差(比率）と推定残差（比率）の比較です。
これが図２７です。
それぞれ振幅のプラス側とマイナス側に分けて、
4種類のデータになります。
横軸はμの値で、縦軸は2種類の残差(比率）です。
図２７
4個の残差（比率）は
μ＝0.9のところが違いがあり
その他も多少の凸凹はありますが
ザックリ桁的には（図２６の丸い囲みのように）
ほぼ同じとみていいでしょう。
この桁、つまり7e-6の各残差(比率）は
図２６では最大値と最小値のちょうど中央あたりになります。
これはどこから来るのか？なぜそうなるのか？
を今までのデータから調べてみました。
図２８がμ=0.9のときのｙ（ｔ）とその推定残差(比率）の関係です。
横軸がｙ（ｔ）でこの最大値は、
右側の縦長の赤い丸印付近にあたります。
また、この位置に相当する縦軸、
つまり横長の赤い丸印付近の値が推定残差(比率）になります、
ということで、この部分を横方向に拡大したものが図２９です。
図２８
図２９のグラフの緑色の丸印の位置がｙ（ｔ）の最大値で、
これに対応する縦軸の値（緑色の丸印）が
残差(比率）になります。
この位置は推定残差(比率）の最大値より
3桁ほど低くなっています。
ということでこの違いが説明できました。
図２９
とここまで書いたところで、
重大な（？）思い違いに気が付きました。
入手した振幅の値（＝理論値）は、
周期解の極大値かつ最大値なのですが、
数値積分した結果の最大値は
0.1秒刻みで計算した中の最大値であり、
必ずしも極大値と一致しているわけではない
ということです。
例えば、
極大値が0.1秒刻みの中間の時間0.05秒にあるかもしれない・・・
極大値こなる保証がまったくないまま
議論していたことになりますね。
なので、
残念ながら
今までの周期解の振幅については真値残差、推定残差の両方とも
意味がないことになります。
現状での計算結果で言えることは、
少なくとも真値残差(比率）の極性（プラスマイナスの符号）は
（最大値が極値に完全に一致しない限り）
必ずマイナスでなければならないはずが、
プラスになっているということくらいです。
これは
おそらく積分の刻み幅を数桁小さくすれば改善されるでしょう。
今までの実績から見て、
PCのマクロ実行時間が１０日以上掛かるとは思いますが・・・
まあ、
推定残差（比率）が三桁とか四桁とか言っているところに
文字通り、桁違いの話をするのもどうかと思うのですが・・・
(^^;
最後の問題、というか課題はやや厄介ですね。
尤も
本来の目的とはずいぶんかけ離れているので
おそらくやらないとは思いますが
もし、やるとすれば
極大値の検出は
現在の方法、つまり
0.1秒の刻み幅表示（＝エクセルのセルへの書き込み）する方法だと
ほぼ不可能です。
基本的なアルゴリズムを変更して
積分処理の最小ループの中に
最大値と最小値の検出処理を入れて
最小の刻み幅単位での検出を毎回実行するように
プログラムを変更する必要がありますね。
おそらくPCのマクロ実行時間が一けた以上増えると思います。
さらに
最大値と最小値の検出は
周期解（＝リミットサイクル）に入ったあとの時間から
行わなければなりません。
初期値を限定すればなんとかなりそうですが、
あまり汎用性のないプログラムにはしたくないので、
この検出方法も考える必要があります。
気が向いて、気分がいいときにやるときが来るか・・・な？？
色々使えるプログラムにするには
想定される問題を事前に明確にしたあとで ←これ大事！
ひとつひとつ解決しないといけないものです。
どこかの国のなんとかカードの問題とよく似ていますね。
