การใช้ Sympy/Python เพื่อแก้สมการและปัญหาทางคณิตศาสตร์

Sympy เป็น package ที่เกี่ยวกับสัญลักษณ์ทางคณิตศาสตร์ ใช้ในการจัดรูปสมการ, กระจายเทอม, การจัด factor, การแก้สมการและการคำนวณทาง calculus

การใช้งานเพียงแต่ import sympy เข้ามา ซึ่งทำได้สามวิธีคือ

1. import sympy as sp จากนั้นทุกๆคำสั่งของ sympy ให้ขึ้นต้นว่า sy.คำสั่งที่ต้องการ วิธีนี้จะไม่มีปัญหาเมื่อใช้คำสั่งที่มีชื่อเหมือนกันจากแพคเกจที่ต่างกัน เป็นวิธีที่ถูกต้องที่สุด

2. from sympy import * วิธีนี้ง่ายที่สุดทำให้เราพิมพ์คำสั่งของ sympy ได้โดยตรงโดยไม่ต้องพิมพ์ sympy. หรือ sy. นำหน้า แต่จะมีปัญหาหาก import หลาย package ที่มีชื่อคำสั่งตรงกันจึงไม่ควรใช้และโปรแกรมจะแจ้งเตือน แต่ถ้าคุณต้องการทดสอบ code สั้นๆเฉพาะกิจก็ใช้ได้

3. from sympy import ชื่อคำสั่ง1, ชื่อคำสั่ง2 เหมือนวิธีที่สองแต่จะ import เฉพาะคำสั่งที่ต้องการใช้โอกาสที่จะสับสนกับแพคเกจอื่นที่มีชื่อคำสั่งเดียวกันจะน้อยลง

เนื่องจากในตอนนี้เรากำลังใช้ sympy ตัวเดียวซึ่งวิธีที่ง่ายที่สุดคือแบบที่สอง แต่ตัวอย่างในบทความนี้จะใช้แบบที่สามเพื่อให้คุณเห็นที่หัวว่าต้องใช้คำสั่งใดบ้าง

วิธีที่ง่ายที่สุดคือทดลองบน jupyter notebook ซึ่งจะให้เอาท์พุทเป็นสมการที่สวยงามอยู่แล้ว คุณเพียงแต่ copy/paste โปรแกรมตามตัวอย่างแล้ว run ดูผลและทดลองเปลี่ยนสมการดูโดยเปลี่ยนคำสั่งแสดงผลจาก print(xxx) เป็น xxx แต่หากคุณรันบน console ก็ไม่ต้องเปลี่ยนอะไรเพียงแต่สมการจะดูไม่สวยงามเท่านั้น หากคุณรันบน spyder ให้ดูตอนท้ายของบทความจะมีวิธีแปลงเอาท์พุทให้อยู่ในรูปสมการที่สวยงามให้

1. สมการตัวแปรเดียว

สมมติว่าเราต้องการป้อนสมการนี้เข้าไปผ่าน sympy

x^2+6x+8
from sympy import Symbol
x = Symbol('x')
expr = x**2 + 6*x + 8
print(expr)

หากรันใน console เอาท์พุทคือ x**2 + 6*x + 8

2. สมการหลายตัวแปร

หากมีตัวแปรหลายตัวให้ใช้คำสั่ง symbols แทนจะสะดวกกว่าที่ต้องมาใช้ Symbol ทีละบรรทัด

from sympy import symbols
x, y = symbols('x y')
expr = x**2 + 6*x + y
print(expr)

$x^2 + 6x + y$

มีตัวแปรเดียวก็ใช้ symbols ได้เช่น x = symbols(‘x’) ก็เทียบเท่ากับ x = Symbol(‘x’) ส่วนที่อยู่ใน argument หมายถึงสัญลักษณ์ที่เราต้องการให้แสดงที่เอาท์พุทซึ่งสามารถใช้อักษรกรีกได้ ตัวแปรทางซ้ายหมายถึงที่ใช้อ้างในสมการ เช่น

from sympy import symbols, sin, pi
om, ld, t = symbols('ω λ t')
expr = (2*pi/ld)*sin(om*t)
print(expr)

$\frac{2 \pi \sin{\left(t ω \right)}}{λ}$

3. สมการซ้อนสมการ

หากเรามีสมการยาวๆหรือซับซ้อนเราสามารถแยกเป็นสมการย่อยๆแล้วใช้วิธีแทนที่ตัวแปร เช่น

from sympy import symbols, sqrt
x, y = symbols('x y')
y = x**2 
y = x/(x+1)*sqrt(y)
print(y)

เราต้องการสมการ $\frac{x \sqrt{x^2}}{x+1}$ แต่เราไม่ใส่เข้าไปตรงๆแต่แยกเป็น y=x2 ก่อนจากนั้นก็แทนลงในสมการ $\frac{x \sqrt{y}}{x+1}$

4. Sympy มีการเรียงลำดับตัวแปรใหม่

sympy มีข้อเสียอย่างหนึ่งตรงที่สมการที่พิมพ์ออกมาจะเรียงตามตัวอักษรของตัวแปร ตามกำลังสูงไปต่ำและตามด้วยตัวเลข เช่น

