การดิฟและอินทีเกรทด้วย python/sympy

Visits:113

Sympy นอกจากจะเป็นแพคเกจเกี่ยวกับสัญลักษณ์ทางคณิตศาสตร์และแก้สมการพื้นๆได้แล้วยังสามารถ diff และ integrate ได้อีกด้วย

ทดลองรันโดยการ copy โปรแกรมไปใส่ใน jupyter notebook จะสะดวกที่สุดโดยคำสั่ง print( expr ) ให้เปลี่ยนเป็น expr แต่หากรันใน editor/tools อื่นๆก็ copy/paste ไปใส่ตรงๆได้เลยแต่ output ที่ออกมาจะดูลำบาก

คุณควรจะอ่านพื้นฐานของ sympy มาก่อนแล้ว ถ้ายังแนะนำให้อ่านก่อน

1. การดิฟ (differentiation) ด้วย sympy

Example 1.1

เริ่มจากฟังก์ชั่นพื้นฐาน $y=x^2$ ซึ่งเรารู้ว่า $\frac{d}{dx}x^2 = 2x$

from sympy import symbols
x = symbols('x')
expr = x**2
expr = expr.diff(x)
print( expr )

เราสามารถสร้าง symbol เป็นตัวแปรค้างไว้ก็ได้เช่น

\frac{d}{dx}(x^{-m}) = - mx^{-m-1}
from sympy import symbols
x, m = symbols('x m')
expr = x**m
expr = expr.diff(x)
print(expr)

แต่ตอน diff ในคำสั่งที่ 4 จะต้องระบุว่า diff เทียบกับอะไร สังเกตว่า sympy ไม่ได้ให้คำตอบในรูปง่ายๆ $-mx^{-m-1}$ ตามที่คาดหวัง แต่ให้คำตอบเป็น $\frac{m x^{- m}}{x}$ แม้ว่าจะ simplify แล้วก็ตาม(exercise)

Example 1.2

สมการต่อไปนี้เมื่อ diff ด้วยมือโดยใช้สูตรผลหาร $\frac{d}{dx}(\frac{u}{v}) = \frac{vdu/dx – udv/dx}{v^2}$ จะได้

\frac{d}{dx}(\frac{2x+5}{3x-2}) = -\frac{19}{\left(3 x - 2\right)^{2}}
from sympy import symbols
x = symbols('x')
expr = (2*x+5)/(3*x-2)
expr = expr.diff(x)
print(expr)

คำตอบคือ

- \frac{3 \left(2 x + 5\right)}{\left(3 x - 2\right)^{2}} + \frac{2}{3 x - 2}

ซึ่งไม่ตรงกับที่คาดหวัง(sympy มีอัลกอริธึมในการคิดเป็นของตัวเอง) แต่หากเราใช้ simplify(eq) หรือ together(eq) ก็จะได้ผลตรงกับที่คาดหวัง (exercise)

Example 1.3

ใช้ diff ผลหารเช่นเดียวกับตัวอย่างก่อนหน้า

\begin{aligned}
\dfrac{d}{dx}(\frac{\sin x}{1-\cos x}) &= \frac{(1-\cos x)(\cos x) - \sin^2 x}{(1- \cos x)^2} \\
&= \frac{\cos x - \cos^2 x -\sin^2 x}{(1- \cos x)^2} \\&= \frac{1}{\cos x - 1}
\end{aligned}
from sympy import symbols, sin, cos
x = symbols('x')
expr = sin(x)/(1-cos(x))
expr = expr.diff(x)
print(expr)

ผลออกมาได้ตรงตามคาดหวังเช่นกัน

\frac{\cos{\left(x \right)}}{1 - \cos{\left(x \right)}} - \frac{\sin^{2}{\left(x \right)}}{\left(1 - \cos{\left(x \right)}\right)^{2}}

ถ้าเราใช้ simplify(eq), trigsimp(eq) ก็จะได้ผลตามความหวังเช่นกัน แต่ together(eq) ลดรูปได้ไม่สุด ลองดู

Example 1.4

ตัวอย่างนี้เรา diff ด้วยมือด้วยสูตรผลคูณ $\frac{d}{dx}(uv) = udv/dx + vdu/dx$

\begin{aligned}
\frac{d}{dx}(\sec x\csc x) &= \sec(x)(-\csc(x)\cot(x)) + \csc(x)(\sec(x)\tan(x))  \\
&= -\frac{1}{\cos(x)}\frac{1}{\sin(x)}\frac{\cos(x)}{\sin(x)} + \frac{1}{\sin(x)}\frac{1}{\cos(x)}\frac{\sin(x)}{\cos(x)} \\
&= -\frac{1}{\sin^2(x)}  + \frac{1}{\cos^2(x)}
\end{aligned}
from sympy import symbols, sec, csc
x = symbols('x')
expr = sec(x)*csc(x)
expr = expr.diff(x)
print(expr)

