2017/05/07(日)だいたい国道2号自転車旅

たまには遊びの日記も書いてみよう。

このGW、9連休だったので、長い休みでないと出来ないことをしたいと思って、大阪発で西へ西へ自転車で向かってみました。
  • 長い休みはGWかお盆休みしかなくて、西の方はお盆休みは暑すぎて無理なので、GWに行くしかない。
  • 以前(11年前)に国道1号で大阪まで自転車で行ったことがあるので、その続き(国道2号)を一度走ってみたかった。
あたりが一応の動機でしょうか。自転車に乗るのは昨年8月のお盆以来なので、タイヤの空気を入れるところから始まります。大丈夫か俺。

4月29日(土)

夕方荷造りして自転車で東京駅へ。新幹線に載せて新大阪へ。夜10時頃新大阪着。

25.97km。

4月30日(日)

新大阪から少し迷いつつ、梅田の国道1号から2号に変わる交差点へ。ここから2号を西へ向かう。

IMG_5227.JPG
IMG_5229.JPG
IMG_5231.JPG
IMG_5232.JPG


体は重いが、だらだら走ってやがて神戸に。震災慰霊碑など。

IMG_5235.JPG
IMG_5236.JPG
IMG_5237.JPG
IMG_5238.JPG
IMG_5239.JPG
IMG_5240.JPG


明石着。明石と言えば東経135度。マンホールにも135って書いてあった。

IMG_5243.JPG
IMG_5247.JPG
IMG_5249.JPG


昼を食べながらこの日のゴールを考えた。ゴールは姫路、相生、赤穂温泉あたりが考えられたが、姫路では少し近すぎと思い、赤穂温泉に宿を予約した。が、どうもgooglemapの操作をミスっていたらしく、赤穂温泉は予想よりも20kmほど遠いことに走りながら気づく。こういう旅のときは基本的に景色を楽しみながら走りたいので夜間走行はしないポリシー。しかし、少し頑張らないと日没前に赤穂温泉に着くか怪しい。

時間が無いので姫路城は1分でスルー。

IMG_5253.JPG


赤穂市は国道2号を離れ国道250号へ。地形は山がちになり、標高100m程度のプチ峠を越え、何とか日没前に到着。

IMG_5256.JPG
IMG_5259.JPG


この日の走行距離は135.74km。体が慣れるまでは抑え気味で行こうと思っていたのに、少し走り過ぎた感。宿が思いの外いい宿で、温泉サイコーで、このままここに滞在したい気分に。

5月1日(月)

いい宿に後ろ髪をひかれつつ、出発。とりあえず市内を少し回った後、先へ向かう。

IMG_5266.JPG
IMG_5267.JPG
IMG_5268.JPG
IMG_5269.JPG
IMG_5270.JPG
IMG_5271.JPG
IMG_5272.JPG


軽く峠を越えて、岡山県に入る。

IMG_5273.JPG
IMG_5274.JPG


岡山市街地で見かけた。点字ブロック発祥の地だそうだ。

IMG_5286.JPG
IMG_5287.JPG
IMG_5288.JPG
IMG_5289.JPG
IMG_5290.JPG
IMG_5291.JPG


この日のゴール、倉敷では近すぎ、福山では少し遠いと思ったが、頑張って福山まで行くことにした。

国道2号の200kmポスト。

IMG_5295.JPG


福山市からは広島県。岡山県は一日で走り抜けてしまった。

IMG_5297.JPG


福山着。駅前のビジホの窓から福山城がとてもよく見えたのには驚いた。

IMG_5298.JPG
IMG_5300.JPG
IMG_5302.JPG
IMG_5303.JPG


この日の走行距離は122.45km。またも少し長め。この日の宿もとても快適だった。

5月2日(火)

朝もきれいな福山城を窓から見て出発。

IMG_5305.JPG


すぐに尾道市に入る。しまなみ海道は以前走ったことがあるのでスルー。