from sympy import symbols
x, b = symbols('x b')
expr = x**2 + 6*x + 5
print(expr)
expr = x**2 + 6*x + b
print(expr)

เอาท์พุทแรกคือ x**2 + 6*x + 5 ซึ่งตรงกับที่เราคาดหวัง แต่เอาท์พุทต่อไปคือ b + x**2 + 6*x ซึ่งผิดไปจากที่เราคาดหวังคือ x**2 + 6*x + b ทั้งนี้เพราะเทอมในสมการจะถูกจัดเรียงใหม่ตามตัวอักษรเพราะ b อยู่ก่อน x เลขยกกำลังของ x เรียงจากสูงไปต่ำจาก 2 ไป 1 ตามเดิม

ตัวอย่างของรูปสมการที่เปลี่ยนไปอีกอันหนึ่ง

from sympy import symbols
x, a = symbols('x a')
expr = (x - a)**2
print(expr)

เอาท์พุทคือ (-a + x)**2 สมการเช่นนี้เรามักจะใช้ในกรณีที่เราต้องการให้ a เป็นตัวแปรไว้ก่อนจากนั้นค่อยมาเปลี่ยน a ทีหลัง แต่ถ้าเราไม่ print สมการออกมาก็จะไม่ใช่ปัญหาอะไร

5. การใช้ subs แทนที่ตัวแปร

5.1 ใช้แทนที่ตัวแปร

หากเรามีสมการเดิมอยู่แล้วและต้องการเปลี่ยนตัวแปรใหม่ก็ทำได้ด้วยคำสั่ง subs (substitute=แทนที่)

from sympy import symbols
x, y = symbols('x y')
expr = x**2 + 6*x + 8
expr = eq.subs(x, y)
print(expr)

ซึ่งเอาท์พุทก็คือ y**2 + 6*y + 8

การแทนที่ไม่จำเป็นต้องเป็นตัวแปรที่นิยามด้วย Symbol หรือ symbols อาจเป็นฟังก์ชั่นคณิตศาสตร์ก็ได้

from sympy import symbols, sin
x, y = symbols('x y')
expr = x**2 + 6*x + 8
expr = expr.subs(x, sin(y))
print(expr)

ซึ่งเอาท์พุทก็คือ $(\sin y)^2 + 6sin(y) + 8$

ข้อสังเกต

subs ไม่ใช่คำสั่งอิสระใน package ของ sympy โดยตรงจึงไม่ต้อง import เข้ามา ตามหลักของภาษาแบบ OOP (Object-Oriented Programming ซึ่ง python เป็นหนึ่งในนั้น) คำสั่งอะไรที่นำหน้าด้วยจุดเรียกว่า method เป็นคำสั่งที่ฝังอยู่ใน object ที่อยู่ก่อนหน้าจุด ในกรณีนี้ eq คือ object ที่เป็นตัวแทนของสมการซึ่งถูกสร้างขึ้นมาจาก object ชื่อ x และ y อีกที คำสั่ง subs(หรือจริงๆควรเรียก method) จะแทนที่ x หรือ y ตามที่เรากำหนด

5.2 ใช้แทนค่าตัวแปรด้วยตัวเลข

หากแทนที่ด้วยตัวเลขก็จะเป็นการแทนค่าลงไปในสมการเพื่อหาคำตอบ

from sympy import symbols
x = Symbol('x')
expr = x**2 + 6*x + 8
print( expr.subs(x, 1) )
print( expr )

เอาท์พุทคือ 15 และ x*2 + 6x*y + 8 ตามลำดับ สังเกตว่า expr ยังคงเป็นสมการเดิมหลังจาก subs ในคำสั่ง 4 เพราะการแทนที่เกิดขึ้นเฉพาะกิจไม่ได้ไปเปลี่ยนแปลงตัวสมการ หากต้องการเปลี่ยนสมการถาวรต้องเขียนเป็น

from sympy import Symbol
x = Symbol('x')
expr = x**2 + 6*x + 8
expr = expr.subs(x, 1)
print( expr )

หรือหากมีตัวแปรสองตัว

from sympy import symbols
x, y = symbols('x y')
expr = x**2 + 6*x*y + 8
print( expr.subs( {x:1, y:2} ) )
print( expr.subs( [ (x, 1), (y, 2) ] ) )
print(expr)

เอาท์พุทคือ 21 และ x*2 + 6x*y + 8 , โดยบรรทัดที่ 4 และ 5 เป็นการแทนค่า x=1 และ y=2 เหมือนกันให้เลือกใช้แบบใดแบบหนึ่ง

