ちょっと興味深い行列式を見つけたので

この行列式は、nを整数とすると、f(n)のようになるようです。

|f(n+1)|/|f(n)|=n!
と
n≡2,n≡3　(mod　4)のときf(n)<0
n≡0,n≡1　(mod　4)のときf(n)>0
であることがわかりました。


まだ導出はできませんでした。できる見通しがあるかどうかもわかりません
途中まで証明しようとして、行き止まりに詰まった感じがしました

偶奇置換とレヴィ・チヴィタ記号

3個の要素を持つ列v1に置換操作σを施してv2にする演算を

v2=σ×v1

としたとき

ε(v2)=sgn(σ)ε(v1)

となるのは、要素数が少なかったからしらみつぶしが楽だったが

n個の要素一般に言えるのかどうかは確かめようがあるのか？

と思っていた。

しかし、4つ以上の足を持つレヴィ・チヴィタ記号(エディントンのイプシロン)の定義が
要素3つのときとは趣がちょっと違っていて
定義からモロに偶置換と奇置換(sgn(σ))だったので

任意の個数nについても同様に

ε(v2)=sgn(σ)ε(v1)

が言えそうだ。



というのも置換操作の例題を見て
列の左から作用させるのに、僕は勝手に右から作用させてしまって
例題の時点で答えと違っていて先が思いやられると一瞬嘆いたりして戸惑ったのだったが

3つの要素の置換を見て左右どちらかというのに気づいたきっかけが
置換したあとの列が、左から作用させても右から作用させてもイプシロンが不変というところからだった。


こうなったらただでは起きないぞと思い
行列の行列式やベクトルのノルムのように
置換に対しても何らかの固有の(数の集合体を代表する単一の)数があるのではないかと思って、

εとsgn(σ)の関係に相当するものがないかと考えてみた次第だった。






＝＝＝＝＝＝＝＝＝＝
添え字が3つのイプシロンは順方向とか逆方向とかで定義できたが
添え字が4つ以上となるとそうはいかず
別のアプローチが必要となる。

たとえば
ε1234=1のとき
ε1234の1ペアだけを入れ替えて、
ε2134とかε1243とかにした場合はε=-1となる。

これは要素が3つだけの
ε123=1のときも
ε213=-1
について、「逆方向」のアプローチのほかにも
「2と1だけを入れ替えた」と見なすことができる



また、置換σは、2つだけを入れ替える置換に分解することができ、
その分解した置換の個数の偶奇で色々なことが評価できるらしい

ウルフラムαでネタバレしたから言えること・・・

有限深さの井戸型ポテンシャルでの、変位xの2乗の期待値を求めようとして、つまずきかけたのね


こういう式が出てくるじゃないですかー(しらんがな)

ここのb→∞の極限を取るところであっさりつまずきやがりまして、ウルフラム先生にご登場願ってしまったわけですけども

どうもロピタル使えばいけるっぽいことを知ってしまって
でもロピタルの発動条件ってどんなんだっけ？って、ねえ？



心から数学してねーな俺って感じでごめんねごめんねー＞人＜

固有値多項式(方程式)とトレース

10行10列




9行9列




任意の複素行列対応ですよ奥さん！

トレースの値が、n次多項式のn-1次の係数の符号反転になるんです！
(scilab)
いやぁすごいものを見させてもらいました。

理屈は余裕が出たら書きますが、僕が見たページでは、余因子展開してたっぽいです。
その余因子展開の方法が、最初の行の1列目から順番に一行全部やってるっぽいですねｗ
σとかSnとかsgnとか、説明もなしに初めて見る記号たちがいましたよ・・・
置換行列だったっけ

ほかの次数のは無理としても、せめてn-1次は出せておきたいね。

あとn-2次と0次な～

少なくともエルミート行列にしなきゃならんし
ノルム=1も定義しないといかんだろうし
たぶん特殊ユニタリ生成子限定の定理になるんじゃねーかな、n-2次の係数は。

固有値方程式 トレース 定理

n次の固有値多項式(方程式)のn-1次の係数がトレースになる
っての、タイトルみたいなググり方で見つけたっぽい！


ノルムがn-2次の係数に一致するのはどうググればいいものかなー


ノルムっつうぐらいだから、「特殊ユニタリ生成子の」ぐらいの制限があるのかも


