kashiの日記

• コンパイル時にマクロKV_USE_TPFMAを定義すると、twoproductの計算にfma命令を使うようになる。
• dka.hpp (Durand Kerner Aberth法) の重解時の安定性が向上。
• 精度保証付き数値計算とkvライブラリの概要を説明するスライドを掲載

• 長い休みはGWかお盆休みしかなくて、西の方はお盆休みは暑すぎて無理なので、GWに行くしかない。
• 以前(11年前)に国道1号で大阪まで自転車で行ったことがあるので、その続き(国道2号)を一度走ってみたかった。

4月29日(土)

4月30日(日)

5月1日(月)

5月2日(火)

5月3日(水)

5月4日(木)

5月5日(金)

(2020年3月31日追記)

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

• 「スタート→歯車アイコン→アプリ→プログラムと機能→Windowsの機能の有効化または無効化」で、「Windows Subsysyem for Linux (Beta)」をチェックする。再起動。
  

  

  
  

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

  

  
  

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

  

  
  

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

  

  
  

lxrun /uninstall /full

• gnuplot
• firefox
• lxterminal
• Xアプリでの日本語入力
• tex

#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);
}

([-1.1313544770895053,-1.131354477089504],[-0.33940634312685164,-0.33940634312685119])
([0.63135447708950442,0.63135447708950499],[0.18940634312685131,0.18940634312685154])

([-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])

([-1.1313544770895053,-1.131354477089504],[-0.33940634312685298,-0.33940634312684986])
([0.6313544770895042,0.63135447708950521],[0.18940634312685095,0.18940634312685187])

([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])

#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++;
	}
}