คำตอบจากโปรแกรมคือ

\tan{\left(x \right)} \csc{\left(x \right)} \sec{\left(x \right)} - \cot{\left(x \right)} \csc{\left(x \right)} \sec{\left(x \right)}

ซึ่งหากใช้ simplify(eq) หรือ trigsimp(eq) ก็จะได้รูปตามที่คาดหวัง มีข้อสังเกตอันหนึ่งหากใช้ together(eq) จะได้

\left(\tan{\left(x \right)} - \cot{\left(x \right)}\right) \csc{\left(x \right)} \sec{\left(x \right)}

ซึ่งเป็นรูปฟอร์มที่ง่ายกว่าคำตอบจากโปรแกรมอันแรกและอาจเป็นรูปที่เราต้องการ

2. การอินทีเกรท (integration) ด้วย sympy

Example 2.1

\int x^2 \;dx = \frac{1}{3}x^3 + C
from sympy import symbols, integrate
x = symbols('x')
expr = x**2
expr = expr.integrate(x)
print(expr)

สังเกตว่าเอาท์พุทที่ได้ไม่มีค่าคงที่

Example 2.2

ถ้าต้องการหา definite integral

\int_0^2 x^2 \;dx = \frac{1}{3}x^3 \bigg |_0^3 = 9
from sympy import symbols
x = symbols('x')
expr = x**2
expr = expr.integrate((x, 0, 3))
print(expr)

Example 2.3

หากมีตัวแปรมากกว่าหนึ่งเช่น

\int (x^2 + y) dy = x^2y + \frac{y^2}{2}
from sympy import symbols, sin
x, y = symbols('x y')
expr = x**2 + y
expr = expr.integrate(y)
print(expr)

ให้ระบุว่าต้องการ integrate เทียบกับ y

Example 2.4

ทำด้วยมือด้วย integrate by part

\begin{aligned}
\int x e^{4x}dx &=\frac{1}{4}\int x \;d(e^{4x})=\frac{1}{4}xe^{4x} - \frac{1}{4}\int e^{4x}dx \\&= \frac{1}{4}xe^{4x}-\frac{1}{16}e^{4x}
\end{aligned}

Exercise คำตอบจากโปรแกรมคือ $\frac{\left(4 x – 1\right) e^{4 x}}{16}$

Exercise แต่หากใช้คำสั่ง cancel() ช่วยแยกเทอมจะได้รูปเดียวกับที่ทำด้วยมือข้างต้น

expr = expr.integrate(x).cancel()

Example 2.5

ค่า limit สามารถติดเป็นตัวแปรได้

\int_a^b e^x dx = e^b-e^a
from sympy import symbols, exp
x, a, b = symbols('x a b')
expr = exp(x)
expr = expr.integrate( (x, a, b) )
print(expr)

Example 2.6

ใช้อักษร “โอ” สองตัว oo แทน $\infty$ หรือ -oo แทน $-\infty$

\int_0^\infty e^{-x^2}dx = \frac{\sqrt{\pi}}{2}
from sympy import symbols, exp, oo
x, a, b = symbols('x a b')
expr = exp(-x**2)
expr = expr.integrate( (x, 0, oo) )
print(expr)

Example 2.7

การอินทีเกรทสองชั้น (double integral)

\int\int xy^2dxdy = \int \frac{1}{2}x^2 y^2 dy = \frac{1}{6}x^2 y^3
from sympy import symbols
x, y = symbols('x y')
expr = x*y**2
expr = expr.integrate( x, y )
print(expr)

แต่หากสลับลำดับการอินทีเกรท

\int\int xy^2dydx = \int \frac{1}{3}x y^3 dx = \frac{1}{6}x^2 y^3

ก็เปลี่ยนคำสั่งที่ 4 เป็น expr.integrate( y, x ) ซึ่งโดยทั่วไปจะเท่ากันยกเว้น definite integral เช่น

\int_0^1\int_0^2 xy^2dxdy = \int_0^1 [ \frac{1}{2}x^2 \Big |_0^2] y^2 dy = \int_0^1 2y^2 dy = \frac{2}{3}y^3 \Big |_0^1 = \frac{2}{3}
from sympy import symbols, exp
x, y = symbols('x y')
expr = x*y**2
expr = expr.integrate( (x, 0, 2), (y, 0 , 1) )
print(expr)

แต่หากสลับลำดับการอินทีเกรทเป็น dydx

\int_0^1 \int_0^2 xy^2dydx = \frac{4}{3}

คำสั่งที่ 4 ต้องเปลี่ยนเป็น expr.integrate( (y, 0, 2), (x, 0 , 1) )

Example 2.8

หา $\int x^3\sqrt{1-x^2}$ ด้วยมือด้วยวิธี substitution