思えばなんでこの沼に足突っ込んじゃったんだっけな
中途半端に経験あるもんだから、途中でちょっと天狗になって、車輪の再生産工場気取ってみたりして、ネタがなくて人にも聞けず、自爆


元々はパウリ行列、ゲルマン行列、クォータニオン界隈のイメージを自分がしたくて
なんかいい例題はないかって考えたのがパウリ行列版オイラーの公式だったんだよな

はぁ・・・初心に帰るべきだよなぁ

成立しなかった定理

であるときでさえ

この2つの固有値群3つ2セットは一致しない

ただし、2次行列の固有値λ1とλ2は以下のように定義されるとする


はい、ダメでした。なんてこったー
scilabに乱数ぶち込んで数値的に確かめてみて、こうもやすやすとダメだったとはね
2次と3次の時点でダメなら、4次以降の任意の次数に拡張とか無理じゃん



それにしても
忍者ブログの不具合治ってよかった～！

いやまじで僕だけに課せられたペナルティかと疑い始めてましたからね！
平和の象徴オールマイトの目が影に隠れてるのは乳首のメタファーとか言ったせいで社会から抹殺されたんじゃないかってビクビクしてたからね！


それか、先日間違えたコレを訂正させないための陰謀論かと思いましたからね！

高次方程式にありがちな、n-1次とn-2次の項の係数(特殊ユニタリ生成子)

n次の特殊ユニタリ生成子の固有値方程式がな？


ちゃんとこうなるのかまだ証明できてねえのよ。

どーにか証明に近づけないか、あがいたり昼寝したりしてたんだけど
今朝、ふといいことを思いついてな


こんな風に、ゲルマン行列の中にパウリ行列みたいのがマトリョーシカになってる構造を利用して
固有値使ってなんか楽できねーかなーって

まあ、中に入ってる2次の行列が、パウリ行列に似てるけど厳密にはちょっと違うから
話はそんな単純じゃないんだけどね

固有値x1とx2を求める方程式は以下のようになるよ


解(固有値)はこうなる

このx1,x2を

にぶち込むと、ちゃんと目当ての

こういう式にたどり着く分けです。


今回は3次のゲルマン行列を例にしたから大したことなかったんだけど、
4次行列の固有値を求めたくなったら、結構楽できそうなわけですよ。

ただなー、これをn次一般に拡張できるかが問題で
たとえば5次の生成子だったら4次方程式を解かにゃならん
ってなったらただの先延ばしだし
6次の生成子だったら5次方程式になってなんのこっちゃってなるわけよ。
というか、5次方程式を出すために5次の行列式で死ぬ思いしそうなんだけど


今考えてるのは、解と係数の関係を使うこと。
これで再帰構造の数珠つなぎになってぱっかーんとスカラーまでフィーバー起こればありがたいんだけど、なるかなあ？

目的はあくまで方程式の解じゃないんだよな、n-1次の係数がゼロになって、n次とn-2次の係数が逆相になるって証明がしたいのよねえ

scilabによるSU(9)の実装

昨日のブログに書いた、9次の特殊ユニタリSU(9)の実装を、scilabで行ってみます。
scilabは、信号処理に特化したオープンソースのプログラミング言語なので
複素行列の扱いが得意です。


まず80個の自由度を持つので、0から1までの間の一様（？）乱数を80個生成しましょう。
1行80列のベクトルsとして生成します。

s=rand(1,80)

次に、このベクトルのノルム(長さ)を1にします。
s=s/norm(s)

できあがったら一応、長さが1になっていることを確認しましょう。

norm(s)


9つの対角成分の入力をします。
x1～x9の9つです。

x1=s(3)+s(8)/sqrt(3)+s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x2=-s(3)+s(8)/sqrt(3)+s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x3=-2*s(8)/sqrt(3)+s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x4=-3*s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x5=-4*s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x6=-5*s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x7=-6*s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x8=-7*s(63)/sqrt(28)+s(80)/sqrt(36)

x9=-8*s(80)/sqrt(36)

ここで、計算ミスがないように、検算をしてみましょう。
x1からs9までの和はゼロになるはずです。



次に、非対角成分のうちの上三角部分を定義しましょう。

a12=s(1)+%i*s(2)

