การใช้งานเลขเชิงซ้อนด้วย python/sympy

เลขเชิงซ้อนหรือ complex number เป็นเรื่องที่พบได้ทั่วไปในการคำนวณทางคณิตศาสตร์และวิทยาศาสตร์ เราเรียนเลขเชิงซ้อนมาตั้งแต่ชั้นมัธยมซึ่งก็ไม่ใช่เรื่องยากอะไรคือเลขที่เป็นผลบวกของเลขจริงกับเลขที่เป็นจำนวนเท่าของเลข i หรือ เช่น 10.2+6i

Python สามารถคำนวณเลขเชิงซ้อนได้โดยตรงโดยไม่ต้องใช้แพคเกจใดๆ โดยใช้สัญลักษณ์ เช่น

expr1 = 1+2j
expr2 = 3+4j
print(exp1*expr2)

ซึ่งเอาท์พุทคือ (-5+10j)

1. การคูณเลขเชิงซ้อนด้วย sympy

แต่ในงานทางด้านscience/engineering เราไม่ได้เกี่ยวข้องกับการคำนวณแบบเครื่องคิดเลขอย่างนี้แต่จะเ­กี่ยวข้องกับสมการที่เต็มไปด้วยตัวแปรที่ต้องมีการจัดรูปหรือแก้สมการออกมาให้ติดในรูปตัวแปรแล้วนำสมการมาเขียนกราฟหรือคำนวณเป็นตัวเลขจริงในขั้นตอนสุดท้ายโดยการแทนค่าตัวแปร

ใน sympyจะแทนเลขเชิงซ้อน $\sqrt{-1}$ ด้วย I (อักษร “ไอ” ตัวใหญ่) ซึ่งจะต้อง import เข้ามา การคำนวณข้างต้นหากเขียนในสไตล์ของ sympy จะเป็น

from sympy import I
expr1 = 1 + 2*I
expr2 = 2 + 3*I
expr = expr1*expr2
print(expr)

คำสั่งที่ 2 และ 3 เราสร้างสัญลักษณ์ expr1 และ expr2 แทนเลขเชิงซ้อน 1+2i และ 2+3i ตามลำดับ คำสั่ง 4 เอาเลขสองตัวมาคูณกันซึ่งเอาท์พุทก็คือ (1+2i)(2+3i) เหมือนไม่ได้มีการคูณกันทั้งนี้เพราะว่า sympy มองทุกอย่างเป็นสัญลักษณ์ หากต้องการคำนวณเป็นตัวเลขสุดท้ายจะต้องใช้คำสั่ง evalf

from sympy import I, evalf, simplify
eq1 = 1 + 2*I
eq2 = 2 + 3*I
eq = eq1*eq2
eq.evalf()
eq.simplify()

คำสั่ง evalf() จะให้คำตอบเป็น -4.0+7.0i ซึ่งหากไม่ระบุอะไรชัดเจน sympy จะมองเลขจำนวนเต็มเป็นเลขทศนิยม เช่น 1+ 2*I จะมองเป็น 1.0 + 2.0*I เมื่อหาผลคูณก็จะให้คำตอบเป็นเลขทศนิยมด้วย หากเราต้องการคำตอบที่ไม่มี .0 ก็ใช้คำสั่ง simplify เข้าช่วย เอาท์พุทจะเป็น -4+7i หรือวิธีที่ง่ายกว่าคือใช้คำสั่ง eq.expand() เพียงคำสั่งเดียวแทนที่คำสั่ง eq.evalf() และ eq.simplify()

เนื่องจาก sympy เป็นเครื่องมือในการจัดการกับสัญลักษณ์ การคูณเลขเชิงซ้อนจึงเขียนในสไตล์ของ sympy ได้เป็น

from sympy import symbols, I
a, b, c, d = symbols('a b c d')
exp = (a+b*I)*(c+d*I)
print(exp)

เอาท์พุทคือ $(a+bi)(c+di)$ ซึ่งยังคงติดในรูปสัญลักษณ์หากต้องการกระจายเทอมจะต้องใช้ expr.expand() จึงจะได้เอาท์พุทเป็น

ac + iad + ibc - bd

ซึ่งจะไม่มีการรวมเทอมที่เป็น i ให้เราต้องใช้คำสั่ง factor(I) ต่อท้ายเพื่อให้จัดองค์ประกอบตาม I คำสั่งจึงต้องเป็น expr.expand().factor(I) เอาท์พุทจะเป็น

ac-bd+i(ad+bc)