เช่นเดียวกับตัวอย่างก่อนหน้าเราแทนที่โดยไม่เอาผลใส่กลับเข้าไปใน expr คือไม่ได้ใช้คำสั่ง
expr = expr.subs( [ (x, 1), (y, 2) ] ดังนั้นสมการเดิมยังคงอยู่

ข้อสังเกต

subs จะคำนวณให้เฉพาะสมการ polynomial เท่านั้นหากในสมการมีฟังก์ชั่น sin, cos, exponential หรืออื่นๆจะไม่คำนวณให้ยังคงค้างเป็นชื่อฟังก์ชั่นเอาไว้ เช่น

from sympy import Symbol, sin, cos, exp
x = Symbol('x')
expr = x**2 + 2*sin(x)*cos(x) + exp(3*x)
print( expr.subs(x, 2) )

เอาท์พุทคือ 2sin(2)cos(2) + 4 + exp(6), สังเกตว่าคำนวณเฉพาะ x**2 = 4 (ซึ่งเป็น polynomial)นอกนั้นแทนค่า x=2 ค้างเอาไว้ ซึ่งหากต้องการคำนวณเป็นตัวเลขสุทธิต้องใช้ evalf หรือ lambdify

6. ใช้ evalf คำนวณ

ต่างจาก subs คำสั่งนี้จะคำนวณค่าจริงๆไม่ใช่เป็นการแทนที่ ใช้ได้กับทุกฟังก์ชั่นไม่จำเป็นต้องเป็น polynomial แต่ควรใช้เป็นคำสั่งเฉพาะกิจเมื่อต้องการคำนวณสั้นๆ

from sympy import sqrt, sin, pi, exp
print( sqrt(2).evalf() )
print( sin(pi).evalf() )
print( exp(2).evalf() )

คำตอบคือ 1.41421356237310 และ 0 และ e2 = 7.38905609893065

ปกติแสดง 15 ตำแหน่ง หากต้องการเปลี่ยนก็ให้ระบุเช่น

from sympy import sqrt
print( sqrt(2).evalf(6) )

คำตอบคือ 1.41421 แสดง 6 ตำแหน่ง

หากต้องการแทนตัวเลขลงในสมการ

from sympy import sqrt, symbols
x, y = symbols('x y')
expr = x + 1
print( expr.evalf( 5, subs={x:2} ) )
expr = sqrt(x + y)
print( expr.evalf( 5, subs={x:1, y:2} ) )

คำตอบคือ 3.0000 และ 1.7320 โดย 5 คือจำนวนตำแหน่ง

ข้อสังเกต

sympy มองทุกอย่างเป็น symbol หรือสัญญลักษณ์ ตามตัวอย่างนี้

from sympy import pi
print( 20/(2*pi) )
print( (20/(2*pi)).evalf() )

เอาท์พุทแรกจะเป็น 10/pi เพราะ sympy จะค้างสัญญลักษณ์เอาไว้ หากต้องการผลลัพท์ต้องใช้ evalf ในคำสั่งที่ 3

evalf เป็น method เช่นเดียวกับ subs จึงไม่ต้อง import เข้ามา

7. ใช้ lambdify คำนวณ

คำสั่งนี้มีประสิทธิภาพมากกว่า evalf เมื่อมีการคำนวณหลายขั้นตอนควรใช้คำสั่งนี้ เราต้องนิยามสมการให้อยู่ในรูปฟังก์ชั่นทางคณิตศาสตร์ก่อนโดยผ่านคำสั่ง lambdify

from sympy import Symbol, lambdify
x = Symbol('x')
expr=x**2
y = lambdify(x, expr)
print( y(3) )

คำสั่งที่ 4 คือการนิยามฟังก์ชั่น y(x) = x**2 โดย x หมายถึงโดเมนของฟังก์ชั่น เราต้องการหาค่า y เมื่อ x=3 ซึ่งก็คือคำสั่ง y(3) เอาท์พุทคือ 9

(หรือเขียนแบบรวบรัดยุบบรรทัดที่ 3, 4 เป็น y=lambdify( x, x**2) ) หรือเขียนอีกแบบ

from sympy import Symbol
x = Symbol('x')
y = lambda x:x**2
print( y(3) )

หากมีโดเมนมากกว่าหนึ่งเช่น x, y

from sympy import symbols, lambdify
x, y = symbols('x y')
eq = x**2 + 3*y
z = lambdify( [x, y], eq )
print( z(1, 2) )

หรือ

from sympy import symbols
x, y = symbols('x y')
z = lambda x, y:x**2 + 3*y
print( z(1, 2) )

เรานิยามฟังก์ชั่น z(x, y) = x2 + 3y แล้วคำนวณค่า z(1, 2) ซึ่งเอาท์พุทคือ 7

จากตัวอย่างสุดท้ายในหัวข้อ 5 ที่ค้างเป็นฟังก์ชั่นเพราะใช้คำสั่ง subs เมื่อเปลี่ยนมาคำนวณด้วย lambdify ก็จะเป็น

from sympy import Symbol, sin, cos, exp, lambdify
x = Symbol('x')
expr = x**2 + 2*sin(x)*cos(x) + exp(3*x)
f = lambdify(x, expr)
print( f(2) )

เอาท์พุทคือ 406.6719909974272, แต่ถ้าเขียนโดยใช้ lambda

from sympy import Symbol, sin, cos, exp
x = Symbol('x')
f = lambda x:x**2 + 2*sin(x)*cos(x) + exp(3*x)
f(2).evalf()

f(2) ยังคงติดในรูป symbols ต้องใช้คำสั่ง evalf ถึงจะได้คำตอบสุดท้าย

8. การแยก factor

Example 8.1

x^2+6x+8=(x+2)(x+4)

ทำโดยคำสั่ง factor

from sympy import Symbol
x = Symbol('x')
eq = x**2 + 6*x + 8
print( eq.factor() )

เอาท์พุทคือ (x + 2)*(x + 4)

Exercise

x^2 + x \xrightarrow[]{.factor()}x(x+1)

Example 8.2 ถ้าต้องการแยกแบบนี้ให้ใช้ .factor_terms()

-x-1 \xrightarrow[]{factor\_terms} -(x+1)
from sympy import Symbol, factor_terms
x = Symbol('x')
expr = -x-1
print( factor_terms(expr) )

Example 8.2

from sympy import symbols
x = symbols('x')
expr = x**(1/2) + x**(5/2)
print( expr.factor() )

หากเลขยกกำลังเป็นเศษส่วนหรือเลขทศนิยม factor จะไม่แยกองค์ประกอบ เอาท์พุทคือ $x^{0.5}+x^{2.5}$ ต้องใช้ Rational(1/2) และ Rational(5/2)

from sympy import symbols, Rational
x = symbols('x')
expr = x**(Rational(1/2)) + x**(Rational(5/2))
print( expr.factor() )

คำสั่งที่ 3 จะบอกให้ sympy มองเป็น $expr=\sqrt{x}+x^{5/2}$ ไม่ใช่ $x=x^{0.5}+x^{2.5}$ และจะแยกองค์ประกอบออกมาให้เป็น $\sqrt{x}(x^2+1)$

Example 8.3

(\cos x)^2-(\sin y)^2
from sympy import symbols, sin, cos
x, y = symbols('x y')
expr = cos(x)**2-sin(y)**2
print( expr.factor() )

เอาท์พุทคือ $(\cos x-\sin y)(\cos x+\sin y)$

9. การกระจายเทอม expand

หากต้องการกระจายเทอมก็ใช้คำสั่ง expand

Example 9.1

(x+2)(x+4)=x^2 + 6x +8
from sympy import Symbol
x = Symbol('x')
expr = (x + 2)*(x + 4)
print( expr.expand() )

เอาท์พุทคือ $x^2 + 6x + 8$

Example 9.2 เปลี่ยนเทอมหารเป็นเทอมบวก

\frac{x+y}{xy} = \frac{1}{x}  + \frac{1}{y}
from sympy import symbols
x, y= symbols('x y')
expr = (x+y)/(x*y)
print( expr.expand() )

เอาท์พุทคือ $\frac{1}{y} + \frac{1}{x}$

10. การทำให้สมการอยู่ในรูปที่ง่ายขึ้นด้วย simplify

คำสั่งนี้เป็นคำสั่งทั่วไปที่ช่วยลดรูปสมการให้ดูง่ายที่สุดเท่าที่ทำได้

Example 10.1

x^2+yx=x(x+y)
from sympy import symbols
x, y = symbols('x y')
expr = x**2 + y*x
print( expr.simplify() )

ในกรณีนี้ simplify มีผลเหมือนคำสั่ง factor คุณลองเปลี่ยน simplify เป็น factor ใน code ตัวอย่างดูก็จะให้ผลเหมือนกัน จุดต่างคือ simplify ทำให้ดูง่ายแต่ไม่ได้มีความสามารถเดียวกับคำสั่ง factor คุณลองเปลี่ยนบรรทัดที่ 3 เป็น x**2 + 6*x + 8 จะเห็นว่า simplify ไม่ได้แยกองค์ประกอบเป็น (x+2)(x+4)

Example 10.2

\frac{x^2+x-2}{x+2} = \frac{(x-1)(x+2)}{x+2} = x-1 
from sympy import symbols
x = symbols('x')
expr = (x**2+x-2)/(x+2)
print( expr.simplify() )

ผลคือ x-1 ลองสลับใช้ simplify กับ factor ดูจะได้เอาท์พุทเหมือนกัน แม้ว่า simplify จะไม่ได้แยกองค์ประกอบเทอมบนโดยตรงแต่เมื่อมีเทอมล่างที่ตรงกันก็ตัดกันไปให้ผลเหมือนกันตัวอย่างนี้แสดงให้เห็นว่าคำสั่ง factor ทำหน้าทั้งการแยกองค์ประกอบแล้ว simplify ไปในตัวคือตัดเทอมที่เหมือนกันออกไป

Example 10.3

\frac{1}{x} + \frac{1}{y} = \frac{x+y}{xy}
from sympy import symbols
x, y= symbols('x y')
expr = 1/x + 1/y
print( expr.simplify() )

เอาท์พุทคือ $\frac{x+y}{xy}$ สังเกตว่ารูปที่ simplify คือรูปที่เกิดจากการคูณไขว้ หากต้องการย้อนกลับก็ใช้คำสั่ง expand ซึ่งก็คือสมการเดียวกับในหัวข้อ 9 ที่ผ่านมา

Example 10.4

แยกตัวแปรร่วมแบบเดียวกับ factor ใน example 8.2 แต่ใช้ simplify แทน

from sympy import symbols
x = symbols('x')
x = x**(1/2) + x**(5/2)
print( simplify(x, rational=True) )

ไม่ว่าจะเขียนเลขยกกำลังจะเป็น 1/2, 5/2 หรือ 0.5, 2.5 เอาท์พุทก็คือ $\sqrt{x}(x^2+1)$ ถ้าไม่มี rational=True เอาท์พุทจะเป็น $x^{0.5} + x^{2.5}$

Example 10.5

Simplify ในฟังก์ชั่นตรีโกณ

from sympy import symbols, sin, cos
x = symbols('x')
expr = sin(x)/cos(x)
print( expr.simplify() )

เอาท์พุทคือ tan(x), ลองเปลี่ยนบรรทัดที่ 3 เป็น expr = sin(x)**2 + cos(x)**2 ดู เอาท์พุทคือ 1

ใช้ simplify แบบเร่งด่วน

ถ้าจุดประสงค์ของคุณคือต้องการทำให้สมการดูง่ายขึ้นเท่านั้นโดยไม่ต้องการคำนวณเป็นตัวเลข คุณไม่จำเป็นต้องใช้คำสั่ง Symbol หรือ symbols เพื่อสร้างสัญลักษณ์ x, y, z, … ขึ้นในสมการ เช่นเราต้องการลดรูป

2sin(x)cos(x)=sin(2x)

เขียนคำสั่งง่ายๆเพียงใส่สมการเข้าไปในคำสั่ง simplify โดยมี ‘ ‘ คร่อม

from sympy import simplify
expr = simplify('2*sin(x)*cos(x)')
print( expr )

เอาท์พุทคือ sin(2*x),

สังเกตว่า หนึ่ง เราไม่ต้องนิยาม Symbol(‘x’) , สอง เราใส่สมการเป็น stringเข้าไปตรงๆใน simplify ได้เลย, สาม เราใช้ฟังก์ชั่น sin, cos แต่เราไม่จำเป็นต้อง import เข้ามา

หากใช้วิธีนี้จะไม่สามารถใช้ subs หรือ lambdify เพื่อทำการคำนวณได้เพราะไม่มีการนิยามสมการขึ้นมาจริงๆ ดังนั้นคำสั่งต่อไปนี้จึงใช้ไม่ได้

from sympy import symbols, simplify
eq = simplify('x**2 + sin(x) + exp(x*y)')
#x, y = symbols('x y')
print( eq.subs( [(x,2), (y,3)] ) )

เอาท์พุทคือ NameError: name ‘x’ is not defined (จริงๆทั้ง x และ y), เพราะการแทนที่หรือคำนวณด้วย subs (หรือ lambdify) จะต้องมีการสร้างตัวแปรด้วยคำสั่ง Symbol หรือ symbols มาก่อน หากเอา # ในคำสั่งที่ 3 ออกก็จะคำนวณได้ตามปกติซึ่งเอาท์พุทคือ sin(2) + 4 + exp(6)

แต่ถ้าต้องการคำตอบเป็นตัวเลขสุทธิก็ใช้ lambdify

from sympy import symbols, simplify, lambdify
eq = simplify('x**2 + sin(x) + exp(x*y)')
x, y = symbols('x y')
f = lambdify( [x, y], eq)
print( f(2, 3) )

คำตอบคือ 408.3380909195608

simplify ไม่ได้ให้ผลที่ตรงกับที่เราคาดหมายเสมอไปตัวอย่างเช่น

from sympy import symbols
x, m = symbols('x m')
expr = 2*x**(-m)/x
expr.simplify()

ผลที่คาดหวังคือ $2x^{-m-1}$ แต่คำตอบคือ $\frac{2 x^{- m}}{x}$ คือไม่ได้ทำอะไร ยกเว้นตัวคูณเป็น 1 ไม่ใช่ 2

Exercise ลอง simplify

\frac{x^{-m}}{x} \xrightarrow[]{simplify} x^{-m-1}

11. การจัดการกับฟังก์ชั่น polynomial

นอกจาก factor, simplify, expand ซึ่งเป็นคำสั่งพื้นๆแล้วยังมีคำสั่งที่ใช้จัดการกับสมการแบบ polymomial เพิ่มเติมอีกเล็กน้อย

11.1 แยก rational function ด้วย apart

\frac{x-1}{x^2+5x+6} = \frac{4}{x+3} + \frac{3}{x+2}
from sympy import apart
x, y = sympy.symbols('x y')
eq = apart( (x-1)/(x**2 + 5*x + 6))
print(eq)

Exercise ลองใช้ apart แยก

\frac{x+2}{x+1} = \frac{x+1+1}{x+1}=1 + \frac{1}{x+1}

11.2 ตัดเทอมเศษและส่วนที่เหมือนกันด้วย cancel

\frac{x+2}{x^2+5x+6} = \frac{x+2}{(x+2)(x+3)} = \frac{1}{x+3}
frpm sympy import cancel
x, y = sympy.symbols('x y')
eq = cancel( (x+2)/(x**2 + 5*x + 6))
print(eq)

ซึ่งให้ผลเช่นเดียวกับ simplify

Exercise ลองใช้ cancel กับ expression ต่อไปนี้

\frac{(4x-1)e^{4x}}{4} \xrightarrow[]{cancel} x e^{4 x} - \frac{e^{4 x}}{4}

สังเกตว่า cancel จะกระจายเทอมให้ด้วย

11.3 ยุบเป็นเทอมเดียวด้วย together

\frac{x}{y+1} + \frac{y}{x+1} = \frac{x(x+1) + y(y+1)}{(x+1)(y+1)}
from sympy import symbols, together
x, y = symbols('x y')
eq = together( x/(y+1) + y/(x+1) )
print(eq)

แต่ถ้าเราต้องการกระจายเทอมเฉพาะตัวตั้งให้เป็น

\frac{x}{y+1} + \frac{y}{x+1} = \frac{x^2 + y^2 + x + y}{(x+1)(y+1)}
from sympy import symbols, together, fraction
x, y = symbols('x y')
eq = together( x/(y+1) + y/(x+1) )
n, d = fraction(eq)
print( expand(n)/d )

fraction ในคำสั่งที่ 4 จะแยกเทอมตั้งและเทอมหารออกมาใส่ตัวแปร n และ d ตามลำดับ จากนั้นเราใช้ expand เฉพาะเทอมตั้ง

11.4 จัดการกับฟังก์ชั่น polynomial ด้วยแพคเกจของ numpy*

ที่ผ่านมาเราใช้คำสั่ง symbol สร้างตัวแปรแล้วนำมาสร้างสมการ แต่บางครั้งหากเราทำงานกับฟังก์ชั่น polynomial อย่างเดียวมีอีกทางเลือกคือใช้แพคเกจ numpy แทน

from numpy import poly1d
eq = poly1d([2, 1, 5, 6])
print( eq )
print( eq(1) )

ผลที่ได้คือ $eq = 2x^3 + x^2 + 5x +6$ เมื่อแทน x=1 คำตอบคือ 14

แต่ถ้าหากต้องการผลลัพท์อย่างเดียวไม่ต้องการสมการ

from numpy import polyval
eq = polyval([2, 1, 5, 6], 1)
print( eq )

คำตอบคือ 14 เช่นเดียวกัน

numpy มีคำสั่งอื่นๆอีกที่เกี่ยวข้องกับ polynomial โดยเฉพาะ เอาไว้โอกาสต่อไปหรือคุณไปค้นคว้าเพิ่มเติมเอง

12. เอกลักษณ์ทางตรีโกณมิติ

คำสั่ง expand ที่กล่าวมาข้างต้นใช้ได้ดีกับสมการแบบ polynomial หากเป็นฟังก์ตรีโกณจะมีคำสั่งเฉพาะ

12.1 กระจายด้วย expand_trig

sympy กระจาย identity (เอกลักษณ์) ของฟังก์ชั่นตรีโกณ(trigonometry) ด้วยคำสั่ง expand_trig

from sympy import symbols, sin, expand_trig
x, y = symbols('x y')
eq = sin(2 * x)
eq = expand_trig(eq)
print(eq)

เอาท์พุทคือ 2*sin(x)*cos(x)

12.2 ลดรูปด้วย trigsimp

ในทางกลับกันเราสามารถลดรูปหรือทำให้สมการง่ายลงด้วยคำสั่ง trigsimp (trigonometry simplify) ซึ่งเมื่อต้องการทำกลับจากตัวอย่างข้างต้น

from sympy import Symbol, sin, cos, trigsimp
x = Symbol('x')
eq = 2*sin(x)*cos(x)
eq = trigsimp(eq)
print(eq)

เอาท์พุทคือ sin(2*x), สังเกตว่าใช้ simplify ตามตัวอย่างก่อนหน้าได้เหมือนกัน

อีกตัวอย่างหนึ่ง เอกลักษณ์ทางตรีโกณพื้นฐานคือ

cos^2(x)+sin^2(x)=1
from sympy import Symbol, sin, cos, trigsimp
x = Symbol('x')
eq = cos(x)**2 + sin(x)**2
eq = trigsimp(eq)
print(eq)

เอาท์พุทคือ 1, คำสั่งนี้ดูเหมือนทำการคำนวณค่า $cos^2(x)+sin^2(x)$ จริงๆไม่ใช่แต่เป็นการลดรูปเทอมดังกล่าวลงให้อยู่ในรูปที่ง่ายที่สุด(เท่าที่ sympy จะทำได้) และมีผลเหมือนกับ simplify ในตัวอย่างก่อนหน้า

ในทางกลับกันคุณอาจคาดหวังว่า expand_trig จะกระจาย 1 กลับไปเป็น $cos^2(x)+sin^2(x)$

from sympy import Symbol, expand_trig
x = Symbol('x')
eq = 1
eq = expand_trig(eq)
print(eq)

แต่เอาท์พุทคือ 1, นั่นคือ sympy ไม่คิดว่า 1 จำเป็นจะต้องเป็น $cos^2(x)+sin^2(x)$ เสมอไป จึงไม่ทำอะไร

12.3 ใช่ว่าจะได้ตามที่ต้องการ

ลองเล่นกับเอกลักษณ์ที่ซับซ้อนขึ้น

\sin(\alpha+\beta) = \sin(\alpha)\cos(\beta) + \cos(\alpha)\sin(\beta)
\cos(\alpha+\beta)=\cos(\alpha)\cos(\beta) - \sin(\alpha)\sin(\beta)
from sympy import symbols, sin, cos, trigsimp
a, b = symbols('a b')
eq = sin(a)*cos(b)+cos(a)*sin(b) + cos(a)*cos(b)-sin(a)*sin(b)
eq = trigsimp(eq)
print(eq)

เอาท์พุทที่เราคาดหวังคือ sin(a + b)+cos(a+b) แต่ผลที่ได้คือ sqrt(2)*sin(a + b + pi/4) หรือ

\sqrt{2}sin(a+b+\frac{\pi}{4})

ซึ่งเท่ากันแต่ไม่ตรงตามที่คาดหวัง แต่หากเปลี่ยนโปรแกรมเป็น

from sympy import trigsimp
eq = trigsimp('sin(a)*cos(b)+cos(a)*sin(b)') + trigsimp('cos(a)*cos(b)-sin(a)*sin(b)')
print(eq)

เอาท์พุทคือ

\sin(a+b) + \cos(a+b)

สังเกตว่าเราใส่สมการเข้าไปเป็นพาราเมตเตอร์ใน trigsimp ได้โดยตรงเหมือน simplify ซึ่งคำสั่ง symbol ก็ไม่มีความจำเป็นและไม่จำเป็นต้อง import symbols, sin, cos

13. การใช้ sympy แก้สมการ

13.1 สมการ Quadratic/Polynomial

สมมติเราต้องการแก้สมการ

x^2 = 4

เรารู้ว่าคำตอบคือ x = +2 หรือ -2 เมื่อต้องการให้ sympy แก้สมการให้ย้ายข้างให้เป็นศูนย์ด้านใดด้านหนึ่ง

x^2-4=0
from sympy import Symbol, solve
x = Symbol('x')
result = solve(x**2 - 4, x)
print(result)

เอาท์พุทคือ [-2, 2] (ต้องทำความเข้าใจว่าเครื่องหมาย [ ] ไม่ใช่หมายถึงช่วงปิดในทางคณิตศาสตร์แต่เป็นเครื่องหมายของ list ใน python) อ่านว่าคำตอบคือ -2 หรือ 2

อีกตัวอย่าง

x^2 + 5x + 6 = 0
from sympy import Symbol, solve
x = Symbol('x')
result = solve(x**2 + 5*x + 6, x)
print(result)

คำตอบคือ [-3, -2] หมายถึง -3 หรือ -2 อยู่ในรูป list หรือวิธีที่ง่ายกว่าคือใช้ numpy แล้วใส่เฉพาะ coefficient*

from numpy import roots
print( roots([1, 5, 6]) )

คำตอบคือ [-3. -2.] ซึ่งเป็นข้อมูลประเภท numpy.ndarray ของ numpy ซึ่งไม่ใช่ list อาจจะดูงงๆหน่อย

อีกตัวอย่างกรณีที่คำตอบเป็นเลข complex

x^2 + 4x + 5 =0
from sympy import Symbol, solve
x = Symbol('x')
result = solve(x**2 + 4*x + 5, x)
print(result)

[ -2-I, -2+I ] โดย I หมายถึง i ที่เป็น complex number หรือใช้ numpy* จะง่ายกว่า

from numpy import roots
roots([1, 4, 5])

คำตอบคือ array( [ -2.+1.j, -2.-1.j ] ) โดย 1.j หมายถึง i

สมการ polynomial ที่สูงกว่ากำลังสองก็ทำแบบเดียวกัน

13.2 ฟังก์ชั่นอื่นๆ

ฟังก์ชั่นอื่นเช่น sin, cos, exp, log ก็ใช้คำสั่ง solve เช่นเดียวกัน

13.2.1 ฟังก์ชั่นตรีโกณ

sympy จะมองโดเมนของฟังก์ชั่นตรีโกณในหน่วยของ $\pi$

from sympy import sin, pi
sin(pi/2).evalf()

คำตอบคือ $sin(\pi / 2) = 1$ และลองเล่นดูคุณจะได้ $sin(0)=0, sin(\pi)=0, sin(3\pi/2)=-1$ (คำสั่ง evalf ดูหัวข้อ 6)

หรือเขียนอีกแบบ

from sympy import symbols, sin, pi
x = symbols('x')
eq = sin(x)
print( sin(x).evalf( subs={x:pi/2} ) )

13.3 อสมการ

ตัวอย่าง 1

x^2-4 < 0

มีคำตอบเป็น

(x-2)(x+2)<0
(x-2)<0 \;\; and \;\; (x+2)>0 \;\;=>\;\; x<2 \;\; and \;\; x>-2\\
or \\
(x-2)>0 \;\; and \;\;(x+2)<0 \;\;=>\;\; x>2 \;\;and\;\;x<-2

คำตอบคือช่วงเปิด (-2, 2) หาคำตอบด้วยคำสั่ง solve_poly_inequality และ Poly

from sympy import Poly, solve_poly_inequality
result = solve_poly_inequality(Poly(x**2-4), '<')
print(result)

เอาท์พุทคือ [Interval.open(-2, 2)] หมายถึงช่วงเปิด (-2,2) หากใช้ตัวแปรอื่นที่ไม่ใช่ x เช่น y จะต้องเพิ่มเติมเล็กน้อย

from sympy import Symbol, Poly, solve_poly_inequality
y = Symbol('y')
result = solve_poly_inequality(Poly(y**2-4, y), '<')
print(result)

บางครั้งอาจจะเจอตัวอย่างที่มีเครื่องเครามากขึ้นตามเว็บไซต์

from sympy import Poly, Symbol, solve_poly_inequality
x = Symbol('x')
eq = x**2 - 4 < 0
lhs = eq.lhs
p = Poly(lhs, x)
rel = eq.rel_op
result = solve_poly_inequality(p, rel)
print( result )

lhs หมายถึง left hand side หรือทางซ้ายของสมการซึ่งก็คือ x**2 – 4 ส่วน rel_op หมายถึง operator ที่ใช้เปรียบเทียบซึ่งก็คือ < คำตอบก็เหมือนกัน

13.4 อสมการของ rational function

\frac{x-2}{x+2}<0

มีคำตอบเป็น

x-2<0 \;\;and\;\; x+2>0 \\
or \\
x-2>0 \;\;and\;\; x+2<0

ซึ่งก็คือช่วงเปิด (-2, 2) เหมือนตัวอย่างก่อนหน้า ใช้คำสัง solve_rational_inequalities และ Poly

from sympy import Poly, Symbol, solve_rational_inequalities
result = solve_rational_inequalities( [[((Poly(x-2), Poly(x+2)), '<')]] )
print( result )

เอาท์พุทคือ [Interval.open(-2, 2)] หมายถึงช่วงเปิด (-2,2) บางครั้งอาจเจอตัวอย่างที่ซับซ้อนกว่าแต่ก็ผลเหมือนกัน

from sympy import Poly, Symbol, solve_rational_inequalities
ieq = ((x-2)/(x+2)) < 0
lhs = ieq.lhs
numer, denom = lhs.as_numer_denom()
pn = Poly(numer)
pd = Poly(denom)
rel = ieq.rel_op
result = solve_rational_inequalities( [[((pn, pd), rel)]] )
print( result )

as_numer_denom หมายถึงแยกตัวตั้ง(numer) และตัวหาร(denom) ออกมาจากเทอมซ้ายของอสมการ(ที่แยกมาด้วย lhs แล้ว) จากนั้นนำมาสร้างเป็นสมการ polynomial ด้วยคำสั่ง Poly

13.5 อสมการที่เป็นฟังก์ชั่นอื่นๆ

e^x>1

คำตอบคือช่วงเปิด (0, ∞)

from sympy import Symbol, solve_univariate_inequality, exp
x = Symbol('x')
ieq = exp(x) > 1
result = solve_univariate_inequality(ieq, x, relational=False)
print( result )

เอาท์พุทคือ Interval.open(0, oo) หากเปลี่ยน relational=True เอาท์พุทจะแสดงเป็น (-oo < x) & (x < oo)

คำสั่ง solve_univariate_inequality ใช้ได้กับฟังก์ชัน polynomial ด้วย

14. แสดงสมการให้ดูเป็นปกติใน Spyder

ถ้าคุณรันโปรแกรมบน jupyter notebook ก็ข้ามหัวข้อนี้ไปได้ หากต้องการรันบน spyder จะดูยากแต่มีวิธีช่วยโดยเราสามารถแปลงเอาท์พุทเหล่านี้ให้อยู่ในรูปสมการปกติโดยเปลี่ยนให้เป็นรูปกราฟฟิคแบบ png ซึ่งสามารถ copy ไปแสดงในที่อื่นๆได้

from sympy import *
from IPython.display import Image, display
from IPython.lib.latextools import latex_to_png
def equToPng(equ, show=False):
     lt = latex(equ)
     data = latex_to_png(lt, wrap=True, color='white')
     if show:
         display(Image(data=data))
#example
x, y = symbols('x y')
eq = x**2 + 6*x*y + 8
equToPng(eq, True)
eq = sqrt(2)*sin(x)**2 + sqrt(x*y**2)/2
equToPng(eq, True)

ให้คุณ copy code ในบรรทัดที่ 1-8 แทรกเข้าไปข้างบนของทุกโปรแกรม จากนั้นเปลี่ยนคำสั่งจาก print(eq) ไปเป็น equToPng(eq, True) ก็จะได้รูปสมการปกติปรากฏใน console ซึ่งเป็นรูป png ที่สามารถ save ไว้ได้

อีกวิธีที่สะดวกกว่าคือทำเป็น package แล้ว import เข้ามา อ่านบทความนี้

15. สรุป

คงจะพอให้เห็นภาพรวมของการใช้ sympy กับงานทางคณิตศาสตร์ ง่ายๆคือเริ่มจากการเขียนสัญลักษณ์(symbolic หรือ sym ตามชื่อ package)ก็จะได้สมการจากนั้นก็ค่อย expand, factor, caculate และ solve ตามต้องการ ให้ลองทดลองเปลี่ยนสมการเล่นดู รายละเอียดที่ซับซ้อนพบได้ในตอนต่อไป

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