a13=s(4)+%i*s(5)
a23=s(6)+%i*s(7)

a14=s(9)+%i*s(10)
a24=s(11)+%i*s(12)
a34=s(13)+%i*s(14)

a15=s(16)+%i*s(17)
a25=s(18)+%i*s(19)
a35=s(20)+%i*s(21)
a45=s(22)+%i*s(23)

a16=s(25)+%i*s(26)
a26=s(27)+%i*s(28)
a36=s(29)+%i*s(30)
a46=s(31)+%i*s(32)
a56=s(33)+%i*s(34)

a17=s(36)+%i*s(37)
a27=s(38)+%i*s(39)
a37=s(40)+%i*s(41)
a47=s(42)+%i*s(43)
a57=s(44)+%i*s(45)
a67=s(46)+%i*s(47)

a18=s(49)+%i*s(50)
a28=s(51)+%i*s(52)
a38=s(53)+%i*s(54)
a48=s(55)+%i*s(56)
a58=s(57)+%i*s(58)
a68=s(59)+%i*s(60)
a78=s(61)+%i*s(62)

a19=s(64)+%i*s(65)
a29=s(66)+%i*s(67)
a39=s(68)+%i*s(69)
a49=s(70)+%i*s(71)
a59=s(72)+%i*s(73)
a69=s(74)+%i*s(75)
a79=s(76)+%i*s(77)
a89=s(78)+%i*s(79)

要素を書きだしたら、これを行列に入れます。

いっぺんにやると凡ミスをしやすいので、1行ずつ入れることにします。
B1=[0,a12,a13,a14,a15,a16,a17,a18,a19]
B2=[0,0,a23,a24,a25,a26,a27,a28,a29]
B3=[0,0,0,a34,a35,a36,a37,a38,a39]
B4=[0,0,0,0,a45,a46,a47,a48,a49]
B5=[0,0,0,0,0,a56,a57,a58,a59]
B6=[0,0,0,0,0,0,a67,a68,a69]
B7=[0,0,0,0,0,0,0,a78,a79]
B8=[0,0,0,0,0,0,0,0,a89]
B9=[0,0,0,0,0,0,0,0,0]

A=[B1;B2;B3;B4;B5;B6;B7;B8;B9]

こうして、上三角部分は埋まりました。

エルミート行列なので、下三角は上三角の複素共役です。

A=A+A'

ここまでで、非対角成分すべてが埋まりました。

念のため、ここでエルミートであることを確認してみましょう。
A-A'=0

対角要素を入れます。
A(1,1)=x1
A(2,2)=x2
A(3,3)=x3

A(4,4)=x4
A(5,5)=x5
A(6,6)=x6
A(7,7)=x7
A(8,8)=x8
A(9,9)=x9


これでエルミート行列Aが完成したので
今一度A-A'=0が成立するのを確かめておきましょう。


ここで、いきなりexpm(%i*A)とやってもいいのですが、せっかくなので寄り道して

Aの固有値方程式と固有値を見てみましょう。

変数xを定義します。
x=poly(0,"x")

次に、固有値方程式Bを見てみましょう。

B=det(A-x*eye(9,9))
detは行列式、eye(n,m)はn行m列の単位行列を意味します。

 0.0000013 + 0.0000969x + 0.0011569x^2 - 0.0025147x^3 - 0.0637862x^4 - 0.1498150x^5 + 0.2886920x^6 + x^7 - x^9

しっかりと、7次の符号が9次の符号の逆相になっていて、8次の係数がゼロになっていますね。
係数がちゃんと実数だけでできていることも確認しておきましょう。

この9次方程式の解はすべて実数です。エルミート行列であることからそういえます。

Aの固有値がこの9次方程式の解なのですが

spec(A)を計算してみるとわかります。

  - 0.6196769  

  - 0.4362829  

  - 0.2728734  

  - 0.2080910  

  - 0.0761649  

  - 0.0165342  

    0.1421505  

    0.4476562  

    1.0398167 

エルミート行列をただの数に例えると実数
歪エルミート行列は純虚数に対応するため
expm(%i*A)は、複素平面上の単位円に相当します。

これがユニタリ＝ユニットな行列と呼ばれる所以です。

C=expm(%i*A)
を計算してみて(行列指数関数などの場合は、関数名の後ろに行列(matrix)を意味するmを付ける必要があります)

