読み書きプログラミング

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

Crandall

(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のアルゴリズムというものを使うそうです。 組み込み関数を使う代わ…

(13) 留数定理を使った複素積分

複数の極を持つがそれ以外では正則な函数の単純閉曲線上の積分は留数定理を使って計算することができます。 留数定理を活用した複素積分の計算をMaximaでやってみましょう。 例は、0" align="middle" />を以下の経路で積分した場合です。 外周の積分は極限で…

(12) Scherk曲面の曲率

で表される曲面の各点の平均曲率Cは以下のように定義することができます。 球面の平均曲率を計算してみましょう。 /* 4.2.13m */ assume(r > 0)$ f(x, y) := sqrt(r^2 - x^2 - y^2)$ fx : diff(f(x, y), x); fy : diff(f(x, y), y)$ den : sqrt(1 + fx^2 + f…

(11) の導出

を導出します。 以下のように周期2の函数を定義します。は実定数です。 のフーリエ変換を求めてみましょう。 /* 4.2.6m, 4.2.7m, 4.2.8m, 4.2.9m */ load(fourie) f(x) := a*x^2 + b*x^4$ f(x) = totalfourier(f(x), x, 1); 使用した機能 totalfourier を代…

(9) テータ定数の恒等式

Jacobiの三番目のテータ定数は以下のように定義されます。(はqの函数ですが、テータ函数の原点の値であることから定数と呼ばれます。) この函数に関して、以下の恒等式が成り立つことが知られています。これを有限項に限定して確認してみましょう。まず、右…

(10) 円周等分多項式

円周等分多項式は、の複素数解のうち、n乗すると初めて1になる解を根に持つ既約多項式です。正確な定義は以下の式になります。円周等分多項式をMaximaでいくつか書き出してみましょう。 /* 4.1.24m */ cyclotomic(k, z) := trigrat(product(if gcd(j, k) = 1…

(8) Poisson和公式を用いた無限級数の計算

整数全体に渡る以下の和を計算します。函数に以下のPoisson和公式を適用するところから始めます。 最後の積分をMaximaに解かせてみましょう。 /* 4.1.10m */ term : integrate(exp(-y)*cos(m*y), y, 0, inf); 使用した機能 integrate inf 従って、以下が成り…

(7) Legendre多項式

Legendre多項式は、以下の母函数を展開することで得られます。この性質から、Legendre多項式は双極子ポテンシャルを距離の逆数で展開する際に現れます。早速展開してみましょう。 /* 4.1.2m */ legendreGen : taylor((1 - 2*z*t + t^2)^(-1/2), t, 0, 4); 使…

(6) 3次元パラメトリックプロット

三つ葉結び目を描いてみましょう。 まずは2次元平面への投影図です。 /* 3.6.4m */ a : 1.0$ b : 0.4$ r(t) := a + b*cos(3*t/2)$ plot2d([parametric, r(t)*cos(t), r(t)*sin(t), [t, 0, 4*%pi + 0.1], [nticks, 200]], [x, -2.0, 2.0], [y, -1.5, 1.5])$ …

(5) パラメトリックプロット

パラメトリックプロットの例として、等角螺旋を描いてみましょう。 /* 3.5.1m */ r(t) := exp(-0.1*t)$ plot2d([parametric, r(t)*cos(t), r(t)*sin(t), [t, 0, 30], [nticks, 500]], [x, -1.5, 1.5], [y, -1, 1]); 螺旋状の点と原点を結んだ線とその点の接…

(4) 3次元剛体アニメーション

正五角錐を回転させるアニメーションを製作しましょう。 /* 3.4.1m.mac */ load(draw)$ xyz2polar(v) := block([r : sqrt(v[1]^2 + v[2]^2 + v[3]^2)], [r, acos(v[3]/r), mod(atan2(v[2], v[1]), 2*%pi)])$ t : makelist(ev([cos(a), sin(a), 0], a = 2*%pi…