振動子は同期する

この記事はrogy Advent Calendar 2020 5日目の記事です。

勢いで登録したはいいんですが、全く書きたいことなくてどうしようかな~、というところで友人氏に要望はないか聞いてみたら、有難いことに答えてくれたので、フーリエ級数展開線形代数っぽい視点から眺めてみるやつを書くか~と思ってちょっと調べてみたら、それっぽいちゃんとした記事(http://blog.physips.com/entry/fourier_orthogonal )が一瞬で出てきてしまったので、こちらのほうを参照頂ければと思います(完)。

さて、では今回は何を書くかというと、共通のガウスノイズを受けた振動子が同期していく様子を見てみたいと思います。(なんで?)

振動子が同期するとは

振動子というのは、読んで字のごとく振動するもののことです。メトロノームとかです。で、振動子が同期する、というのは、最初全然違うふるまいをしていた複数の振動子が、時間がたつにつれて同じようなふるまいをするようになることです。

メトロノームが同期する動画を見たことある人もいるんじゃないでしょうか。

www.youtube.com

あと、ホタルのチカチカとか。

www.youtube.com

これらは相互作用によって同期しているのですが、これから紹介するのは、少しだけ毛色の違う同期をします。

共通のガウスノイズを受ける振動子

はい。では、今回考えていく振動子に登場していただきます。じゃん

f:id:conserva:20201205002914g:plain
FitzHugh-南雲モデル

これはFitzHugh-南雲モデルというやつ(のパラメータを適当に設定したやつ)で、ごらんのとおり2次元平面をグルグル回ります。神経細胞の発火活動のモデルに使われたりするらしいです。 で、違う初期値から始まるやつをもう一つ用意します。

f:id:conserva:20201205134513g:plain
2つの振動子たち

2つの振動子たちは今お互いに全然関係ない二人なので、バラバラです。ここに、共通した弱いノイズ(ガウスノイズ)を与えてみます。すると

f:id:conserva:20201205134240g:plain
同じノイズを受けたふたり

gifをループさせているので途中で分裂してるみたいに見えることに気づきましたが、ちゃんと最後同期してます。ちょっとプルプルしながら2つの〇が近づいていく様子がうかがえると思います。

2つの振動子は何ら相互作用をしていないのに、世界から同じ揺れを与えられたら同期してしまったわけです。

何で同期するのか

FitzHugh-南雲モデルに限らず、こういう風に放っておくとある軌道をグルグル回りだすような系はたくさんあります。わっかの上をグルグル回るので、振動子の状態をそのわっかの上の角度みたいなものを使って考えられそうです。その角度のようなものを位相(phase)と呼びます。2つの振動子のうち一つの位相を\theta_1、もう一つの位相を\theta_2とし、\phi = \theta_1 - \theta_2とすると、\phiが0に収束すれば、同期したといえます。詳細については参考文献を後ろに乗せるので良ければ参照していただきたいですが、リアプノフ指数というのを(近似的に)計算することによって、共通のノイズを受ける同じ特性を持つ2つの振動子は同期する傾向を示すということが言えます。

使ったコード(MATLAB R2020b)

a = -0.1;                % FitzHugh-Nagumoのパラメータ
epsilon = 0.008;          % FitzHugh-Nagumoのパラメータ
gamma = 2;               % FitzHugh-Nagumoのパラメータ
varepsilon = 0.003;       % ノイズの強さ

dt = 0.01;
t = 0:dt:1000;
x0 = [0.1;0];           % 振動子1の初期値
y0 = [1;0.05];             % 振動子2の初期値
x = zeros(2,length(t)); % 振動子1の状態
y = x;                  % 振動子2の状態
x(:,1) = x0;
y(:,1) = y0;
z = x;   %軌道プロット用

for k = 2:length(t)
    noise = varepsilon*sqrt(dt)*randn*[1;1];            % [1;1]方向にノイズ
    x(:,k) = fitz(x(:,k-1),a,epsilon,gamma,dt) + noise; % Euler-Maruyama法
    y(:,k) = fitz(y(:,k-1),a,epsilon,gamma,dt) + noise;
    z(:,k) = fitz(z(:,k-1),a,epsilon,gamma,dt);
end
%plot(t,[x;y]); % 状態のプロット

%% plot(animation)
n = 48;
f = @(t)plot([x(1,round(n*t/dt)+1);y(1,round(n*t/dt)+1)]...
    ,[x(2,round(n*t/dt)+1);y(2,round(n*t/dt)+1)],'o');
plot(z(1,round(end/2):end),z(2,round(end/2):end)); % 軌道のわっかを表示
hold on
fanimator(f,'AnimationRange',[0 20])
axis([min(x(1,:))-0.1 max(x(1,:))+0.1 min(x(2,:))-0.1 max(x(2,:))+0.1])
hold off

%% playAnimation
playAnimation

function y = fitz(x,a,epsilon,gamma,dt)
    y = [x(1) + (-x(1)*(x(1)-1)*(x(1)-a)-x(2))*dt;x(2) + epsilon*(x(1)-gamma*x(2))*dt];
end

セクションごとに実行するとfigureが動き出します。

参考文献

ネットワーク・カオス- 非線形ダイナミクス,複雑系と情報ネットワーク - (情報ネットワーク科学シリーズ 第 4巻), 中尾 裕也 (著), 長谷川 幹雄 (著), 合原 一幸 (著), 電子情報通信学会 (監修) , Amazonのリンク

市川雛菜さんに関する妄想

はじめに

今から、先ほどシャワーを浴びている最中に唐突に降ってきた妄想を垂れ流します。

本文

市川雛菜さん宅の隣の家の男の子は、毎朝笑顔で雛菜さんに挨拶してくるが、服や髪の隙間から垣間見える痣や火傷の痕、時おり隣家から聞こえてくる怒鳴り声などが、ご近所の井戸端会議の話題に上ったり上らなかったり。

雛菜さんは特に気にせず暮らしている。


ある日、樋口円香さんが、恐らくは何らかの書類を届けるとか、そういったことのために雛菜さん宅にやってくる。

たまたま件の男の子に遭遇し、やはり笑顔で挨拶される。目に怪我でもしたのか、眼帯をしている。円香さんは運がいいので、風が吹いた拍子とかにたばこによってついたと思しき火傷痕を見つけちゃう。

とりあえずその場では特に気にしないことにしたが、雛菜さんと話してる途中に怒鳴り声とか聞こえてきたので流石に気になって、雛菜さんにそれとなく訊いたら、

「ん~、でもヨソはヨソだからね~」

と答えられて空気が少し凍る。

円香さんはとても優しいけど、なかなか一歩が踏み出せないタイプだった(個人の見解)ので、すごく気にしてしまうものの特に行動は起こせない。

しかし、アイドルとしての日々を通して性格がちょっと変化して、大事な時には一歩を踏み出せるようになっていた(個人の見解)円香さんは、数日後勇気を出して適切な機関に通報する。(すごい!)


そしていろいろ順調に言った結果、お隣のご家族は離れて暮らすようになる。

円香さんはそれを知って、態度には出さないものの心のどこかでは少し満足感を覚えているし、この経験を糧にもっと成長していく。

しかし男の子は決して有難いなんてことは思っていないし、それによって幸せになることもない(なぜなら僕がシャニマス世界に放った刺客なので)。

雛菜さんは男の子が全然幸せになってないことに気付いているけど、やっぱり気にしないで今まで通り日々を送る。

おわりに

今回のオススメ作品は、岡崎京子氏の『pink』です。

ではまた~ノシ

近況

はじめに

ごきげんよう。正直ブログ面倒になってきておりますが、「何か書きませんか」みたいなメールがはてなブログからメチャメチャ来るので近況でも書いておきます。

最近は新型肺炎のあれやこれやで自粛ムード気味(?)ですが、そのせいかは知りませんが僕はいつもより活力に満ち満ちているような心地です。 例年の春休みにはほとんど外出しませんでしたが、今年はお陰様でいろいろとお出かけし、先週には長らくの夢であった人吉(夏目友人帳とまいてつの聖地巡礼)に行きました。 人吉でのあれやこれやについてもいずれまとめて書きたいところではありますが、ものぐさ故旅行してもあまり写真が残っていないので多分寂しい感じになってしまうでしょう。別にいいけど。

けんきう室

実は3月の初めに研究室に所属する予定だったのですが、これも新型肺炎でどこかへ行きました。どこまで行ってしまったのでしょう。例年通りならあと4か月ほどで卒論(?)の発表があるはずですが。

マァ、とにかく唐突に謎のモラトリアムが生まれたので好きなことを好きなだけやってます。唐突な休みがこんなに楽しいとは。 早速図書館からイロイロ借りてきましたが、今はそのすべてを放っておいて、何故かOne of 積読であるところの『線型代数学』(佐武一郎著)を読んでいます(なんで?)。

線型代数

いい本ですね。今まで積読にしておいてナンですが。 お陰で学部一年で通り抜けたはずの線形代数の基礎的な部分が結構抜けていることが分かりました。まだ埋まってないですが...。 いちおー初めてでも分かりやすいように御馴染みの数ベクトルから始まってるんですが、後の方でやる公理から組み立てるベクトル空間にもすぐに馴染めるような説明になっていると思います。ただ、それでも一度どこかで線形代数にちょろっと触っておいた方が分かりやすいかもしれません。

おわりに

ま~他にも細かいことは結構いろいろありましたが、こんなものですかね。

前回「これからも色々アニメとか布教する記事書くよ~」みたいなこと書いた気がするんですが、内容にあまり触れずオススメするのメチャメチャ難しいことに気づいたので止めます。 代わりに記事書いたら最後の方にチョロっと題名だけ書くみたいな感じで行こうかと思います。

というわけで今回の作品は...

『忘念のザムド』です。

ではまた~ノシ

ガラクタ通りのステインはいいぞ

まえがき

ブログ始めてみたはいいものの余りにも書くことがないので、とりあえずこれからは一週間に一度くらいの頻度でアニメやらの布教記事を出して行こうという方針です。マァ布教といっても見て頂かないことには魅力が伝わらないと思うので、見始める動機にでもなれば。

本題

兎にも角にも、記念すべき第一回は「ガラクタ通りのステイン」というアニメとなります。 このアニメは一話7分かそこらですし、調べたら(その法的な正当性は分かりませんが)Youtubeやらニコ動やらに一話が転がっていたので見始めやすいと思います(dアニメストアにもあります)。とりあえず再生して頂きたい。

先ずオープニングが始まります。このオープニングが気に入っていただけたなら、恐らくその後の本編も気に入っていただけると思います。僕は初めてこのオープニングを見たとき、その音楽も大変気に入ったのですが、なによりこの部分でストーンと恋に落ちました。

f:id:conserva:20200209154621p:plain

そして本編に突入し、第一話「卵」が始まります。あとは実際に見てお楽しみください。本当にオススメです。

因みに僕は第四話「極楽鳥」と第八話「リモコンロボット」が特に好きです。

ブログ始めました!!!!!!!!!!!!!!

この記事はrogy Advent Calender 201920日目の記事です。

こんにちは。痕鯖といいます。

これを機にブログをはじめようかと思い立って、実際にはじめました。おわり。

と行きたいところですが、折角なのでMATLABの小ネタを二つくらい披露したいと思います。

歌からメインボーカルだけ取り出す

先ず曲を読み込みます。 onvocalにはメインボーカルが入ったやつ、offvocalには入ってないやつを入れます。

[onvocal,fs] = audioread('01 秘密のトワレ.flac');
offvocal = audioread('03 秘密のトワレ (オリジナル・カラオケ).flac');

下の関数に入れます。

function y = vocalcut(onvocal,offvocal,fs)

t = 2*fs; %最初2秒分だけ見る
a = onvocal(1:t,1);
b = offvocal(1:t,1);
[c,d] = xcorr(a,b);
[~,ind] = max(abs(c));
lag = d(ind)
if lag>0 %ズレの正負でonvocalとoffvocalのどっちをずらすか
    y1 = onvocal(lag+1:end,:);
    if length(y1)>length(offvocal)
        y1 = y1(1:length(offvocal),:);
        y2 = offvocal;
    else
        y2 = offvocal(1:length(y1),:);
    end
else
    y2 = offvocal(abs(lag)+1:end,:);
    if length(y2)>length(onvocal)
        y2 = y2(1:length(onvocal),:);
        y1 = onvocal;
    else
        y1 = onvocal(1:length(y2),:);
    end
end
y = y1-y2;

end

ボーカルだけ取り出すのにvocalcutってどうなの?と思ったときには作り終わってたので気に入らなかったら変えてください。

出てきたやつを再生します。

sound(y,fs); % clear soundでとまります。

やってることとしては、onvocalとoffvocalのイントロの部分(ここでは簡単のため最初の2秒分だけ)を比較してズレを計算して、onvocalからoffvocalを引けば多分ボーカルだけ取り出せるだろうという感じです。

あまり考えていないので、例えば最初アカペラから始まる曲とか多分失敗します。

多項式の根をグレブナー基底を使って調べてみる

グレブナー基底って何?みたいな詳しいことは書きません(書けません)。何やってるか知りたいよ、という方には『グレブナ基底と代数多様体入門』という本をおすすめしておきます。

今回は解きたい問題として「辺の長さの合計が2の三角形の中で一番面積が大きいのは何?」というのを考えます。

とりあえずコードをはります。a, b, cは各辺の長さでlambdaはラグランジュの未定乗数です。

ここではただ多項式の根を調べたいだけなので、何かしらわからなくてもとりあえず「こーゆう風に多項式の根が分かるんだ~」という感じでお願いします。

syms a b c lambda real % 変数たち
f = (1-a)*(1-b)*(1-c); % 三角形の面積(ヘロンの公式)
g = a+b+c-2; % 辺の長さは合計で2という制約(g=0)
h = [gradient(f)-lambda*gradient(g); g] % 解きたい多項式(全部"左辺=0"のカタチの左辺)
G = gbasis(h.', 'MonomialOrder', 'lex') % (簡約された)グレブナー基底
kotae = solve(G); % 上の行で出てきた多項式=0を解く(実はsolve(h)で解ける)
disp(kotae.a);disp(kotae.b);disp(kotae.c);disp(kotae.lambda); % コタエを表示

大事なのは5行目で、Gは多分hよりも簡単になったやつ(一番下はlambdaだけの式)になってるはずです。

おわりに

ここまで読んでいただきありがとうございました。これからも日記的なものをチョビチョビ更新していけたらな~と思っています。