誤差関数によるステップ関数でGibbs現象は起こるか?
ステップ関数等で不連続に打ち切られた関数をフーリエ変換しようとすると、どんなに頑張っても振動が残る。
これはギブス現象として知られている。
ギブズ現象 - Wikipedia
では、ステップ関数の代わりに誤差関数で滑らかにしたら、どれくらい収束が良いのか?
ちょっと気になったので、やってみた。
まず、誤差関数で作る擬ステップ関数等を定義する。
import numpy as np from scipy.special import erf from scipy.integrate import simps from matplotlib import pyplot as plt # pseudo step function def pstep( x, x0, hdx): return 0.5 * ( 1 + erf( ( x - x0 ) / ( 0.5 * hdx ) ) ) # pseudo rectangular def prect( x, L ): hdx = L / 16 return pstep( x, L / 4, hdx ) - pstep( x, 3/4 * L, hdx ) # rectangular def rect( x, L ): func = np.zeros( len( x ) ) for i, xx in enumerate( x ): if L / 4 <= xx <= 3/4 * L: func[ i ] = 1.0 return func
滑らかな矩形波をフーリエ変換する。
フーリエ変換はsimpson積分で自前で頑張る。
L = np.pi N = 120 x = np.linspace( 0, L, N ) k = 2 * np.pi / L * np.arange( -100, 100 + 1 ) func = prect( x, L ) # func = rect( x, L ) fn = [] for kk in k: integrand = func * np.exp( -1.j * kk * x ) fn.append( complex( simps( np.real( integrand ), x ), simps( np.imag( integrand ), x ) ) / L ) plt.plot( x, func ) fourier = 0. for i, f in enumerate( fn ): fourier += f * np.exp( 1.j * k[ i ] * x ) plt.plot( x, fourier.real ) # fourier.imag is almost zero. plt.show()
結果元の関数をほぼ完全に再現している。
程度なので、そんなに大変じゃない。
通常の矩形波に対して同じ条件でやってみた。収束する気配は無い。
ちなみにを10倍に増やしても、振動が細かくなるだけで、消えることはない。
滑らかにする大切さを知りました。
Juliaで一次元井戸型ポテンシャル
以下のサイトの下の方に、Juliaで一次元のシュレーディンガー方程式を解くPDFが紹介されている。
物理ノートby永井
Juliaの練習としてやってみた。
PDF内では、無限の井戸の中に斥力ポテンシャルを入れた場合をやっているが、ここでは引力ポテンシャルに対してやってみる。
解き方としては、中心差分を用いて二階微分を近似した有限要素法に対応する(と思っている)。
とりあえず、引力ポテンシャルは滅茶苦茶深くした。
using Plots using LinearAlgebra function calc_V( N, a, V0 ) vec_V = zeros( Float64, N ) center = N * a / 2 L = N * a / 2 for i in 1 : N x = i * a if center - L / 2 <= x <= center + L / 2 vec_V[ i ] = - V0 end end return vec_V end function make_H1dv( N, a, V0 ) mat_H = zeros( Float64, N, N ) vec_V = calc_V( N, a, V0 ) for i in 1 : N for dx in -1 : 1 j = i + dx v = -1 / ( 2 * a^2 ) if dx == 0 v = 1 / a^2 + vec_V[ i ] end if 1 <= j <= N mat_H[ i, j ] = v end end end return mat_H end function main( N, a, V0 ) N = 100 a = 1 / N V0 = 1e4 mat_H = make_H1dv( N, a, V0 ) energy, phi = eigen( mat_H ) return energy, phi end energy, phi = main() plot( energy ) savefig( "energy.png" ) plot( phi[ :, 1:5 ] ) savefig( "phi.png" )
最初の固有値25個くらいは値が負なので、束縛状態であることがわかる。
最初の5個くらいの波動関数を確認すると、しっかり中心の有限井戸に束縛されている。
プロットしてみるとわかるが、energyの図の真ん中らへんは、無限井戸も含めた全体に広がる解で、さらに高いエネルギーの解は有限井戸を避けて片側に強く局在した状態になる。
片側だけに寄るのは対称性的におかしいが、そこはまぁ数値解の御愛嬌という感じだろうか。
とりあえず、デカイ箱を用意しておいて、その中に引力ポテンシャルを入れれば、束縛解をパッと出せることがわかった。
単振動をT行列を使って解く。
前回、Green関数を使って古典単振動の軌跡を求めた。
koideforest.hatenadiary.com
今回は、無限級数の別表現として、行列を使ってみる。
行列は、Green関数を用いて次のように定義出来る。
この行列を用いて、前回の式を書き直すと、 だから、
ここで、Green関数のインパルス応答としての役割を考えると、作用、外力に対して、
と現わせることから、外力を行列で表すと、
外力を(意味わからんが)「(外力が無い時の)位置に対する応答」としてみなせそうな雰囲気を感じる。
そろそろ実際に行列を求めてみる。
Green関数は前回と同様に、以下のものを使用する。
したがって、
は露わには含まれていないが、 にこの条件が暗に含まれているため、明記化しておく。
以上をまとめると、
前回の境界条件を引き続き採用すると、であるから、
よって、強制振動が外力として加わっているとみなすことが出来る。
無限級数の行列による置き換えは、以下のような関係を求めたのと同様であるから、この時点で方程式は解けていると言っても良い。
よって、解が得られた。
もちろん、真面目にGreen関数との積を積分して求めても良い。
計算量は多いが、統一的に扱えるのは理解の助けになると思われる。
単振動を(非摂動)Green関数を使って解く。
前回、Green関数を使って古典力学の基礎問題を解いた。
koideforest.hatenadiary.com
koideforest.hatenadiary.com
今回は古典単振動の軌跡をGreen関数を使って求める。
自由落下の時は、非斉次項が定数だったが、単振動では求めたい関数自身が含まれているため、固有値問題のような形になっている。
そのため、量子力学での摂動展開のような形式で解くため、Green関数のイメージを掴むのにとても良いと思う。
単振動の方程式は以下の通り。
答えが分かっているので、安心して問題を解き進められる。
Green関数は、境界条件から、自由落下の時と同様のものを使えば良いことがわかる。
よって、形式的な解は、
注目するべきは、解くべき が積分の中で再登場している点である。
量子力学で言えば、リップマン・シュウィンガー方程式と同様の形をしている。
これを解くために、 に右辺をまるごとそのまま再代入する。
代入を繰り返すと、無限級数が得られる。
ここで、右辺第二項を計算してみると、
よって、少なくともの二次以上になる。
境界条件はで与えられるので、積分を含む項はで消えてしまい、の決定に寄与しないことがわかる。
これにより、に含まれる係数が求まり、は時間に依存しないことがわかる。
これらの結果の整理すると、
したがって、問題はGreen関数の積の積分ということになる。
いくつか計算してみる。
したがって、以下のような一般化が予想できる。
よって、 と求めることが出来た。
自由落下をGreen関数で解く。
二階微分だけの演算子に対するGreen関数を求めた。
koideforest.hatenadiary.com
求めたと言っても、斉次解が含まれていないので、Green関数の一般解にはなっていない。
斉次解を含めると、
ここで、自由落下の問題を考える。
この時、Green関数を用いると、
一応、確認すると、
Green関数の境界条件を以下のように設定すると、
本当にこれで良いのかは知らない。。。
これを用いると、
だから、
よって、ちゃんと自由落下の軌跡である が得られた。
二階微分だけの演算子のGreen関数(フーリエ変換経由)
前回、単振動方程式におけるGreen関数を導出した。
koideforest.hatenadiary.com
二階微分だけとなると と置き換えることに対応するのは明らかである。
実は、ここから直接求めようとすると、二位の極をまともに扱わないといけないので、死にかけた(というか死んだ)。
そのため、単振動のGreen関数から の極限を取ってあげると、
と素直に求まる。
一応、確認すると、
「一度、困難を外す」というのは時に重要だと感じた次第である。