IMG_5309.JPG
IMG_5311.JPG
IMG_5312.JPG
IMG_5323.JPG
IMG_5324.JPG


原田知世の「時をかける少女」の撮影に使われた学校らしい。

IMG_5316.JPG
IMG_5317.JPG
IMG_5318.JPG
IMG_5319.JPG
IMG_5320.JPG
IMG_5321.JPG


出発して以来ずっと、国道2号走りづらいなあと感じていたところで、こんな看板を見て海沿いを行った方が楽しいだろうと予想、国道185号で呉方面に向かうことにする。

IMG_5329.JPG


国道185号は予想通りとても走りやすく、景色も最高。

IMG_5330.JPG
IMG_5331.JPG
IMG_5332.JPG


竹原市に入る。竹原には江戸時代の町並みを保存した地区があり、「時をかける少女」の大半は尾道ではなく竹原で撮影されたらしい。確かにこんな感じだったような気がする。

IMG_5347.JPG
IMG_5348.JPG
IMG_5349.JPG
IMG_5350.JPG
IMG_5351.JPG
IMG_5352.JPG
IMG_5353.JPG
IMG_5354.JPG
IMG_5355.JPG
IMG_5356.JPG
IMG_5357.JPG
IMG_5358.JPG
IMG_5359.JPG
IMG_5360.JPG
IMG_5362.JPG
IMG_5363.JPG
IMG_5364.JPG
IMG_5365.JPG
IMG_5366.JPG
IMG_5367.JPG
IMG_5368.JPG
IMG_5369.JPG
IMG_5370.JPG
IMG_5371.JPG
IMG_5372.JPG
IMG_5373.JPG
IMG_5374.JPG


最後は1700mもの長いトンネルを通って、呉に到着。

IMG_5379.JPG
IMG_5381.JPG
IMG_5383.JPG


そして呉と言えば、「この世界の片隅に」の舞台。駅ビルの中を探索してみたところ、あちこちで押してるのが分かる。

IMG_5385.JPG
IMG_5386.JPG
IMG_5393.JPG
IMG_5397.JPG
IMG_5399.JPG
IMG_5411.JPG
IMG_5412.JPG
IMG_5413.JPG
IMG_5414.JPG


この日の走行距離は109.70km。割と観光モードでした。脳内BGMはずっと原田知世の「時かけ」。

5月3日(水)

さて、1日2日が休みでなかった人にとってはこの日から本格的にGW。で、検索してみると、先へ進もうにも山口県内の宿が全く見つからない。仕方がないので今日は広島市街地までとし、疲労も溜まっているので休養も兼ねて観光モードに。

まずは「この世界の片隅に」聖地巡礼。

灰ヶ峰は街のどこからでも見える。かつてここに高射砲があった。

IMG_5445.JPG


小春橋。ただの橋である。

IMG_5446.JPG
IMG_5447.JPG
IMG_5448.JPG
IMG_5449.JPG


このへんも映画に出てきた気がする。

IMG_5454.JPG
IMG_5455.JPG
IMG_5456.JPG
IMG_5457.JPG
IMG_5458.JPG


歴史の見える丘。港を見下ろせ、戦艦大和を建造したというドックも見える。

IMG_5463.JPG
IMG_5464.JPG
IMG_5465.JPG
IMG_5466.JPG
IMG_5467.JPG
IMG_5468.JPG
IMG_5469.JPG
IMG_5474.JPG
IMG_5475.JPG


その歴史の見える丘の近く、映画中で最も悲しい場面はここで起きたという設定らしい。

IMG_5461.JPG
IMG_5476.JPG
IMG_5477.JPG


大和ミュージアムと海上自衛隊史料館。

IMG_5480.JPG
IMG_5481.JPG
IMG_5482.JPG
IMG_5484.JPG
IMG_5485.JPG


大和の甲板の1/4と同じ大きさの広場らしい。

IMG_5486.JPG
IMG_5489.JPG
IMG_5490.JPG
IMG_5492.JPG
IMG_5493.JPG
IMG_5495.JPG