spec(C)を計算してみますと

    0.5063784 + 0.8623114i  

    0.9014641 + 0.4328538i  

    0.8140662 - 0.5807721i  

    0.9899136 + 0.1416723i  

    0.9063287 - 0.4225735i  

    0.9630005 - 0.2694997i  

    0.9784271 - 0.2065925i  

    0.9998633 - 0.0165335i  

    0.9971009 - 0.0760913i  
このようになり、固有値がすべて複素平面における単位円周上にあることが確認できます。


また、このユニタリ行列は「特殊」ユニタリ行列といって、
abs(det(C))とやるまでもなく
det(C)そのものが1になるという性質があります。

つまり個の場合は9つすべての固有値の積が、ちょうど実正数の1になるということです。



なお、今回乱数にすべて正の数を用いたのは、
複素共役を見やすくするためでした。
Aの上三角の虚部がすべて正の数で、下三角の虚部がすべて負の数になります。

9次の特殊ユニタリSU(9)

なにをやっとるんだおれは



9行9列の行列なので、自由度は9×9-1=80個

なので、要素が80個あるベクトルsを用意し
2乗和（ノルム）を1にします。

対角要素をx1～x9とし、これらは実数
m行目n列目の上三角非対角成分をamnとして、実数ベクトルsから作った任意の複素数とします。

1から9までのxの和はゼロになるようにします。
つまり、作ろうとする行列Aのトレースはゼロです。

また、Aはエルミート行列なので、
下三角の非対角成分は、上三角の複素共役になります。


9次エルミート行列Aに虚数単位iをかけて、行列指数関数にすると
9次の特殊ユニタリ行列Uができます。

このようにして作ったユニタリ行列は「特殊」ユニタリと呼ばれ
行列式をとった時点で正の実数値1を取ります。

なお、一般的なユニタリ行列は行列式をとってさらにその絶対値を取らないと1になるとは限りません
つまり、一般的なユニタリUの行列式detUは、複素平面内では原点を中心とした半径1の円周上のどこかにいます。
ユニタリ行列の固有値も同様に、単位円上のどこかにいます。
固有値に関しては特殊ユニタリも一般的なユニタリも、正の実数値1とは限りません。



とまあ、ここまで、力技でSU(9)の計算の仕方を説明してきましたが
対応する物理的現象がおそらくまだ見つかっていないため(いつか見つかるの？)
とても虚しいです。
たとえていうなら、高次元の球である超球よりも虚しい
高濃度の無限と同じくらい虚しい感じでしょうか(よくは知らないんですが)


もっぱらの興味の対象は、4次ぐらいまでのユニタリのようですし
ほかにも特徴的な行列がたくさん定義されているようです。
僕は彼らの顔と名前がまだ全然一致していません



なお、対角要素xを形成するためのsの番号はよく
nを整数としてn^2-1番目と定義されるようです。

また、重みを付けるための分母のルートの中身は1から始まり
1,3,6,10,15,21・・・という数列で
差が1つずつ増える数列となります。
確か階差数列という名前でしたっけね


生成子Aの固有値方程式がですね
x^9-x^7+a6x^6+a5x^5+a4x^4+a3x^3+a2x^2+a1x+a0=0
(9つの実数解)
(あ、係数anはもちろん、全部実数です)

という形になって
SU(n)だとn-1次の係数がゼロ、n-2次の係数がn次の係数の逆相になるということをですね
なかなか証明できないでいるんです。めっちゃややこしそうで、なんか手が出ない
とりあえずテキトーなn=9ぐらいで確認したかったんですね

自動でSU(n)の生成子の固有値方程式をしらみつぶししてくれるプログラム、書けるかなあ

コサインとサイン・虚数と実数とごはんとおかず

今ググってもそうじゃないんだけど、前に確か、実フーリエ変換とかいうワードでググったら離散コサイン変換ってワードが関係づけられて、なんのこっちゃ？って思った覚えがあった気がする


じゃあ離散サイン変換ってないんか？って今ググったら、離散コサイン変化のwikiに亜種として載ってた


どういうわけか、三角関数のサインとコサインって
最初はサインをごはん、コサインをおかずとして扱うんだけど
いつの間にかコサインがごはん、サインがおかずになってるんだよな。

