Loading [MathJax]/jax/output/HTML-CSS/jax.js

2019년 9월 14일 토요일

파이썬 SymPy 기호수학 - 7. 퓨리어 급수와 변환

파이썬 Sympy 기호수학 - 7. 퓨리어 급수와 변환

SymPy 사용법

SymPy 라이브러리를 불러오고, 기호 변수를 선언한다. 출력을 LaTex 수식으로 나타내고, 맷플롯립 모듈을 불러온다.

In [1]:
from sympy import *
          
x, y, z, t = symbols('x y z t')
f, g, h = symbols('f, g, h', cls=Function)

init_printing()
%matplotlib inline 

퓨리어 급수

주기가 2π인 함수 f(x)를 퓨리어 급수로 나타내면

f(x)=a0+n=1(ancosnx+bnsinnx)

a0=12πππf(x)dxan=1πππf(x)cosnxdxbn=1πππf(x)sinnxdx

[예제] Heaviside 계단 함수를 퓨리어 급수로 표현하기

In [2]:
Heaviside(x)
Out[2]:
θ(x)
In [3]:
plot( Heaviside(x) )
Out[3]:
<sympy.plotting.plot.Plot at 0x222f62525c0>

Heaviside 계단 함수를 (π,π) 구간의 주기 함수로서 퓨리어 급수를 전개한다.

In [4]:
fourier_series( Heaviside(x), (x,-pi,pi) )
Out[4]:
2sin(x)π+2sin(3x)3π+12+

처음 4개의 항을 취하여 그래프로 나타내보면

In [5]:
s = _
s.truncate(4)
Out[5]:
2sin(x)π+2sin(3x)3π+2sin(5x)5π+12
In [6]:
plot( s.truncate(4), (x,-2*pi,2*pi) )
Out[6]:
<sympy.plotting.plot.Plot at 0x222f6743b70>

10개의 항을 취하면, 계단함수에 더욱 가까워진다.

In [7]:
plot( s.truncate(10), (x,-2*pi,2*pi) )
Out[7]:
<sympy.plotting.plot.Plot at 0x222f65a3da0>

[예제] 함수 f(x)=x(π,π) 구간의 주기 함수로서 퓨리어 급수를 전개한다.

In [8]:
s = fourier_series( x, (x,-pi,pi) )
s
Out[8]:
2sin(x)sin(2x)+2sin(3x)3+
In [9]:
plot( s.truncate(10), (x,-2*pi,2*pi) )
Out[9]:
<sympy.plotting.plot.Plot at 0x222f79b4c18>

퓨리어 변환

함수 f(x) 의 퓨리어 변환 ˆf(w)는 다음과 같이 정의된다. [Kreyszig]

ˆf(w)=12πf(x)eiwxdxf(x)=12πˆf(w)eiwxdw

SymPy에서 퓨리어 변환 F(k)는 다음과 같이 정의된다. [Wikipedia]

F(k)=f(x)e2πikxdxf(x)=F(k)e2πikxdk

퓨리어 변환에 관한 표현들이 서로 조금 다름에 유의하고, 두 정의를 비교하면 다음과 같다.

w=2πk

ˆf(w)=12πF(k)

In [10]:
k, w = symbols('k w')
a = Symbol('a', positive=True)

Gaussian 함수 eax2a>0

In [11]:
Fk = fourier_transform( exp(-a*x**2), x, k )
Fk
Out[11]:
πeπ2k2aa

퓨리어 변환 ˆf(w) 로 나타내면

In [12]:
fw = Fk.subs( k, w/(2*pi) ) / sqrt(2*pi)
fw
Out[12]:
2ew24a2a

지수함수 (좌우대칭) ea|x|a>0

In [13]:
Fk = fourier_transform( exp(-a*abs(x)), x, k )
Fk
Out[13]:
2aa2+4π2k2
In [14]:
fw = Fk.subs( k, w/(2*pi) ) / sqrt(2*pi)
fw
Out[14]:
2aπ(a2+w2)

지수함수 (양의 축) eaxθ(x)a>0,x>0

In [15]:
Fk = fourier_transform( exp(-a*x)*Heaviside(x), x, k )
Fk
Out[15]:
1a+2iπk
In [16]:
fw = Fk.subs( k, w/(2*pi) ) / sqrt(2*pi)
fw
Out[16]:
22π(a+iw)

델타함수 δ(x)

In [17]:
Fk = fourier_transform( DiracDelta(x), x, k )
Fk
Out[17]:
1
In [18]:
fw = Fk.subs( k, w/(2*pi) ) / sqrt(2*pi)
fw
Out[18]:
22π
In [19]:
fourier_transform( 1, x, k )
Out[19]:
Fx[1](k)

상수 1의 퓨리어 변환은 델타함수 δ(k)이다.

SymPy에서 현재 지원되지 않고 있다.

[참고] 다른 형태의 표기법도 흔히 사용된다.

함수 f(x) 의 퓨리어 변환은 다음과 같이 정의된다. [Boas]

F(α)=12πf(x)eiαxdxf(x)=F(α)eiαxdα

함수 f(x) 의 퓨리어 변환은 다음과 같이 정의된다. [Arfken]

g(ω)=12πf(x)eiωxdxf(x)=12πg(ω)eiωxdω

Numeric Analysis 4 - Numeric Linear Algebra

Numeric Analysis 4 - Numeric Linear Algebra Numeric Linear Algebra ¶ ...