映画に出てきた三ッ蔵。

IMG_5501.JPG
IMG_5502.JPG
IMG_5503.JPG
IMG_5504.JPG
IMG_5505.JPG


山の上に登って行って、主人公の家の最寄りのバス停「辰川」を目指す。

IMG_5506.JPG
IMG_5507.JPG
IMG_5508.JPG
IMG_5509.JPG
IMG_5510.JPG
IMG_5511.JPG
IMG_5512.JPG


今はバス停から港を見下ろすことはできない。

IMG_5515.JPG


そこで、更に上の住宅地を探索して、港を見下ろせる場所を探してみた。

IMG_5514.JPG
IMG_5523.JPG
IMG_5524.JPG
IMG_5525.JPG
IMG_5527.JPG


さて、呉観光はこのくらいにして、広島市街地へ向かうことにした。30kmほど。

IMG_5532.JPG
IMG_5533.JPG
IMG_5534.JPG
IMG_5535.JPG


広島城。

IMG_5536.JPG
IMG_5538.JPG


原爆ドーム。

IMG_5542.JPG
IMG_5543.JPG
IMG_5544.JPG
IMG_5550.JPG


相生橋。

IMG_5545.JPG
IMG_5546.JPG
IMG_5547.JPG
IMG_5548.JPG
IMG_5549.JPG


平和記念公園。

IMG_5551.JPG
IMG_5552.JPG
IMG_5553.JPG
IMG_5554.JPG
IMG_5555.JPG
IMG_5556.JPG


カープ戦のある日だけ提供されるという超絶カロリー食、ミートソースカツスパゲッティを食べた。

IMG_5564.JPG
IMG_5566.JPG


この日の走行距離は54.99km。完全観光モードの日でした。

5月4日(木)

前日宿が無くて広島停滞を余儀なくされたが、この日は検索してみると宇部に宿を発見。しかし、ここから宇部までは180kmもある。ここまでの実績を考えると日没前に180kmは厳しそうだが、仕方が無いので気持ちを切り替えて頑張ってみることにする。まとまった食事はとらず、2時間に一度計画的にコンビニでパンとコーラでカロリーを補給する作戦で。

岩国で山口県に入る。有名な錦帯橋を見る余裕なし。

IMG_5570.JPG


かなり距離は長くなるが、国道2号より快適に走れるであろう海沿いの国道188号を行くことにする。大変すばらしい道で、テンション上げて快走できた。昔よく走っていた頃の感覚が少し蘇ってきた感じ。野生解放(笑)。

IMG_5571.JPG
IMG_5573.JPG
IMG_5574.JPG
IMG_5575.JPG
IMG_5576.JPG
IMG_5577.JPG
IMG_5578.JPG


関東では全くみないけど西日本ではやたらとみる印象的な看板。

IMG_5582.JPG
IMG_5583.JPG
IMG_5584.JPG


いよいよゴールの北九州まで100kmを切った。

IMG_5588.JPG


国道190号へ分岐し、何とか宇部に到着。

IMG_5593.JPG
IMG_5594.JPG


この日の走行距離は実に189.20km。200kmオーバーは過去に何度かあるが、みんな夜間走行を含んでいるので、多分日没前にこれだけ走ったのは初めて。サイコンのaveも22km/hを超え、最近の中ではそれなりに速かった。走り出しが7:30、到着が18:30なのでちょうど11時間。いわゆるブルベの200kmは13時間半が制限時間と聞いているので、このくらいのペースなら楽に完走できる、ということか。途中のコンビニ休憩は10:00, 12:00, 14:00, 16:00と2時間おきだった。

5月5日(金)

最終日は距離が短いのでとても気楽。どんどん関門橋が近づいてくる。

IMG_5599.JPG
IMG_5601.JPG
IMG_5602.JPG
IMG_5603.JPG
IMG_5604.JPG
IMG_5605.JPG
IMG_5606.JPG
IMG_5608.JPG


自転車は橋は通れないので、人道トンネルで九州に渡る。