2. การหารเลขเชิงซ้อนด้วย sympy

การหารมีขั้นตอนมากกว่าคูณ อาศัยสูตรการหารเลขเชิงซ้อน

\frac{a+bi}{c+di} = \frac{(a+bi)(c-di)}{c^2+d^2} = \frac{ac+bd-i(ad-bc)}{c^2+d^2}

โปรแกรมจะเป็น

from sympy import symbols, I
a, b, c, d = symbols('a b c d')
expr1 = a + I*b
expr2 = c + I*d
number = ( (expr1*expr2).conjugate() ).expand().factor(I)
denom = abs(expr2)**2
print(numer/denom)

ซึ่งเอาท์พุทคือ $\frac{ac+bd-i(ad-bc)}{c^2+d^2}$

expr2.conjugate() หมายถึงให้หา complex conjugate ของ expr2,
abs(expr2) หมายถึงขนาดของ expr2 ซึ่งเท่ากับ $\sqrt{c^2+d^2}$
abs(expr2)**2 หมายถึงเอามายกกำลังสองอีกที ซึ่งเท่ากับ $c^2+d^2$
หาก expr1 = 1+2i และ expr2 = 3+4i ก็จะใช้คำสั่ง

(numer/denom).subs({a:1, b:2, c:3, d:4})

เอาท์พุทคือ $\frac{11}{25} + \frac{2}{25}i$ แต่หากต้องการผลหารสุทธิก็ใช้คำสั่ง

(numer/denom).subs({a:1, b:2, c:3, d:4}).evalf()

เอาท์พุทจะเป็น $0.44+0.08i$

ทำไมต้องทำให้ยุ่งยาก? เขียนแค่นี้ก็ได้ผลเมือนกัน

from sympy import I
expr1 = 1 + 2*I
expr2 = 2 + 3*I
expr = expr1/expr2
print(expr)

เอาท์พุทคือ $0.44+0.08j$

ขอให้เข้าใจว่า sympy เล่นกับ symbol หรือที่เรียกว่า derivation หรือการ derive สมการ ซึ่งเกี่ยวข้องกับติดค่าเป็นตัวแปรไม่ใช่การคำนวณโดยตรง เรากำลังเกี่ยวข้องกับสูตรที่การแทนค่าตัวแปรเป็นขั้นตอนสุดท้าย

3. เลขเชิงซ้อนในรูป exponential หรือ Euler Identity

e^\theta = \cos\theta + i\sin\theta

ทางซ้ายบางทีเรียกว่า polar form ทางขวาเรียกว่า rectangular form

from sympy import exp, I
theta = symbols('theta', real=True)
eq = exp(theta*I)
print(eq)

ซึ่งเอาท์พุทคือ $e^{i\theta}$ หากต้องการกระจายเป็น rectangular form ก็ใช้คำสั่ง expand(complex=True)

from sympy import exp, I
theta = symbols('theta', real=True)
eq = exp(theta*I)
eq.expand(complex=True)

เอาท์พุทคือ $i\sin\theta + \cos\theta$

คำสั่งที่ 2 ต้องใช้ real=True กำกับเพื่อบอกว่าสัญลักษณ์ เป็นเลขจำนวนจริง ปกติหากเอาสัญลักษณ์ที่นิยามโดยคำสั่ง symbol มาคำนวณร่วมกับตัวเลขเชิงซ้อนที่มี i อยู่ด้วย sympy จะถือว่าสัญลักษณ์นั้นเป็นเลขเชิงซ้อนไปด้วย คำสั่ง expand(complex=True) จะกระจายเลขเชิงซ้อนที่อยู่ในรูป $e^{i\theta}$ ให้อยู่ในรูปผลบวก real + i*real

ลองเอา real=True ในคำสั่งที่ 2 ออกจะได้เอาท์พุทประหลาดๆ

ie^{-im(\theta)}\sin(re(\theta))+e^{-im(\theta)}\cos(re(\theta))

ทั้งนี้ก็เพราะว่า $\theta$ ถูกมองว่าเป็นเลขเชิงซ้อนเมื่อเราเอามา expand

แต่หากเราไม่ใช้ symbol หรือใช้ตัวแปรก็ไม่มีปัญหาอะไรเช่น

from sympy import exp
eq = exp(2*I)
eq.expand(complex=True)

เอาท์พุทคือ $\cos(2)+i\sin(2)$

4. การแปลงจาก rectangular form ไปเป็น polar form

