Octave – DECODE

Langlands Program

admin22 — Wed, 17 Feb 2016 15:24:08 +0000

エドワード・フレンケル「数学の大統一に挑む」(Love and Math -The Heart of Hidden Reality-)の中で語られているある数式にとても興味を持ち、プログラムによって確認してみました。
ハリウッド俳優のようなルックスの著者は、数学の教授で、NHKの白熱教室でもこの数式に関してとても興奮して語っていたのがとても印象的でした。
フェルマーの最終定理の証明することとなった、志村・谷山・ヴェイユ予想はラングランズ・ブログラムの特殊ケースといわれ、ある方程式を解くことで、その不思議なつながりを体験することができると解釈しました。
番組ではこの日本人数学者のエピソードも語られ、その偉大さを初めて知りました。

演算に関してだけ簡単に説明すると、

y^2 + y = x^3 – x^2

この方程式と

q * (1-q)^2 * (1-q^11)^2 * (1-q^2)^2 * (1-q^22)^2 * (1-q^3)^2 * (1-q^33)^2 * (1-q^4)^2 * (1-q^44)^2 …..

という規則性のある式は、元々全く関係のないものですが、前者の式の解の個数と後者の式を展開した時の係数が一致するというものです。
プログラムでは、前者を素数を法とする数で総当たりに解を導きます。後者は数式を展開し係数を取得します。
前者は、数が大きくなればなるほど、計算コストがかかるのに、後者はそれを数式の展開という形で導き出してしまう、ということがとても不思議な部分です。
前者はSwiftで、後者はOctaveで計算しました。

Swiftでは、まず素数列を取得しそれから各素数に対して、方程式を解きます。
素数から解の個数を引いた値と、Octaveで演算した素数番目の係数と見比べます。

素数
3,5,7,11,13,17,19,23,29,31,37

解の個数
4,4,9,10,9,19,19,24,29,24,34

(素数)-(解の個数)
-1,1,-2,1,4,-2,0,-1,0,7,3

Octaveでは、

a1 = [1 0];
function ret = insz(s)
  A = [-1];
  for i = 1:s
    A = [A, 0];
  endfor
  A = [A, 1];
  ret = A;
endfunction
for i = 1:37
        r = insz(i-1);
        a1 = conv(a1, r);
        a1 = conv(a1, r);
        r2 = insz(i*11 - 1);
        a1 = conv(a1, r2);
        a1 = conv(a1, r2);
endfor
disp(a1)

このプログラムにより、係数だけの演算をします。insz関数はゼロを配列の真ん中に挿入しています。
[1 1] -> x + 1
[1 0 1] -> x^2 + 1
[1 0 0 1] -> x^3 + 1
といった具合です。

しくみがシンプルなだけにこの結果に驚かされます。
このような現象が、まだまだたくさん隠されているということなのでしようね。
– Hidden Reality –

(数学者って憧れます)

Differential Equations

admin22 — Sun, 20 Jul 2014 16:02:13 +0000

微分方程式の解を求めるOctaveのlsode関数を使って、微分方程式の復習と、任意データへの適用を試みてみました。

function ydot = f1(y, t)
	ydot = -y;
endfunction

t = linspace(0, 10, 100);
y = lsode("f1", 1, t);
plot(t, y)           #1
plot(t, exp(-t))     #1

plot(t, f1(t, 0))    #2

まずシンプルな微分方程式の解をプロットしてみました。
#1の一つ目は、lsodeで解いたもの、二つ目は、式で解いたものをですが、同じものになることを確認しました。

dy/dt = -y
1/y dy = -1 dt
log y = -t + C
y = K exp(-t) (t=0, y=1)

式はこのように解き、#2も一応確認しました。

tについての微分も試してみました。

function ydot = f2(y, t)
	ydot = -t;
endfunction

y = lsode("f2", 0, t);
plot(t, y)                  #3
plot(t, -1 / 2 .* t .* t)   #3

plot(t, f2(0, t)); #2

#2はf1のときと同じになります。

y= – 1/2 t^2
dy/dx = -t

次が今回やりたかった任意の離散データを対象にしてみました。

function ydot = f3(y, t)
	x1=[0 2 3 5 5.5 6 7 8 9 10];
	y1=[1 8 3 2 5   9 4 6 7 0];

	ydot = spline(x1,y1,t)
endfunction

y = lsode("f3", 1, t);
plot(t, y) #4
plot(t, f3(1, t)); #5

splineを使い、連続関数のようにしてテストしてみました。

適当につくったデータですが#5を積分すると#4、なんとなくイメージできます。しかしspline関数、なかなか綺麗に描いてくれますね。

あと、便宜上関数定義を他と同じ個所に書いていますが、Mファイルに同名ファイル(例 f1.m)を作って、別コンソールでエディタを開き編集しながらテストすると、便利です。

Octave Transfer Function

admin22 — Sat, 07 Jun 2014 16:21:29 +0000

Matlabというシミュレーションアプリケーションの巨大パッケージがありますが、それの互換フリーソフトの位置づけにあるのが、Octaveです。

このブログでも最近はデータ処理の観点からRというアプリケーションをよく使用しますが、もともとは10年ほど前に組み込み系の仕事でMatlabを使っていた経緯から、このRがとてもとっつきの良いものになっています。

データビジュアライゼーションについても、最近よくとりあげられていますが、MatlabやMathmaticaなどのソフトでかなり以前から実現されていました。

ここでOctaveがでてきた経緯としては、行列の演算などの処理でもともとMatabでやってきたことを、振り返る意味でとりあげました。Matlabは大変高価なアプリでとても個人では手がでませんので、Octaveを使用します。

久しぶりにOctaveをさわってみましたが、まだ現在進行形でサポートされていることと(2014, version 3.8.1)、伝達関数(tf)がMatlab互換の方法で使えるようになったこと(Octave Control Systems Toolbox : 以前はOctave独自の方法だった気が・・MatlabではSimulinkでやっていた)がとてもうれしいです。
ということで、よく利用した「一次おくれ」のテストをしてみました。

環境 : Ubuntu 14.04

インストール方法 :

apt-get install octave
apt-get install octave-control

テスト : 矩形波を伝達関数を通してPlotします。

コマンド:

[u,t] = gensig('square',4,10,0.1);
H = [tf([1],[1]) ; tf([1],[1 1])]
lsim(H,u,t)

便宜上、グラフを二つ表示するために、一つ目を何もしない伝達関数としました。

うーん。なかなかわかりやすいグラフになりました。

制御系のアプリですが、Octaveを統計データの処理をする上で使うのもありかなと思っています。

このあたり、なかなか興味深いです。

参考 : http://www.mathworks.co.jp/jp/help/ident/ref/lsim.html