IMG_5609.JPG
IMG_5611.JPG
IMG_5612.JPG
IMG_5615.JPG
IMG_5618.JPG
IMG_5619.JPG
IMG_5620.JPG
IMG_5621.JPG
IMG_5622.JPG
IMG_5623.JPG
IMG_5625.JPG
IMG_5626.JPG
IMG_5627.JPG


国道2号終点の交差点へ。いつかこの続きの国道3号も走るのだろうか。

IMG_5631.JPG
IMG_5632.JPG
IMG_5633.JPG


小倉駅で自転車を畳んで新幹線に載せ、帰宅。

IMG_5639.JPG


ここまでの走行距離は57.54km。


東京駅から自宅まで自走。26.77km。

以上、天気に恵まれ(全く雨に遭わなかった)、充実した旅でした。走行距離は、新大阪-小倉間で669.62km、自宅と東京駅の間の往復も合計すると、722.36km。走ってるうちに少しずつ体力が戻ってくる感じは楽しかったけど、どうやって現実に帰ろうかのう。

2017/04/17(月)bash on Ubuntu on Windowsで遊んでみる (Creators Update編)

bash on ubuntu on windowsで遊んでみるの記事で、bash on Ubuntu on Windowsの紹介をしましたが、それがこの4月のWindows 10 Creators Updateで大きく改善したとのことで、再度入れて遊んでみました。

前回(Anniversary Update時)との差分は、
  • ubuntuのバージョンが14.04から16.04になった。
  • ターミナルが改善し、日本語が文字化けせずきちんと扱えるようになった。WindowsのIMEでの日本語入力も可能に。
  • ubuntuからwindows内のコマンドを呼び出したり、windowsからubuntu内のコマンドを呼び出したりできるようになった。
  • ifconfigやpingなどのネットワーク系のコマンドが動作するようになった。
あたりのようです。

条件はWindows10の64bitで、Creators Updateになっていることです。「スタート→歯車アイコン→システム→バージョン情報」で
version1703.png

のようにバージョンが1703になっていればCreators Updateになっています。まだだった場合でも強制的にアップデートすることもできます。

インストール方法は前回とほとんど変わっていませんが一応書いておきます。
  • 「設定→アプリ→プログラムと機能→Windowsの機能の有効化または無効化」で、「Windows Subsysyem for Linux (Beta)」をチェックする。再起動。
    check-wsl.png

  • 「スタート→歯車アイコン→更新とセキュリティ→開発者向け」で、「開発者モードを使う」を、「サイドロードアプリ」から「開発者モード」に変更。再起動。
    devel-mode.png

  • 「スタートを右クリック→Windows PowerShell(管理者)」でPower Shellを起動し、bashとタイプ。"y"でダウンロードとインストールが始まります。数分かかります。ユーザIDとパスワードを聞かれてインストール完了。ここで入れたパスワードは"sudo"実行時に使います。
    install.png

  • 次回以降は、スタートメニューに「Bash on Ubuntu on Windows」が登録されているのでそこから起動できます。
    startmenu.png

なお、前回のAnniversary Update時にbash on Ubuntu on Windowsを入れていた場合は、自動で新しいものにアップデートされることはない模様です。必要なファイルがあればバックアップした上で、いったんPower Shell(管理者)で
lxrun /uninstall /full
で全削除し、「bashとタイプ」のところからやりなおすのが確実です。
uninstall.png

使い方は、前回のAnniversary Updateのときとほとんど変わっていません。日本語が文字化けしなくなり、またターミナルからwindowsのIMEを使って日本語入力ができるようになったので、ちょっとしたプログラミングの演習程度ならこのままで十分に使えるようになりました。

また、X serverをwindowsに入れてGUIを使うソフトウェアを動かすこともできます。bash on ubuntu on windowsで遊んでみるの記事で紹介した
  • gnuplot
  • firefox
  • lxterminal
  • Xアプリでの日本語入力
  • tex