\begin{aligned}
u &= 1-x^2, \;\;\; du = -2xdx, \\
dx &= -\frac{du}{2x}=-\frac{du}{2(1-u)^{1/2}}\\
x^2 &= 1-u, \;\;\; x^3 = (1-u)^{3/2} \\
-\frac{1}{2}\int \frac{(1-u)^{3/2}u^{1/2}}{(1-u)^{1/2}}du &= -\frac{1}{2}\int(1-u)u^{1/2}du \\
&= (\frac{1}{5}u^{5/2} - \frac{1}{3}u^{3/2}) + C \\
&= (\frac{1}{5}(1-x^2)^{5/2} - \frac{1}{3}(1-x^2)^{3/2}) + C
\end{aligned}
from sympy import sqrt, symbols
x = symbols('x')
expr = x**3*sqrt(1-x**2)
expr = expr.integrate(x)
print(expr)
\frac{x^{4} \sqrt{1 - x^{2}}}{5} - \frac{x^{2} \sqrt{1 - x^{2}}}{15} - \frac{2 \sqrt{1 - x^{2}}}{15}

เอาท์พุทจาก sympy ต่างจากที่ทำด้วยมือมากแต่หาก

Exercise ลองใช้ expr.integrate(x).simplify() เพื่อจัดรูปเอาท์พุทใหม่ ซึ่งจะได้

 \frac{\sqrt{1 - x^{2}} \left(3 x^{4} - x^{2} - 2\right)}{15}

Exercise ลอง simplify คำตอบที่ทำด้วยมือด้วย simplify(‘(1/5)*(1-x**2)**(5/2) – (1/3)*(1-x**2)**(3/2)’) ก็จะได้ผลออกมาตรงกับที่ simplify เอาท์พุทของ sympy (เทอมขวา)

หรือลองทดสอบโดยการ diff กลับก็จะได้ integrand ดั้งเดิม $x^3\sqrt{1-x^2}$

eq.integrate(x).diff(x).simplify()

ที่ยกตัวอย่างมาอยากให้เข้าใจว่ารูปฟอร์มของสมการที่ออกมาจาก sympy ไม่จำเป็นต้องตรงกับที่ทำด้วยมือเพราะ sympy ไม่ได้ใช้วิธีคิดแบบเดียวกับที่เราเรียนในตำราแคลคูลัสหรือมีอัลกอริธึมต่างกัน แต่ถ้าลอง simplify ดูก็จะได้ผลตรงกัน ส่วนเราจะใช้แบบฟอร์มใดขึ้นกับงานที่เราทำอยู่

Example 2.8

ลองใช้ sympy หาคำตอบทั่วไปของการอินทีเกรทว่าตรงกับในสูตรตามตำราหรือไม่

\begin{aligned}
\int (ax+b)^n dx &= \frac{\left(a x + b\right)^{n + 1}}{a(n + 1)}  + C \\
&= \frac{1}{a}\ln (a x + b)
\end{aligned}
from sympy import symbols, integrate
a, b, x, n = symbols('a b x n')
eq = (a*x + b)**n
eq.integrate(x)
\frac{\begin{cases} \frac{\left(a x + b\right)^{n + 1}}{n + 1} & \text{for}\: n \neq -1 \\\log{\left(a x + b \right)} & \text{otherwise} \end{cases}}{a}

3. การอินทีเกรท polynomial ด้วย numpy

สุดท้ายขอจบด้วยการอินทีเกรทฟังก์ชั่น polynomial ด้วยแพคเกจ numpy ซึ่งบางครั้งสะดวกกว่าการใช้ sympy เพราะสามารถสร้างสมการโดยคำสั่ง poly1d ซึ่งใส่เพียง coefficient ของสมการเข้าไปตรงๆซึ่งสะดวกในการเขียนมากกว่า เช่น

\int (3x^2 + 2x +7)dx = x^3 + x + 7x + C
from numpy import poly1d
p = poly1d([3, 2, 7])
print(p)
print( p(2) )

หมายถึงสร้างตัวแปร p ให้แทนสมการ $3x^2 + 2x + 7$, คำสั่งที่ 4 p(2) หมายถึงหาค่า $3(2)^2 + 2(2) + 7 = 23$

หากต้องการหาค่า definite integral ให้ใช้ polyint

\int_1^2 (3x^2 + 2x +7)dx = result[2] -result[1]
from numpy import poly1d, polyint
p = poly1d([3, 2, 7])
result = polyint(p)
print( result )
print( result(2) - result(1) )

คำสั่ง 4 จะเป็นผลจากการอินทีเกรทซึ่ง result ก็คือ $x^3 + x^2 + 7x$

คำสั่ง 5 หากต้องการหา definite integral ก็เพียงแต่เอาผลจากคำสั่ง 4 มาลบกัน

Facebook Comments