積分の時あたりまではサインがごはんなんだよ
積分した結果がアークサインになることはあっても
アークコサインとして表現することはめったにない

でも、複素数が出てきたあたりでごはんとおかずは入れ替わる。
コサインがごはんになるんだ。
たぶん、実数に相当するのがコサインのほうって定義からそうなるんだと思う
(純虚数がサインに相当)

興味深いことに、その実数と純虚数も、クォータニオンのあたりで
それまで実数のほうがごはんだったのに
いつの間にか純虚数のほうがごはんになるんだ

まあそりゃぁ3対1で虚数勢のほうが多くなっちゃうからね。数の暴力だよね。数学だけに。

2つのMADを最大幸福で同時に見る方法

過去ブログ参照
今日のやつ
今日のやつ２

やりたいのはだいたいこんな感じのことです
ざっくり9分割しかしてないので遠近斜め感が薄く、わかりづらいかもしれません

とにかく
上と下それぞれの水平ラインを水平に保ちつつ、遠近感を持たせることで
2つの画面をめいっぱい見せたいのです。水平ラインを上下中央ではなく、上端や下端にしたいのです
こういう機能実装してるソフトもうありますかねえ
もしかしたらExcelにすでに実装済みで、僕の経験が浅くてまだ見つけていないだけかもしれません

いきなりうまくいった

さっきの日記の続きなんだが

例題を作ったら即座にうまくいってしまった。
しかしこれはあくまで数値解だから一般性がほしい。
解析解を解くことにしよう。

いやそれとも、この枠に滑りこませられれば、一般性なんかどうでもいいのかもしれない。
要はフィルムをx,y：-1～1の間(z=0)の正方形にはめ込めばいいんだよ。
長方形バージョンを後から追加すれば近道できそうだな

フィルムに写る中身の映像はまあ、ブラウン運動あるいはジェットコースターみたいのでいいんじゃないかな

ああそうだ横の位置な。センタリングしないといかんよな

備忘録：消失水平線

前にもここで書いた気がするんだけど

こんな風に斜め奥に配置したスクリーンに投影する描写方法あるじゃんか。
フィルムというか写真というか、そういう紙媒体みたいなのを、
左は手前に、右は奥にって感じで描写する感じで
その絵が動くわけです。それがEDの演出になって、フィルムには本編の映像やらが写るわけ。
(シャンゼEDとかゼーガEDとか)


さっき久々にけもフレのMADを聞きながら散歩してて少女終末旅行を思い出してたんだけど
1つの音的なMAD作品に対して、2つの映像的MAD作品がある場合
なるべくデータ容量がかさばらないように保存したいとしたら
上の図みたいに、
左が手前で右が奥の画像を上に配置して
右が手前で左が奥の画像を下に配置すると、省スペースになるよね。

その際、上の画像では、フィルムの台形の上の辺が水平に
下の画像ではフィルムの台形の下の辺が水平になるように配置すると、最大幸福度になると思うわけですよ

じゃあその場合、適切な水平線を維持するためには、どこに消失点を配置するべきか
というのをずっと脳裏でモヤモヤ考えていまして
そもそもこの実験をする際、僕のいつも使ってるExcelの3Dビューア技術はどのくらい役に立つのだろうか
と思っていたわけです。


ここから備忘録なので、書いてる人にしかわからないことが多めになるかと思いますが
横=Ax/(z+z0)
縦=Ay/(z+z0)
の遠近法を取ったとして

この横を新たなx
縦を新たなy
z=1とか0とかで固定して回転させ、さらにその遠近法を取れば、Excel内だけで実験が閉じれるわけじゃないですか。

というのを思いつきました。

問題は、フィルム内の画像に用いるネタですよね。
わかりやすい例があればいいんですけどねえ
変換方法がわかったので、グラフ状になってしまえば動かすことは容易になりました。



このテストが完了したら、そうですね、マイクラや紙Excelみたく
縦と横のセル幅や高さをすんげーちっさくして
画像を取りこむマクロでも組むと面白いかもしれない。

VB（VBA？どっちだっけ？)も習いながら、「操作を記憶させる」マクロにもいい加減手を出したいな