など試してみましたが全て同じように動作しました。安定度はずっと高まっているように感じました。

2017/03/21(火)Henon mapで遊んでみた

先週宮古島でINVA2017という研究集会に行ってました。有名なTucker先生をお呼びして講演をしていただいたのですが、Hénon(エノン)mapに関する講演で、面白そうな話だったので講演中いろいろ計算して遊んでいました。

Hénon mapは、
\begin{pmatrix} x_{n+1} \\ y_{n+1} \end{pmatrix} = \begin{pmatrix} 1 + y_n - a x_n^2 \\ b x_n \end{pmatrix}
みたいな漸化式で定義される2次元離散力学系です。a,bの値によってはカオスになることで有名な力学系で、a=1.4、b=0.3がカオスになる有名な値です。その値を使って(x,y)=(0,0)から10000点計算してxy平面にプロットしてみると、

henong.png


のような図形が得られました。

さて、適当な点から出発してn回反復したとき元の点に戻るような軌道をn周期解と言います。周期解には、一度その軌道に入ったらずっとそこから出ない安定な周期解と、わずかな摂動で軌道から外れてしまう不安定な周期解があります。一般にカオスになるときは不安定な周期解しか無いのが普通ですが、よく使われるa=1.4, b=0.3でなく、そこからわずかに値を変化させると安定な周期解が現れることが知られていて、そのような安定な周期解をどうやって探索するか、が講演の内容でした(多分)。

で、とりあえずkvライブラリを使って周期解の探索をしてみよう、というのが以下の内容です。とりあえずa=1.4, b=0.3で、(x,y)の探索範囲は、この論文にあった式
r = \frac{1+|b|+\sqrt{(1+|b|)^2+4a}}{2}
を使って(x,y)\in([-r,r],[-r,r])で探索することにしました。

n回反復する写像の不動点を全解探索するプログラムを素直に書くと、
#include <kv/allsol.hpp>

namespace ub = boost::numeric::ublas;

template <class TT>
struct Henon {
	TT a, b;
	Henon(TT a, TT b) : a(a), b(b) {}
	template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
		ub::vector<T> r(2);

		r(0) = 1 + x(1) - a * x(0) * x(0);
		r(1) = b * x(0);

		return r;
	}
};

template <class F>
struct Repeat_n {
	F f;
	int n;
	Repeat_n(F f, int n): f(f), n(n) {}
	template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
		int i;
		ub::vector<T> tmp;

		tmp = x;
		for (i=0; i<n; i++) tmp = f(tmp);
		return tmp;
	}
};

template <class F>
struct FixedPoint {
	F f;
	FixedPoint(F f): f(f) {}
	template <class T> T operator() (const T& x){
		return f(x) - x;
	}
};

typedef kv::interval<double> itv;

