読者です 読者をやめる 読者になる 読者になる

読み書きプログラミング

日常のプログラミングで気づいたことを綴っています

(61) Newtonの重力理論と楕円軌道

Newtonの万有引力理論では、二体問題の場合、軌道は楕円になります。 これを順を追って考えます。 角運動量とエネルギーが保存することから以下の式が時間に依存しなくなります。 について解きましょう。 /* 9.4.10m */ derivs : solve([(m/2)*(rDot^2 + r^2…

(62) 水星の近日点移動

一般相対論的な重力場での運動の場合、軌道は楕円にならず、近日点は移動します。これを導くため、Schwarzschild 解から出発します。Schwarzschild 解は以下のような重力場の方程式の球対称解です。この計量で測地線が満たさなければいけない方程式は以下に…

(60) ピンポン球の重力

半径aの薄い球体殻が生じるNewton重力は、殻の外の重力ポテンシャルに関して中心にある質点の重力と等しいことを確認しましょう。 /* 9.4.5m */ assume(not equal((r - a)*(r + a), 0))$ integral : integrate(a * (2 * %pi * a * sin(u)) / (r ^ 2 + a ^ 2 …

(59) Riemannゼータ函数の周辺(パス)

"Crandall"では、Riemannゼータ函数の零点の数値的観察の後、 有限級数と誤差評価による非零点領域の確認 量子振動子の波動函数と余弦積分変換による零点導出 Riemann予想が成り立つ場合のπ(x)の近似からの素数の数の見積もり と続きますが、スクリプトは移…

(58) Riemannゼータ函数の零点

Riemannゼータ函数の絶対値の2乗の逆数をプロットしてみましょう。 /* 9.3.8m */ plot3d(1/cabs(zeta(x +%i*y))^2, [x, 0, 1], [y, 10, 24], [grid, 50, 50], [z, 0, 40])$ x=0.5上のy=14, y=21付近に発散していく場所があることが見えます。実部を0.5に固定…

(57) ゼータ函数のEuler積

以下の恒等式が成り立ちます。右辺がEuler積と呼ばれる形です。これをMaximaを使って実感してみましょう。表示をシンプルにするため、sに-sを代入し、Euler積の中身をMaclaurin展開して、検証する次数まで取り出します。 /* 9.3.3m */ proofLimit : 15$ f(x,…

(56) Vandiverの判定基準

Fermat予想に対して、Kummerに続いて、Vandiverが、非正則素数pに関して以下を証明したそうです。 のうち、分子が素数pで割り切れるものをとします。 pより大きくp(p-1)より小さい素数で、rp+1(rは自然数)の形に書けるものPがあり、かつ、となる自然数uが存…

(55) 素数の正則性チェック(要加筆)

Fermatの最終定理が証明されるまでの歴史のなかで、Kummerが正則素数に限定したものを証明したそうです。正則素数pは、Bernoulli数の偶数項の分子のいずれも割り切らない素数という性質と同値であることがKummerによって示されたそうです。 では、素数の正則…

(50) スペクトル密度

波形とスペクトル密度の関係を見ます。通常は、パワースペクトルを見ますが、振幅の小さなところの特徴も見るため、振幅スペクトルを見ます。まず、振幅スペクトルを計算する関数とそれをプロットする関数を用意しましょう。 /* 8.2.9m */ load("fft")$ matc…

(52) 無限インパルス応答

以下のような無限インパルス応答(IIR)ディジタルフィルタについて考えます。c=0.95でgを0から0.9まで変化させた時のゲインをプロットしてみましょう。 /* 8.3.6m */ c : 0.95$ a : [2 * c * g, -c ^ 2]$ p : length(a)$ b : [-(2 * c * g) / (1 + c ^ 2)]$ q…

(54) 画像処理(パス)

"Crandall"には画像データのエッジ検出がありますが、Maximaで画像データを扱うライブラリを存じないので、パスします。 参考文献 Crandall, Mathematica―理工系ツールとしての (アジソン ウェスレイ・トッパン情報科学シリーズ) p.261-p.265

(53) 有限インパルス応答

有限インパルス応答(FIR)フィルタの場合、バンドパスフィルタは、周波数と共鳴幅を決めるとほぼ決ります。伝達函数は以下のような感じです。特定の周波数に関して、この伝達函数での周波数応答をプロットしてみましょう。 /* 8.3.13m */ sinc(x) := ( local(…

(51) 音声処理(パス)

"Crandall"には音声データのソノグラム可視化がありますが、Maximaで音声データを扱うライブラリを存じないので、パスします。 参考文献 Crandall, Mathematica―理工系ツールとしての (アジソン ウェスレイ・トッパン情報科学シリーズ) p.247-p.252

(49) 論理回路のKarnaugh図

以下の論理回路を例に、考えます。 この回路のKarnaugh図を生成してみましょう。 /* 8.1.14m */ nand(x, y) := not (x and y)$ logic(x, y) := nand(nand(x, y), nand(not x, not y))$ apply(matrix, makelist(makelist([a, b, logic(a>0, b>0)], a, 0, 1), …

(48) 定電圧源回路の静特性

以下の最も簡単な定電圧源回路を考えます。これは非線形回路の例になります。 対応する方程式は以下です。二番目の式は、Shockleyのダイオード方程式です。は飽和漏れ電流、qは素電荷、Tは絶対温度、kはBoltzmann定数です。室温(T=298K)、抵抗100Ωの条件で、…

(47) 能動フィルタを題材とした線形回路分析

能動フィルタのインピーダンスを求めます。 具体例として以下の回路を考えます。回路図エディタとして、回路シミュレータのPSIMを利用しました。PSIMは素子数と機能が限定されたデモ版がフリーで使うことができます。PSIMから上記回路のネットリストをファイ…

(46) 線形回路の定常解析

以下のようなLCR回路の周波数応答を考えます。 V1からV2を求める問題です。方程式を書き下すと、周波数領域でMaximaに解かせてみましょう。 /* 8.1.3m */ output : solve([%e ^ (%i * w * t) - v2 * %e ^ (%i * w * t) = i0 * %e ^ (%i * w * t) * r, i0 = i…

(45) Hodgkin-Huxley系の伝搬速度

ニューロンの活動電位の現象論的方程式Hodgkin-Huxley系を用いた軸索上の電位の伝搬速度について考えます。 出発点は、Hodgkin-Huxley系とケーブル理論になります。 Hodgikin-Huxley系はイオンチャンネルを加味した細胞膜の活動電位の方程式系です。: 膜電位…

(44) 種の絶滅

遺伝子ドリフトと類似の問題である種の絶滅について考えます。 2種類の種A, Bがあり、もし全個体数が一定nで、次の世代の占める割合の確率がそれぞれの個体数比に比例するとします。 Aの個体数がiの時、次の世代の個体数それぞれの確率は、Tは遷移確率行列と…

(43) 遺伝子ドリフト

ある遺伝子座に2種類の遺伝子A, aがある場合を考えます。この形質Aaの表現型はAA, Aa, aaになります。 個体数nの中の遺伝子Aの頻度をfA, 遺伝子aの頻度をfa(=1-fA)とし、この集団からランダムに配偶子を選んで次世代が個体数n生まれるとします。遺伝子Aの頻…

(42) 遺伝的組み換え

二組の遺伝子座で組み換えが起こる確率を組み換え価と呼びます。 以下の組み換え価のテーブルを例に、組み換え価から遺伝的距離を求めることを考えます。 遺伝子座の組 組み換え価(%) a-b 50 a-c 15 a-d 38 a-e 8 b-c 50 b-d 13 b-e 50 c-d 50 c-e 7 d-e 45 …

(41) ヘリウム原子の基底状態

ヘリウム原子の電子で、球対称な解について考えます。 ハミルトニアンはこの場合、となります。単位系は式がシンプルになるように取りました。 基底状態のエネルギーを求めるのに、以下の試行波動関数から出発します。は規格化されていると考えます。基底状…

(41) 水素原子の波動関数の有限要素法解法

水素原子のSchrödinger方程式は厳密解が知られていますが、より複雑な分子を検討する際の足場として数値解法を見ておくことは有意義です。 ここでは、角運動量が0の時の水素原子のSchrödinger方程式を有限要素法で解きます。対象となる微分方程式は、以下の…

(40) 化学平衡

従兄のジョーは化学工場の研究員である。彼は価値のある衛生商品を作り出す上で不可欠の材料であるSuperzapという物資の分析を任されている。SuperzapはZapとRegant Bから以下の反応によって合成される。 ZapとRegant Bよりずっと高価である。この反応定数K…

(39) 化学反応式

メタンが燃える反応では以下の変化が起こります。反応式の係数を求めるスクリプトを作りましょう。 /* 7.1.4m, 7.1.8m */ reactants : [[C, 1, H, 4], [O, 2]]$ products : [[C, 1, O, 2], [H, 2, O, 1]]$ atoms : []$ lexp : 0$ numvars : 1$ for r : 1 thr…

(38) Sierpinski三角形

自己相似な三角形であるSierpinski三角形は、色々な方法で生成することができます。 ここでは、カオスゲーム(ランダム反復函数系)を使ってみます。 ある点に対して、ランダムに(0, 0), (1, 0), (0, 1)のいずれかを選び、ある点と選ばれた点の中点を与える写…

(37) Mandelbrot集合

Mandelbrot集合は、以下の二次写像の反復函数系について、有界な解を与える複素数パラメータcの集合と定義されます。MaximaにはMandelbrot集合を扱うライブラリがいくつか用意されています。それらを使ってみましょう。 load("dynamics")$ load("draw")$ man…

(36) 分岐とFeigenbaum定数

離散力学系について考えます。(本質的にロジスティック写像と同じ)以下の写像を対象にします。例えば、c=1.23の時、初期値x=0から始めるこの繰り返し写像は2点からなる周期解に収束します。 これを見ておきましょう。 /* 6.3.2m */ load("dynamics")$ evolut…

(35) N-ソリトン解

KdV方程式の続きです。 KdV方程式の初期条件をポテンシャルとする定常Schrödinger方程式を解くと、KdV方程式の厳密解が求められることが知られています。 初期条件に以下を仮定します。これをポテンシャルとする定常Schrödinger方程式の解は以下になります。…

(34) ソリトン、数値解、厳密解

Korteweg-de Vries方程式を調べます。これは、非線形微分方程式ですが、可積分であることが知られています。 1ソリトンの厳密解として、以下が知られています。確認してみましょう。 /* 6.2.3m */ u : -2*sech(x - 4*t)^2$ trigsimp(diff(u, t) - 6*u*diff(u…

(33) 共振

強制力と減衰力を受けた調和振動子を考えます。運動方程式は以下の通りです。定常解の周波数は強制力の周波数ω/2πになるはずです。この時の振幅を求めてみましょう。 /* 6.1.21m */ x(t) := A*exp(%i*w*t)$ sol : solve(m*diff(x(t), t, 2) + 2*m*g*diff(x(t…

(32) 膜振動

円形の薄膜の振動を考えます。 極座標を使うと、運動方程式は以下になります。ここでvは音速です。この方程式は、時間と偏角に関して線形なので、変数分離して周波数領域で解くことができます。すなわち、として一般性を失いません。 運動方程式の解となるは…

(31) 固有振動

以下のように3つのバネで構成される2自由度系を考えます。 質量mのそれぞれの質点の釣り合いの位置からのずれを、バネ係数をそれぞれk,jとすると、運動方程式は以下となります。線形問題なので周波数領域で考えることが有効です。と置けます。すると運動方程…

(28) 均一場内の相対論的運動の厳密解

前回の問題の厳密解を求めます。 LagrangianからEuler-Lagrange方程式を導出し、積分ましょう。 /* 5.3.14m, 5.3.16m */ L : -m*c^2*sqrt(1 - (v/c)^2)-q*E*z$ sol : solve((diff(L, v) = (integrate(diff(L, z), t) + c0))^2, v)$ z : integrate(radcan(v),…

(30) Lorentz計量とCompton散乱

二つの慣性座標系がx軸方向に速度vですれ違っているとします。 ラピディティgを以下のように定義します。すると、二つの慣性座標系間のLorentz変換はと表されます。 これは、Lorentz計量を不変に保ちます。 これをMaximaで確認してみましょう。 /* 5.3.31m *…

(29) 一般相対性理論を垣間見る

「Mathematica―理工系ツールとしての (アジソン ウェスレイ・トッパン情報科学シリーズ)」での均一場内の相対論的運動の元々の問いかけは、 均一な重力場の中で、固有時間が最大となる質点の運動は何か。 でした。 一般相対性理論では「均一な重力場=一定加…

(27) 均一場内の相対論的運動の摂動解

z方向に一様な電場の中で時刻と時刻で同じ場所に到達する相対論的電荷qの運動を考えます。簡単のため空間は1次元とします。 二つの時刻の位置が運動の条件になっている場合、Lagrangianと最小作用の原理で考えるのが適しています。 作用は以下のように定義さ…

(26) 相対論的エネルギー

質量mの物体が速度vで移動している時の相対論的総エネルギーは以下のように書けます。総エネルギーを速度で展開してみましょう。 /* 5.3.2m */ E(v) := m*c^2/sqrt(1 - (v/c)^2)$ taylor(E(v), v, 0, 8); 初項には、光速度の二乗を係数に持つ最も美しい物理…

(25) 軌道狙い撃ち法

調和振動子の定常状態のSchrödinger方程式は、なるg(x)を使うことで、形式的に以下の1階微分方程式に変換できます。ψ(x)からg(x)からへの変換式を見ると、ψ(x)が偶函数の時g(0)=0、ψ(x)が奇函数の時であることがわかります。また、ポテンシャルV(x)が調和振…

(24) 摂動計算

位置エネルギーを以下のように量子調和振動子からの1次の摂動で考えます。エネルギー順位の摂動展開は、以下のようになります。ここで、は Hermite多項式について以下の積分公式が知られています。ここで、。を計算するために、をHermite多項式で展開しまし…

(23) 量子調和振動子の離散固有値

1次元の量子調和振動子を考えます。単位系をうまく選んだ定常状態のSchrödinger方程式は以下の通りです。中心から遠いところの振る舞いを考慮すれば、と置き換えるのが妥当です。その結果、Schrödinger方程式は以下になります。級数展開を使って、g(x)の偶函…

(22) トンネル効果

1次元のSchrödinger方程式は簡単な単位系を選択すると以下のように書けます。ここで、d)" />という0" align="middle" />の位置エネルギーを考えると、が粒子の運動エネルギーより大きな場合がトンネル効果が現れる場合になります。 この条件では、波動関数は…

(21) 振り子の周期

振り子の運動について考えます。振り子のHamiltonianは以下のように書けます。振り子の正準運動方程式はHamiltonianから以下の通りに求められます。ここで、lは角運動量、θは角度、mは質量、gは重力加速度、Lは振り子の長さです。 Euler法でMaximaに周期を計…

(20) 運動方程式の級数展開による近似解法

調和振動子の運動方程式は以下の通りです。ここで、pは運動量、qは位置、mは質量、kはバネ係数です。 これは解析的に解くことができます。 初期条件をp(0)=0, q(0)=1として、Maximaに解かせてみましょう。 assume(m > 0, k > 0)$ atvalue(q(t), t = 0, 1)$ a…

(19) 中国の剰余定理アルゴリズムの高速版(要加筆)

互いに素な数で割った余りがそれぞれ、であるようなnはを法として一意に存在します。(中国の剰余定理)解nは以下の式で求めることができます。ここで、、、はを法とするの逆元です。すなわち、を満たします。 これをMaximaで素直に書き下してみましょう。 p :…

(18) Newton法による多項式の逆数のMaclaurin展開

Newton法を使って逆数を求めるには、以下の漸化式を使います。これは多項式にも適用することができます。具体的な例として、Bernouli数を計算してみます。Bernouli数は以下の式で定義される数列です。微分してt=0を代入すると不定になるので代数的に求めるこ…

(17) 大きな整数の高速乗算

例えば、8ビット*8ビットの乗算器を持つマイコンを想像してください。このマイコンで2つのmバイトの整数の掛け算を行うには、普通、回、乗算器を使います。 この掛け算を、のオーダーの乗算回数で行うアルゴリズムを考えます。 一般化して整数x, yを基数wで…

(16) 素数テスト

Maximaには素数テストの組み込み関数primepがあります。これはMiller-Rabinの確率的アルゴリズムを用いたものです。 決定的アルゴリズムであるLucasの素数テストを実装して、の素数性を証明しましょう。 /* 4.3.15m */ p : (10^23 - 1)/9$ factors : ifactor…

(15) Lenstraの楕円曲線法による素因数分解アルゴリズム

6番目のFermat数の因数を楕円曲線法で探してみましょう。 /* 4.3.13m */ /* 因数発見を示すグローバルスイッチ */ notFoundYet : true$ /* Maximaの組み込み関数gcdexは、最大公約数の符号が不定なので、正となるように規格化する関数を用意。 */ extendedGC…

(14) Pollardの素因数分解アルゴリズム

Maximaは素因数分解を実行する関数が用意されています。使ってみましょう。 /* 4.3.1m */ factor(54444439); 使用した機能 factor マニュアルによると、factorはデフォルトでは、Berlekampのアルゴリズムというものを使うそうです。 組み込み関数を使う代わ…