たとえばそうだな、手始めに連続キャプチャした画像を、複数のシートに同じように分配したのち、同じようなトリミングをほどこして、一定時間置きにシートを移動するマクロとかできたらいいよね
これまで全部手動マクロだったからなー、腱鞘炎とか腰痛とかになりそうなんだよ


マイクラだか紙Excelだかはそのあとだなー
どうやって画像情報取りこむのかさっぱりだし
取りこんだら画像処理したいよなーそりゃあなー



やっぱりやることいっぱいだよなー
とりあえず散歩中にやることメモれてよかった。
ブログに書いたからって、スマホのメモはクローズしないでおこう。
記事が多すぎるから、記憶をブログだけに頼るとネタが迷子になる。


まあ焦らずゆっくり、技術や僕自身が老衰で死んだりしない程度にのんびり頑張ろう
もしかしたらAviUtlにそういうのもう全部揃ってるかもわからんけどね
車輪の再発明は僕の夢っつーか憧れっつーか、もうすでにライフワークだから。

めまいカットで理解してみる、遠近法の近似としての等角投影法

このブログではいつもやってることなのですが

たとえばこのように、x,y,z表示で、3軸ともに0から1のキューブのデータを作ったとします

これのx,yの平面グラフを作ると、そのまま等角投影法の立体としての表現もされるのですが

これを遠近法にするために、zにz0という下駄(カメラの引き具合ともいう)を履かせてから
倍率Aをxとyにそれぞれ掛け算して、各zにz0を足したもので割り算

つまり

このように計算することで、以下のように簡易的な遠近法表現を得ることができます。

ここで、先日の結果を用いて、めまいカットを適用してみましょう。

「カメラを近づけたり遠ざけたりしながら、ズームインしたりズームアウトしつつ
一定の距離の対象物の大きさだけを固定する」
というのがめまいカットの定義だと思うのですが

z00を固定したい対象物への奥行き距離
Cを任意の定数としますと

手前(奥行きz00=0)を固定した際はこのように

奥(奥行きz00=1)を固定した際はこのようになりますよね。



ここで、
キューブのある点(x,y,z)、ここでは(1,1,1)と(1,1,0)座標を指すのですが
zが同じな2セットのxとyがあったとして


遠近法めまいカットのグラフに表示する横と縦a、bは

こうなります。縦軸は横軸と同じなので略すとして
zだけが異なる2セットのxとyについて
横軸2-横軸1(z2とz1)がゼロになるz0を考えてみましょう。


計算して整理するとこうなりますが、
ここで分母分子をz0の2乗で割り算してみますと以下のようになって

この状態でz0→∞の極限を取ると
横軸2-横軸1=0となることがわかるかと思います。
縦軸についても同様です。

これは一体なにかといいますと、等角投影法そのものです。

つまり、無限遠点からめまいカットをやると、遠近法は等角投影法に近似されるのです

めまいカットの公式

僕はよく、Excelで簡易3Dビューアを作る際に

2Dグラフに
横=A*x/(z+z0)、
縦=A*y/(z+z0)

というアルゴリズムで遠近法にしてるんですけど

Aとz0というのが、物理的な意味としては
Aがカメラのズーム倍率で、z0がカメラがバックした位置なわけです。


そこで、バックしながらズームすると、ヒッチコックが考案したとされる「めまいカット」になるわけですが

奥行きzの対象物の大きさを固定したまま、Aとz0を連動して動かしたいですよね。


そこで、横あるいは縦を、z0で微分して、ゼロになる条件を探すわけです。
(縦の式は同じなので略します)



展開して整理したら、こんな微分方程式ができますよね。

じゃあこれを解きます。
任意の定数をCとして

このような式が導かれました。


では実装してみましょう。
x,y,zが-2から2までの範囲に、1ずつの目盛で格子を作り
その中に螺旋を入れてみました。


横の右向きがプラスx、縦の上向きがプラスy、手前から奥がプラスzです。


z=0を固定対象とすると、このように、手前の物体と奥の物体が、z=0の平面に双方近づいたり遠ざかったりします。

z=-2だと、奥の物体が一方的に近づいたり遠ざかったりし


逆にz=2だと、手前の物体しかないので、手前の物体だけが近づいたり遠ざかったりします。


z=-1だと、手前より奥の物体の方が少し多い感じで、近づいたり遠ざかったりします。