int main()
{
	int i, n;
	ub::vector<itv> x;
	itv a, b;
	itv::base_type r;

	std::cout.precision(17);

	a = "1.4";
	b = "0.3";
	n = 4;

	Henon<itv> f(a, b);
	Repeat_n< Henon<itv> > g(f, n);
	FixedPoint< Repeat_n< Henon<itv> > > h(g);

	r = ((1+abs(b)+sqrt(pow(1+abs(b), 2)+4*a))/2).upper();

	x.resize(2);
	x(0) = itv(-r, r);
	x(1) = itv(-r, r);

	kv::allsol(h, x);
}
のようになりました(henon0.zip)。関数オブジェクトを使っていて、fがHénon写像、gがそれをn回繰り返したもの、hが不動点形式(g(x)-x)です。このプログラムをnを変えながら実行してみます。n=1では、
([-1.1313544770895053,-1.131354477089504],[-0.33940634312685164,-0.33940634312685119])
([0.63135447708950442,0.63135447708950499],[0.18940634312685131,0.18940634312685154])
の2つの1周期解(要するに不動点)が得られました。n=2にすると、
([-0.47580005117505659,-0.47580005117505597],[0.29274001535251664,0.29274001535251721])
([0.97580005117505597,0.97580005117505653],[-0.14274001535251713,-0.14274001535251665])
([0.63135447708950431,0.6313544770895051],[0.18940634312685117,0.18940634312685171])
([-1.1313544770895053,-1.131354477089504],[-0.33940634312685198,-0.33940634312685086])
の4つの解が得られましたが、よく見るとそのうち2つは先程得られた不動点です。更に、残りの2つの解も実質同じ軌道(p→q→pが2周期解ならq→p→qも当然2周期解)なので、2周期解は実質1つだけということが分かります。n=3だと、
([-1.1313544770895053,-1.131354477089504],[-0.33940634312685298,-0.33940634312684986])
([0.6313544770895042,0.63135447708950521],[0.18940634312685095,0.18940634312685187])
が得られますが、これは不動点なので、3周期解は存在しないことが分かりました。n=4の結果を示すと、
([0.63135447708950431,0.6313544770895051],[0.18940634312685045,0.1894063431268524])
([0.63819399262715537,0.63819399262715638],[-0.21203003316582356,-0.2120300331658213])
([-0.70676677721940862,-0.70676677721940761],[0.33752098096070709,0.33752098096070832])
([0.2177617657186289,0.21776176571863363],[0.19145819778814535,0.19145819778814813])
([-0.47580005117505797,-0.47580005117505475],[0.2927400153525157,0.29274001535251815])
([-1.1313544770895064,-1.1313544770895032],[-0.33940634312685631,-0.33940634312684653])
([0.97580005117505508,0.97580005117505731],[-0.14274001535251838,-0.1427400153525154])
([1.1250699365356913,1.1250699365356934],[0.065328529715587863,0.065328529715590903])
と8つの解が得られますが、2つは不動点、2つは2周期解で、4周期解は位相をずらした4つが得られるので実質1つの4周期解があることが分かります。

さて、このように人間の目で選り分けるのは大変なので、より短い周期解の繰り返しや位相がずれただけのものを排除する部分を付け加えて見ました。ついでに、その周期解の安定性をその写像の解のところでのヤコビ行列を計算しその固有値の絶対値の最大値を調べることによって判定する部分も付けてみました。
#include <kv/allsol.hpp>
#include <kv/eig.hpp>

namespace ub = boost::numeric::ublas;

template <class TT>
struct Henon {
	TT a, b;
	Henon(TT a, TT b) : a(a), b(b) {}
	template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
		ub::vector<T> r(2);

		r(0) = 1 + x(1) - a * x(0) * x(0);
		r(1) = b * x(0);

		return r;
	}
};

template <class F>
struct Repeat_n {
	F f;
	int n;
	Repeat_n(F f, int n): f(f), n(n) {}
	template <class T> ub::vector<T> operator() (const ub::vector<T>& x){
		int i;
		ub::vector<T> tmp;

		tmp = x;
		for (i=0; i<n; i++) tmp = f(tmp);
		return tmp;
	}
};

template <class F>
struct FixedPoint {
	F f;
	FixedPoint(F f): f(f) {}
	template <class T> T operator() (const T& x){
		return f(x) - x;
	}
};

typedef kv::interval<double> itv;