Sympy ไม่มีคำสั่งในการแปลงระหว่าง polar form และ rectangular form โดยตรง ต้องเขียนเป็นฟังก์ชั่นช่วย เช่น 1+i เมื่อเปลี่ยนไปเป็น polar form จะได้

1+i = \sqrt{2}e^{i\frac{\pi}{4}}

$\sqrt{2}$ คือขนาดของเลขเชิงซ้อนซึ่งก็คือ $\sqrt{1^2+1^2}$ และ $\frac{\pi}{4}$ คือมุมระหว่าง real part และ imaginary part ซึ่งก็คือ $arctan\frac{imaginary-part}{real-part}$ ซึ่งมุมจริงจะต้องมาบวกลบ $\pi$ หรือลบ $2\pi$ ตามแต่ตำแหน่ง quandrant

สูตรการแปลงจาก ไปเป็น polar form คือ

(\sqrt{a^2+b^2})e^{i\theta}

โดย $\theta = tan^{-1}\frac{b}{a} = arctan\frac{b}{a}$ แต่ว่า $\theta$ เปลี่ยนไปตาม quadrant ตามรูป

เขียนเป็นฟังก์ชั่น toPolar() เพื่อเปลี่ยนเลขเชิงซ้อนเป็น polar form

from sympy import exp, I, atan, re, im, pi

def toPolar(expr):
     
    #case re(expr)>0 and im(expr)>0
    angle = atan( im(expr)/re(expr) )
    theta = angle
    
    if re(expr)<0 and im(expr)>0:
        theta = pi + angle
        
    if re(expr)<0 and im(expr)<0:
        theta = pi + angle
        
    if re(expr)>0 and im(expr)<0:
        theta = 2*pi + angle
         
    return abs(expr)*exp(I*theta)
        
def toRectangular(expr):
    return expr.expand(complex=True).factor(I)

expr = -1+1*I
expr = toPolar(expr)
#toRectangular(expr)เป็นการแปลง -1+i ไปเป็น polar form ซึ่งเอาท์พุทคือ
from sympy import exp, I, atan, re, im, pi

def toPolar(expr):
     
    #case re(expr)>0 and im(expr)>0
    angle = atan( im(expr)/re(expr) )
    theta = angle
    
    if re(expr)<0 and im(expr)>0:
        theta = pi + angle
        
    if re(expr)<0 and im(expr)<0:
        theta = pi + angle
        
    if re(expr)>0 and im(expr)<0:
        theta = 2*pi + angle
         
    return abs(expr)*exp(I*theta)
        
def toRectangular(expr):
    return expr.expand(complex=True).factor(I)

expr = -1+1*I
expr = toPolar(expr)
#toRectangular(expr)

คำสั่งที่ใช้คือ re, im และ abs ซึ่งใช้หา real part, imaginary part และขนาดของเลขเชิงซ้อนตามลำดับ

เป็นการแปลง -1+i ไปเป็น polar form ซึ่งเอาท์พุทคือ $\sqrt{2}e^{\frac{3i\pi}{4}}$

ลองเปลี่ยนตัวเลขในคำสั่ง 23 ดู และ run คำสั่ง toPolar() เพื่อแปลงไปเป็น polar form คำสั่ง toRectangular() ใช้เปลี่ยน polar form ไปเป็น rectangular form อีกที ซึ่งในโปรแกรมใช้สำหรับ check กลับ เอาท์พุทก็จะควรจะได้ตรงกับเลขเลขซ้อนที่ใส่ในคำสั่ง 23

4. Note

1. Sympy ใช้ atan แทน arctan
2. ลองแปลง toPolar(-1-1*I) และ toPolar(1-1*I) ซึ่งเป็น complex number ใน quadrant 3และ 4 สังเกตว่าเอาท์พุทคือ $\sqrt{2}e^{-\frac{3i\pi}{4}}$ และ $\sqrt{2}e^{-\frac{i\pi}{4}}$ ตามลำดับ
3. สังเกตว่าใน quadrant 3 และ 4 sympy จะเปลี่ยน ให้เป็นมุมติดลบซึ่งเทียบเท่ากับมุมบวก $\frac{5\pi}{4}$ และ $\frac{7\pi}{4}$ ตามลำดับ

จรัสพรรณ เปรมปรีบุตร

แต่ก่อนไม่ชอบ python โตมากับ C/Java ตอนนี้หลงรัก มีเล่าเรื่องฟิสิกส์ เรื่องจักรวาลในมุมมองแบบชาวบ้านประกอบสำหรับคนชอบ