int main()
{
	int i, n;
	std::list< ub::vector<itv> > result, result2;
	std::list< ub::vector<itv> >::iterator p, p2;
	ub::vector<itv> x, x2;
	bool flag;
	itv a, b;
	itv::base_type r;

	std::cout.precision(17);

	a = "1.4";
	b = "0.3";
	n = 4;

	Henon<itv> f(a, b);
	Repeat_n< Henon<itv> > g(f, n);
	FixedPoint< Repeat_n< Henon<itv> > > h(g);

	r = ((1+abs(b)+sqrt(pow(1+abs(b), 2)+4*a))/2).upper();

	x.resize(2);
	x(0) = itv(-r, r);
	x(1) = itv(-r, r);

	result = kv::allsol(h, x);

	p = result.begin();
	while (p != result.end()) {
		x = *p;
		flag = true;
		x2 = x;
		for (i=0; i<n-1; i++) {
			x2 = f(x2);
			// check the solution is truly n-pediodic or not
			if (overlap(x, x2)) {
				flag = false;
				break;
			}
			// check the solution is new or not
			p2 = result2.begin();
			while (p2 != result2.end()) {
				if (overlap(*p2, x2)) {
					flag = false;
					break;
				}
				p2++;
			}
			if (flag == false){
				break;
			}
		}

		if (flag) {
			result2.push_back(x);
			std::cout << "true pediodic soluion " << result2.size() << "\n";
			x2 = x;
			for (i=0; i<n; i++) {
				std::cout << x2 << "\n";
				x2 = f(x2);
			}
			std::cout << "\n";

			// calculate maximum eigenvalue of Jacobian
			ub::vector<itv> v1;
			ub::vector< kv::complex<itv> > l;
			itv lm, tmp;
			ub::matrix<itv> m1;
			kv::autodif<itv>::split(g(kv::autodif<itv>::init(x)), v1, m1);
			kv::veig(m1, l);
			std::cout << l << "\n";
			lm = 0;
			for (i=0; i<2; i++) {
				tmp = abs(l(i));
				lm.lower() = std::max(lm.lower(), tmp.lower());
				lm.upper() = std::max(lm.upper(), tmp.upper());
			}
			std::cout << lm << "\n";
			if (lm < 1) {
				std::cout << "stable\n";
			} else if (lm > 1) {
				std::cout << "unstable\n";
			} else {
				std::cout << "stability unknown\n";
			}
		}
		p++;
	}
}
(henon1.zip) これで調べた周期解とその安定性は次のとおりです。
nxy安定性
1[-1.1313544770895053,-1.131354477089504][-0.33940634312685164,-0.33940634312685119]unstable
1[0.63135447708950442,0.63135447708950499][0.18940634312685131,0.18940634312685154]unstable
2[-0.47580005117505659,-0.47580005117505597][0.29274001535251664,0.29274001535251721]unstable
4[0.63819399262715537,0.63819399262715638][-0.21203003316582356,-0.2120300331658213]unstable
6[0.44190995135922078,0.44190995135922451][-0.24126597029523553,-0.24126597029522911]unstable
6[1.0380595354868291,1.0380595354868339][0.093435694641492983,0.093435694641505183]unstable
7[-1.0872860458490461,-1.0872860458490447][0.36967525887581315,0.36967525887581654]unstable
7[-1.0466775735267937,-1.0466775735267894][0.35496793781757529,0.35496793781758157]unstable
7[-0.92617327728451871,-0.92617327728451703][0.34236913808086111,0.3423691380808675]unstable
7[0.81803519601224594,0.81803519601224873][0.15474440787832016,0.15474440787832872]unstable
8[-1.1493133062560075,-1.1493133062560048][0.36560087176430355,0.36560087176431478]unstable
8[-0.8274542278695698,-0.82745422786956712][-0.37646388947765531,-0.37646388947764258]unstable
8[0.90002978920450349,0.90002978920451305][0.13188647717963458,0.13188647717967495]unstable
8[-0.44837319092732131,-0.44837319092731264][-0.34223569372342872,-0.34223569372341905]unstable
8[-0.83680827077054099,-0.83680827077053698][0.35013340845532586,0.35013340845533498]unstable
8[-0.80876720316291007,-0.80876720316289762][0.32983594103956154,0.32983594103958375]unstable
8[0.60111180966072053,0.6011118096607403][-0.21508292460527792,-0.21508292460523506]unstable
n=13まで計算したので、解の数のみ示します。全て不安定解でした。
n周期解の数
12
21
30
41
50
62
74
87
96
1010
1114
1219
1332
全解探索に成功したので、各nに対してこれらで全ての周期解を尽くしていることの数学的証明になっているはずです。

さて、この論文によると、a=1.3866414978735625,b=0.3だと、8周期の安定解が存在するのだそうです。実際、aの値をそれに変えて計算してみると、
nxy安定性
8[0.8742591352635014,0.87425913526350741][0.14575249177868274,0.14575249177870631]unstable
8[-0.84160148826306525,-0.84160148826306058][0.35209093466818297,0.35209093466819364]unstable
8[-0.4853729323948971,-0.48537293239488293][-0.34686239645463785,-0.34686239645462463]unstable
8[0.60580464136231371,0.6058046413623337][-0.21568623282216013,-0.21568623282211638]unstable
8[1.1065521562926103,1.1065521562926253][-0.11533020216315528,-0.11533020216311613]unstable
8[0.42199129139879743,0.42199129139889974][0.1755174868908467,0.17551748689086294]stable
8[0.93092246076173934,0.93092246076183872][0.12609408826805895,0.12609408826806834]unstable
と、6つの不安定解と1つの安定解が得られました。8回繰り返す写像の安定解の場所でのヤコビ行列は、
\begin{pmatrix} [-0.2071075038203655,-0.20710750356848622] & [1.3046842272497074,1.3046842274628959]) \\ [0.046141607262891293,0.046141607302409162] & [-0.29098818813819278,-0.29098818810610671] \end{pmatrix}
で、その固有値の絶対値の最大値は[0.49796393435494867,0.49796393621891983]でした。1より小さいのでこの解は安定なことが分かります。

さて、もう少し遊んでみたりもしたのですが記事も長くなってきたのでこのへんまでにしましょう。力学系方面の知識が無いのでこういう軌道を計算することがどう面白いのかはあまり分からないのですが、いろいろ楽しめました。力学系方面で精度保証付き数値計算をしたいという需要はありそうな気がするので、共同研究なんかに繋がるといいなあ。

2017/03/10(金)kv-0.4.41

kvライブラリを0.4.41にアップデートしました。

今回の変更は、キャストし忘れでode solverの内部使用型をddとしたときにコンパイル出来なくなっていた修正だけです。

また、ドキュメントにkvライブラリの応用例のページを追加しました。まだ2つだけしか載せていませんが、手元には載せたい応用例がたくさんあるので今後増やしていきたいと思います。また、このページでは数式を表示するプラグインのMathJaxを使ってみました。何か不具合があればお知らせください。

2017/01/27(金)kv-0.4.40

kvライブラリを0.4.40にアップデートしました。

今回は、T大のT先生の依頼によるODE solver仕様変更が含まれています。その他、psaのドキュメントを最新にしたり、kraw_approxの1変数版を追加したりしました。

ODE solverの仕様変更は以下のようなものです。常微分方程式の計算を行なうodelong_maffineなどのodelong系の関数では、計算しながら何かをさせるためのcallback関数を指定する機能があります (正確に言うと、()演算子を持つcallback用のオブジェクト)。これの戻り値は今まではvoidだったのですが、それをboolに変更してその値をチェックするようにし、もしfalseなら計算を中止するようにしました。これにより、終了条件として時刻でなく、「値が何かの条件を満たしたら終了」のようなことが出来るようになりました。これは便利になっているのですが、もし今までのversionのkvを利用していてcallback関数を継承によりカスタマイズする機能を使っていた場合、コンパイルエラーになってしまうという申し訳ないことになってしまいました。修正は簡単で、operator()を定義する関数の戻り値をboolに変えて末尾にreturn true;を補うだけです。

こういうことはなるべく起きないように仕様を考えるときは気をつけていたのですが、今回はメリットが仕様変更のデメリットを上回ると考えてこうすることに決心しました。恐らくこの機能を実際に使っているのはT大のT先生だけだと思いますが、もし他に使っている方がいらっしゃれば大変申し訳ありません。

また、githubの方のソースコードも最新になっています。結局ローカルの開発体制をgitに移行するには至らず、version up時にまとめてgithubを更新するような何とも無意味な状態になっています。新しいものに馴染むのは難しい…
OK キャンセル 確認 その他