การหาผลเฉลยของสมการเชิงอนุพันธ์สามัญด้วยวิธีการแปลงลาปลาซและเขียนโค้ด Python#

บทความนี้กล่าวถึง ตัวอย่างการเขียนโค้ด Python เพื่อหาผลเฉลยของสมการเชิงอนุพันธ์สามัญด้วยวิธีการแปลงลาปลาซ

Keywords: Python, SymPy, Ordinary Differential Equations, ODE Solving, Symbolic Methods, Laplace Transforms, RC Circuit

---

▷ สมการเชิงอนุพันธ์สามัญและการหาผลเฉลย#

การวิเคราะห์วงจรไฟฟ้าพื้นฐานที่ประกอบด้วยตัวต้านทาน (Resistor: R) ตัวเก็บประจุ (Capacitor: C) ตัวเหนี่ยวนำไฟฟ้า (Inductor: L) ซึ่งเรามักจะเรียกวงจรประเภทนี้ว่า วงจร RLC

สำหรับวงจรประเภทนี้ ซึ่งถือว่าเป็นวงจรประเภทเชิงเส้น เราสามารถใช้สมการเชิงอนุพันธ์สามัญ (ODE) เพื่อสร้างแบบจำลองทางคณิตศาสตร์ สามารถอธิบายพฤติกรรมของวงจรไฟฟ้า หรือวิเคราะห์ดูการเปลี่ยนแปลงของปริมาณไฟฟ้าในเชิงเวลาได้ (Transient Analysis) เช่น แรงดันไฟฟ้าหรือกระแสไฟฟ้าในวงจร

ถ้าต้องการหาผลเฉลยของสมการเชิงอนุพันธ์สามัญที่มีคุณสมบัติต่อไปนี้ เราสามารถใช้วิธีที่เรียกว่า "การแปลงลาปลาซ" (Laplace Transform) ในการหาผลเฉลยได้

1. สมการเชิงอนุพันธ์สามัญนั้นเป็นแบบเชิงเส้นและมีสัมประสิทธิ์คงตัว (Linear ODE with Constant Coefficients)
2. มีการกำหนดเงื่อนไขเริ่มต้น (Initial Conditions) สำหรับฟังก์ชันที่ใช้เป็นผลเฉลยของสมการ ODE
3. ฟังก์ชันเจาะจงในสมการเชิงอนุพันธ์สามัญนั้นเป็นฟังก์ชันที่สามารถหาผลการแปลงลาปลาซได้ เช่น ฟังก์ชันที่เป็นค่าคงที่ ฟังก์ชันขั้นบันได ฟังก์ชันเลขชี้กำลัง-เอกซ์โพเนนเชียล ฟังก์ชันไซน์หรือและโคไซน์ เป็นต้น

ฟังก์ชันและอนุพันธ์ของฟังก์ชัน (Derivatives) ในโดเมนของตัวแปร $t$ จะถูกแปลงให้อยู่ในรูปของฟังก์ชันในโดเมน $s$ โดยทั่วไปแล้วจะอยู่ในรูปของเศษส่วนพหุนาม (Rational Polynomials) ของตัวแปร $s$ จากนั้นจึงสามารถใช้เทคนิคต่าง ๆ ทางพีชคณิต (Algebraic Manipulation) กับผลการแปลงลาปลาซดังกล่าว แล้วนำไปหาผลการแปลงลาปลาซผกผัน (Inverse Laplace Transform) เพื่อให้ได้ฟังก์ชันที่เป็นผลเฉลยในโดเมน $t$

ลองมาดูตัวอย่าง สมการเชิงอนุพันธ์สามัญอันดับที่ $n$ และเป็นไปตามข้อ 1

$$a_n y^{(n)}(t) + a_{n-1} y^{(n-1)}(t) + ... + a_1 y'(t) + a_0 y(0) = x(t) \\ $$

• $a_n, ..., a_0$ เป็นค่าคงที่
• $t$ เป็นตัวแปรอิสระ (Independent Variable) หมายถึง เวลา (Time)
• $y(t)$ เป็นฟังก์ชันที่ขึ้นอยู่กับตัวแปร $t$ และเป็นผลเฉลยของ ODE หมายถึง ปริมาณทางไฟฟ้าในวงจร
• $y'(t)$, ..., $y^{(n)}(t)$ เป็นอนุพันธ์อันดับที่ 1 ถึง $n$ ของฟังก์ชัน $y(t)$ ตามลำดับ
• $x(t)$ เป็นฟังก์ชันที่ขึ้นอยู่กับตัวแปร $t$ ตามรูปแบบของแหล่งจ่ายในวงจร และสามารถหาผลการแปลงลาปลาซได้

และปัญหานี้มีเงื่อนไขเริ่มต้นสำหรับ $t=0$ ตามข้อ 2 ดังนี้

$$y^{(n-1)}(0) = y_{n-1}, ..., y'(0) = y_1, y(0) = y_0$$

 

---

▷ การแปลงลาปลาซด้านเดียว#

การแปลงลาปลาซของฟังก์ชันที่มีตัวแปรอิสระ $t$ เป็นการดำเนินการทางคณิตศาสตร์ โดยอาศัยการอินทิเกรต และได้ผลลัพธ์จากการหาอินทิกรัลเป็นฟังก์ชันที่มีตัวแปรอิสระ $s$ (เป็นตัวแปรเชิงซ้อน) สำหรับช่วงเวลาตั้งแต่ $t=0^-$ ไปจนถึง $t=\infty$ ดังนี้

$$\mathcal{L}\{\, f(t) \,\} = F(s) \\[1em] F(s) = \int_{0^-}^\infty f(t)\, e^{-s t} dt, \; t \geq 0 \\ $$

หรือกล่าวได้ว่า ฟังก์ชัน $f(t)$ และ $F(s)$ เป็นคู่ของการแปลงลาปลาซ (Laplace Transform Pair)

การหาผลการแปลงลาปลาซในรูปแบบนี้ เป็นแบบด้านเดียว (One-sided Laplace Transform) ซึ่งหมายถึง ช่วงเวลาตั้งแต่ $t \geq 0$ เป็นต้นไป

ถ้าให้ $f(t)=1$ ซึ่งเป็นฟังก์ชันที่ให้ค่าคงที่เท่ากับ 1 สำหรับค่าใด ๆ ของตัวแปร $t$ ผลการแปลงลาปลาซด้านเดียวของฟังก์ชันนี้ สามารถคำนวณได้ดังนี้

$$\begin{align} F(s) &= \int_{0^-}^\infty f(t)\, e^{-s t} dt \\ &= \int_{0^-}^\infty 1\cdot e^{-s t} dt \\ &= \Big[ \frac{1}{-s} e^{-s t} \Big]_{t=0^-}^{t=\infty}\\ &= \lim_{t\rightarrow \infty} \frac{1}{-s} e^{-s t} - \lim_{t\rightarrow 0^-} \frac{1}{-s} e^{-s t} \\ &= \frac{1}{s}, \;\; \text{ROC: $Re\{s\} > 0$ } \\ \end{align}$$

จากตัวอย่างการคำนวณ ถ้ากำหนดขอบเขตของการลู่เข้า (Region of Convergence: ROC) สำหรับตัวแปร s เช่น ถ้าให้ส่วนที่เป็นจำนวนจริงของตัวแปร $s$ หรือ $Re\{s\}$ มีค่ามากกว่า 0 ก็จะทำให้สามารถหาค่าลิมิตได้ในกรณีที่ $t \rightarrow \infty$

 

ฟังก์ชัน $u(t)$ ที่มีลักษณะเป็นขั้นบันไดหนึ่งหน่วย $u(t)$ (Unit-step Function) หรือเรียกว่า ฟังก์ชันเฮฟวีไซด์ $\theta(t)$ (Heaviside Function) มีค่าแบ่งเป็นช่วงดังนี้

$$u(t) = \begin{cases} 1, & \; t \geq 0 \\ 0, & \; t < 0 \\ \end{cases}\\ \theta(t) = \begin{cases} 1, & \; t > 0 \\ \frac{1}{2}, & \; t=0 \\ 0, & \; t < 0 \\ \end{cases}$$

ผลการแปลงลาปลาซด้านเดียวของฟังก์ชัน $f(t) = u(t)$ ให้ผลเหมือนกับกรณี $f(t) = 1$

$$\begin{align} F(s) &= \int_{0^-}^\infty f(t)\, e^{-s t} dt \\ &= \int_{0^-}^\infty u(t) \cdot e^{-s t} dt \\ &= \int_{0^-}^\infty 1 \cdot e^{-s t} dt \\ &= \frac{1}{s}, \;\; \text{ROC: $Re\{s\} > 0$} \\ \end{align}$$

 

ถ้าให้ $f(t) = e^{-a t}$ เป็นฟังก์ชันเอกซ์โพเนนเชียล (Exponential Function) และให้ $a$ เป็นค่าคงที่ ผลการแปลงลาปลาซด้านเดียวสามารถคำนวณได้ดังนี้

$$\begin{align} F(s) &= \int_{0^-}^\infty f(t)\, e^{-s t} dt \\ &= \int_{0^-}^\infty e^{-a t} \cdot e^{-s t} dt \\ &= \int_{0^-}^\infty e^{-(s+a)t} dt \\ &= \frac{1}{(s+a)}, \;\; \text{ROC: $Re\{s\} > -a$} \\ \end{align}$$

 

ลองมาดูตัวอย่างของคู่การแปลงลาปลาซสำหรับฟังก์ชันพื้นฐานทางคณิตศาสตร์ที่พบเห็นได้บ่อย โดยกำหนดให้ $c$ และ $a$ เป็นค่าคงที่ใด ๆ $\omega_0$ และ $t_0$ เป็นค่าคงที่และเป็นจำนวนจริงที่เป็นบวก

f(t)	F(s)	Description
$\delta(t)$	$1$	Dirac (unit-impulse) function
$c$	$\frac{c}{s}$	constant function
$u(t)$	$\frac{1}{s}$	Heaviside (unit-step) function
$u(t)-u(t-t_0)$	$\frac{1-e^{-s t_0}}{s}$	Pulse function
$e^{-at}$	$\frac{1}{s+a}$	Exponential function
$sin(\omega_0 t)$	$\frac{\omega_0}{s^2 + \omega_0^2}$	Sine function
$cos(\omega_0 t)$	$\frac{s}{s^2 + \omega_0^2}$	Cosine function
$e^{-at} sin(\omega_0t)$	$\frac{\omega_0}{(s+a)^2 + \omega_0^2}$	Exponentially damped sine
$e^{-at} cos(\omega_0t)$	$\frac{(s+a)}{(s+a)^2 + \omega_0^2}$	Exponentially damped cosine
$t\cdot u(t)$	$\frac{1}{s^2}$	Ramp function
$a\cdot f(t) \pm b\cdot g(t)$	$a\cdot F(s) \pm b\cdot G(s)$	Linearity property
$e^{\mp a t} \cdot f(t)$	$F(s\pm a)$	Shifting in s domain
$f'(t)$	$s F(s) - f(0)$	The first derivative of $f(t)$
$f''(t)$	$s^2 F(s) - s f(0) - f'(0)$	The second derivative of $f(t)$
$f^{(n)}(t)$	$s^n F(s) - \sum_{i=0}^{n-1} s^{n-i-1}\, f^{(i)}(0)$	The $n$-th derivative of $f(t)$

โดยทั่วไปแล้ว ตารางของคู่การแปลงลาปลาซ และคุณสมบัติต่าง ๆ ของการแปลงลาปลาซ มีประโยชน์และช่วยให้ง่ายขึ้นเมื่อต้องการหาผลการแปลงลาปลาซ และการแปลงลาปลาซผกผัน

 

---

▷ การเขียนโค้ด Python: การอนุพันธ์และอินทิกรัลสำหรับฟังก์ชัน#

ลองมาดูการเขียนโค้ด Python สำหรับ SymPy (v1.10.0+) ตรวจสอบเวอร์ชันของ SymPy โดยทำคำสั่งต่อไปนี้

import sympy
print(sympy.__version__)

 

SymPy มีคำสั่ง sympy.diff() และ sympy.integrate() เพื่อคำนวณหาอนุพันธ์ (Derivative) และการหาปริพันธ์หรืออินทิกรัลแบบไม่จำกัดเขต (Indefinite Integral)

ลองพิจารณาฟังก์ชัน เช่น $f(t) = e^{-a t}$ และ $f(t) = sin(\omega t)$ แล้วหาอนุพันธ์ และปริพันธ์ของฟังก์ชันเหล่านี้

$$\int e^{-a t} dt = -\frac{e^{-a t}}{a} + C \\ \frac{d}{dt} e^{-a t}= -a e^{-a t} \\ \int sin(\omega t) dt = -\frac{cos(\omega t)}{\omega} + C \\ \frac{d}{dt} sin(\omega t) = \omega cos(\omega t) \\ \frac{d^2}{dt^2} sin(\omega t) = -\omega^2 sin(\omega t) \\ $$

from sympy import *
from IPython.display import display
# define 't' as a symbol 
t = Symbol('t')
# define a symbol for function f(t) 
f = symbols('f', cls = Function)
# define 'a' (non-zero) as a symbol
a = Symbol( 'a', nonzero=True )
# f is an exponential function
f = exp(-a*t)
# compute the integral of f(t)
display( integrate( f, t ) )
# compute the first derivative of f(t)
display( diff( f, t) )
# define 'w' (positive) as a symbol 
w = Symbol( 'omega', positive=True )
# f is a sine function
f = sin(w*t)
# compute the integral of f(t)
display( integrate( f, t ) )
# compute the first and second derivatives of f(t)
for i in range(1,3):
    display( diff(f, t, i) )

 

หากต้องการทดสอบการทำงานของโค้ดตัวอย่าง และแสดงผลลัพธ์ที่ได้ที่เป็นข้อความทางคณิตศาสตร์ ก็สามารถใช้ JupyterLab ผ่านหน้าเว็บเบราว์เซอร์ได้ (ติดตั้งใช้งานด้วยคำสั่ง เช่น pip install jupyterlab)

รูป: การใช้งาน JupyterLab

รูป: เอาต์พุตสำหรับโค้ดตัวอย่างเมื่อรันโค้ดด้วย JupyterLab

 

ฟังก์ชันขั้นบันไดหนึ่งหน่วย (Unit-Step Function) หรือ "ฟังก์ชันเฮฟวิไซด์" (Heaviside Function) ซึ่งใช้สัญลักษณ์ $\theta(t)$ และ "ฟังก์ชันดิแรกเดลตา" (Dirac-Delta Function) ซึ่งใช้สัญลักษณ์ $\delta(t)$ ถือว่ามีความสำคัญในการวิเคราะห์ระบบและวงจรไฟฟ้า ลองมาทบทวนคุณสมบัติของฟังก์ชันดังนี้

$$\theta(t) = \begin{cases} 1 & t > 0 \\ \frac{1}{2} & t = 0 \\ 0 & t < 0 \\ \end{cases}\\[1em] \delta(t) = \begin{cases} \infty & t \neq 0 \\ 0 & t = 0 \\ \end{cases} \\[1em] \int_{-\infty}^{\infty} \delta(t) dt = 1 \\ \int_{-\infty}^{\infty} f(t) \delta(t) dt = \int_{-\infty}^{\infty} f(0) \delta(t) dt = f(0) \\ \int_{-\infty}^{t} \delta(\tau)d\tau = \theta(t) \\ \delta(t) = \frac{d}{dt} \theta(t) \\ $$

ในโค้ดตัวอย่างจะเห็นได้ว่า คำสั่งสำหรับการใช้ฟังก์ชัน $\delta(t)$ คือ DiracDelta() และคำสั่ง Heaviside() สำหรับฟังก์ชันเฮฟวิไซด์ $\theta(t)$ นอกจากนั้นยังมีตัวอย่างการใช้คำสั่ง Piecewise() เพื่อสร้างฟังก์ชันที่ให้ค่าแบ่งเป็นช่วง

การหาอนุพันธ์ของฟังก์ชัน $\theta(t)$ จะได้ฟังก์ชัน $\delta(t)$ เป็นผลลัพธ์ แต่ถ้าอินทิเกรต $\delta(t)$ สำหรับ $t\in (-\infty,\infty)$ จะได้เท่ากับ $1$ (ในการเขียนโค้ดด้วย SymPy ถ้าใช้สัญลักษณ์ sym.oo จะหมายถึง ค่าอนันต์ $\infty$)

from sympy import *
from IPython.display import display
# define a symbol for variable 't'
t = Symbol('t')
# display the Heaviside function
display( Heaviside(t) )
# compute the value of the Heaviside function at t=0
print( Heaviside(0) )
# integrate the Heaviside function
display( integrate( Heaviside(t) ) )
# find the derivative of the Heaviside function
display( diff(Heaviside(t),t ) )
# rewrite Heaviside function as a piecewise function
display( Heaviside(t).rewrite(Piecewise) )
a = Symbol( 'a', nonzero=True )
display( limit( integrate( DiracDelta(t),(t,-a,a)), a,oo) )
# define 'a' as a nonzero symbol 
a = Symbol( 'a', nonzero=True )
# integrate the dirac-delta function for t=(-inf,+inf)
display( limit( integrate( DiracDelta(t),(t,-a,a)), a,oo) )
# define f as a piecewise function 
f = Piecewise( (1,t>0), (0,t<0), (1/2,True) )
display( f )

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง

 

---

▷ การเขียนโค้ด Python: การหาผลการแปลงลาปลาซ#

ถัดไปเป็นตัวอย่างการคำนวณหาผลการแปลงลาปลาซของฟังก์ชัน $\delta(t)$, $\theta(t)$ และ $\theta(t-t_0), t_0 > 0$ ตามลำดับ โดยการใช้คำสั่ง integrate() เพื่ออินทิเกรตฟังก์ชันในช่วง $t\in [0^-,\infty)$ และใช้คำสั่ง laplace_transform() เพื่อการแปลงลาปลาซโดยเฉพาะ

from sympy import *
from IPython.display import display
# define symbols for variables: t and s
t = Symbol('t', real=True)
s = Symbol('s')
a = Symbol( 'a', real=True, negative=True )
# compute the Laplace transform using definite integral
F = integrate( DiracDelta(t)*exp(-s*t), (t,a,oo) )
display( simplify( limit( F, a, 0,'-') ) ) 
# compute the Laplace transform of \delta(t) using built-in function
# return the Laplace transform only
result = laplace_transform( DiracDelta(t), t, s, noconds=True )
display( result )
# return the Laplace transform and the convergence condition
lt,cond,_ = laplace_transform( DiracDelta(t), t, s )
display( lt )
print( f'ROC: Re(s) > {cond}' ) 
# find the Laplace transform of the Heaviside function
lt = laplace_transform( Heaviside(t), t, s, noconds=True )
display( lt )
# find the Laplace transform of the time-shifted Heaviside function
t0 = Symbol( 't_0', real=True, positive=True )
F = integrate( Heaviside(t-t0)*exp(-s*t), (t,a,oo) )
display( simplify( limit( F, a, 0,'-') ) ) 
lt,cond,_ = laplace_transform( Heaviside(t-t0), t, s )
display( lt )
print( f'ROC: Re(s) > {cond}' ) 

จากผลการรันโค้ดจะเห็นได้ว่า ผลการแปลงลาปลาซของฟังก์ชัน $\delta(t)$ จะได้เท่ากับ $1$

$$\begin{align} \mathcal{L}\{ \delta(t)\} &= \int_{0^-}^\infty \delta(t) e^{-s t} dt \\ &= \int_{0^-}^\infty \delta(t) \cdot e^{-s t}|_{t=0}\, dt \\ &= \int_{0^-}^\infty \delta(t) \cdot 1\, dt = 1 \\ \end{align}$$ และมีขอบเขตการลู่เข้า (ROC) สำหรับตัวแปรเชิงซ้อน $Re(s) > -\infty$ (ในกรณีนี้หมายถึง ทั้งหมดของระนาบเชิงซ้อน s-plane)

ผลการแปลงลาปลาซของฟังก์ชัน $\theta(t)$ จะได้เท่ากับ $1/s$

$$\begin{align} \mathcal{L}\{ \theta(t)\} &= \int_{0^-}^\infty \theta(t) e^{-s t} dt \\ &= \int_{0^-}^\infty 1 \cdot e^{-s t} dt \\ &= \frac{1}{s} \qquad (\text{ROC: } Re\{s\} > 0 )\\ \end{align}$$ และมีขอบเขตการลู่เข้า (ROC) สำหรับตัวแปรเชิงซ้อน $Re(s) > 0$ (ในกรณีนี้หมายถึง ครึ่งหนึ่งทางด้านซ้ายของระนาบเชิงซ้อน Half s-plane)

คำสั่ง laplace_transform() เมื่อเรียกใช้คำสั่งนี้ และกำหนดให้อาร์กิวเมนต์ noconds เป็น False (หรือไม่มีการระบุไว้) จะให้ผลลัพธ์เป็นฟังก์ชันที่ได้จากการแปลงลาปลาซ กและเงื่อนไขการลู่เข้าสำหรับการอินทิเกรต แต่ถ้าให้อาร์กิวเมนต์นี้มีค่าเป็น True จะได้เฉพาะผลการแปลงลาปลาซเท่านั้น

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง (ต่อจากส่วนแรก)

ข้อสังเกต: $|arg(s)| < \frac{\pi}{2}$ ซึ่งเป็นเงื่อนไขการลู่เข้าของการอินทิเกรต สำหรับผลการแปลงลาปลาซ หมายความว่าอย่างไร ?

คำสั่ง $arg(s)$ สำหรับ SymPy หมายถึง การหามุมของตัวแปร $s$ ที่เป็นเลขเชิงซ้อน

$$s = Re(s) + j\, Im(s) \\ arg(s) = tan^{-1}\Big(\frac{Im(s)}{Re(s)}\Big) = atan2(Im(s),Re(s))$$

ดังนั้น $|arg(s)| < \frac{\pi}{2}$ ก็หมายความว่า $Re(s)$ จะต้องเป็นบวก หรือ $Re(s) > 0$

 

ตัวอย่างการเขียนโค้ด Python เพื่อแสดงคู่ของฟังก์ชันและผลการแปลงลาปลาซที่พบเห็นได้บ่อยมีดังนี้

from sympy import *
from IPython.display import display, Math, Latex
# define symbols
t = Symbol('t', real=True)
s = Symbol('s')
t0, omega, a = symbols('t_0 omega_0 a', real=True)
# define a list of basic functions
functions = [ Heaviside(t),
    Heaviside(t) - Heaviside(t-t0),
    exp(-a*t), t*exp(-a*t),  t**2*exp(-a*t),
    sin(omega*t), cos(omega*t),
    exp(-a*t)*sin(omega*t), exp(-a*t)*cos(omega*t) ]
# find Laplace transforms of a function list
LT = lambda f: laplace_transform(f,t,s,noconds=True)
results = [ LT(f) for f in functions]
# display Laplace transform pairs in Latex style
for f,F in zip(functions,results):
    ltx = f'$$\\mathcal{{L}} \\{{ {latex(f)} \\}}'
    ltx += '\\quad\\rightarrow\quad'
    ltx += latex(F) + '$$'
    display( Math(ltx) )

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง

 

---

▷ การแยกฟังก์ชันเศษส่วนออกเป็นเศษส่วนย่อย#

เราจะพบว่า ผลการแปลงลาปลาซมักอยู่ในรูปของเศษส่วนย่อย (Rational Polynomials) เช่น $F(s) = N(s)/D(s)$ โดยที่ $N(s)$ และ $D(s)$ เป็นฟังก์ชันพหุนามที่มีตัวแปร $s$

ในกรณีนี้ เราสามารถเขียนรูปใหม่ของฟังก์ชัน $F(s)$ ให้เป็นผลรวมของฟังก์ชันเศษส่วนย่อย หรือ การแยกฟังก์ชันเศษส่วนออกเป็นเศษส่วนย่อย (Partial Fraction Expansion)

ตัวอย่างเช่น

$$\begin{align} F(s) &= \frac{(5s - 5)}{ (s^2+4)(s+1) } \\ &= \frac{2s + 3}{s^2+4} - \frac{2}{s+1} \\ &= 2\Big(\frac{s}{s^2+2^2}\Big) + \frac{3}{2}\Big(\frac{2}{s^2+2^2}\Big) - 2\Big(\frac{1}{s+1}\Big), \;\; (t \geq 0) \\ \end{align}$$

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

$$f(t) = 2 cos(2 t) + \frac{3}{2} sin(2 t)- 2 e^{-t}$$

 

Python ก็มีคำสั่ง apart() สำหรับการแยกเป็นเศษส่วนย่อย และนำมาใช้กับฟังก์ชัน เศษส่วนที่มีตัวแปร $s$ แล้วนำไปหาผลการแปลงลาปลาซผกผันของแต่ละพจน์ย่อยได้ ตามตัวอย่างต่อไปนี้

from sympy import *
from IPython.display import display, Math
# define symbols
t = Symbol('t', real=True)
s = Symbol('s')
# define a rational polynomial function F(s) 
F = 5*(s-1)/((s**2 + 4)*(s+1))
# perform partial fraction expansion of F(s)
Fs = apart(F)
# show the result of partial fraction expansion
display( Math( '$$F(s) = $$' + latex(Fs) ) )
# split the expression into a list of terms
terms = Fs.args
# perform inverse Laplace transform for each term in the list
ILT = lambda term: inverse_laplace_transform( term, s, t )
f = 0
for term in terms:
    f += ILT(term)
# show the inverse Laplace tranform of F(s) --> f(t)
display( Math('$$f(t) = $$' + latex(f)) )

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง

 

---

▷ การหาผลเฉลยของสมการอนุพันธ์สามัญด้วยวิธีการแปลงลาปลาซ#

ถัดไปเป็นตัวอย่างของสมการ ODE อันดับสอง มีสัมประสิทธิ์คงตัวและเป็นสมการเชิงเส้น

$$y'' + 2y' + y = sin(2t) \\ y(0)=−1, y'(0)=1 \\ $$

การหาผลเฉลยด้วยการแปลงลาปลาซมีขั้นตอนดังนี้ เริ่มต้นให้หาผลการแปลงลาปลาซสำหรับฟังก์ชัน $y(t)$ และอนุพันธ์ของฟังก์ชันที่มีอยู่ในสมการในเชิงสัญลักษณ์ รวมถึงฟังก์ชัน $sin(2t)$ ที่อยู่ทางขวามือของสมการ

$$\begin{align} \mathcal{L}\{ y(t) \} &= Y(s) \\ \mathcal{L}\{ y'(t) \} &= s Y(s) - y(0) \\ \mathcal{L}\{ y''(t) \} &= s^2 Y(s) - s y(0) - y'(0) \\ \mathcal{L}\{ sin(2 t) \} &= \frac{2}{s^2 + 4}\\ \end{align}$$

แล้วนำไปแทนลงในสมการ ODE จะได้เป็น

$$\begin{align} [s^2 Y(s) - s y(0) - y'(0)] + 2[s Y(s) - y(0)] + Y(s) &= \frac{2}{s^2 + 4} \\ Y(s)(s^2 + 2s + 1) - s y(0) - 2 y(0) - y'(0) &= \frac{2}{s^2 + 4} \\ \end{align}$$

หากพิจารณาเฉพาะปัญหาที่เป็นเอกพันธ์ (Homogeneous Problem) ด้านขวามือของสมการ ODE จะเป็น 0 ดังนี้

$$Y_h(s)(s^2 + 2s + 1) - s y(0) - 2 y(0) - y'(0) = 0 \\ Y_h(s) = \frac{ s y(0) + 2 y(0) + y'(0) }{ (s+1)^2 } \\ $$ โดยที่ $Y_h(s)$ เป็นผลการแปลงลาปลาซของผลเฉลย $y_h(t)$ สำหรับสมการที่เป็นเอกพันธ์

ถ้าลองจัดรูปใหม่สำหรับ $Y_h(s)$ จะได้ผลรวมของพจน์ย่อยต่อไปนี้ $$\begin{align} Y_h(s) &= \frac{ s y(0) + 2 y(0) + y'(0) }{ (s+1)^2 } \\ &= \frac{ (s+1)y(0) + y(0) + y'(0) }{ (s+1)^2 } \\ &= \frac{ y(0)}{(s+1)} + \frac{y(0) + y'(0) }{ (s+1)^2 } \\ \end{align}$$

และนำไปหาผลการแปลงลาปลาซผกผันได้

$$\begin{align} y_h(t) &= \underbrace{y(0)}_{C_1}\, e^{-t} + \underbrace{(y(0) + y'(0))}_{C_2}\, t\, e^{-t} \\ &= C_1\, e^{-t} + {C_2}\, t\, e^{-t} \\ \end{align}$$ โดยที่ $C_1$ และ $C_2$ เป็นค่าคงที่ใด ๆ ซึ่งขึ้นอยู่กับเงื่อนไขเริ่มต้นของ ODE

 

ในกรณีที่สมการไม่เป็นเอกพันธ์ (Nonhomogeneous Problems) แต่ให้เงื่อนไขเริ่มต้นเป็น 0 ซึ่งกำหนดให้ $y(0) = 0, y'(0) = 0$ หลังจากจัดรูปใหม่ของสมการ จะได้ $Y_p(s)$ ดังนี้

$$Y_p(s) = \frac{2}{ (s+1)^2 (s^2 + 4) } \\ $$ โดยที่ $Y_p(s)$ เป็นผลการแปลงลาปลาซของผลเฉลยเจาะจง $y_p(t)$

หลังจากได้แยกฟังก์ชันเศษส่วน $Y_p(s)$ ออกเป็นเศษส่วนย่อย จะได้เป็น

$$\begin{align} Y_p(s) &= \frac{ 2/5 }{(s+1)^2} + \frac{ 4/25 }{(s+1)} - \frac{ (4s + 6)/25 }{(s^2 + 4)} \\ \end{align}$$

 

ผลเฉลยทั่วไป (General Solution) ของปัญหานี้คือ $Y(s)$ ซึ่งได้จากผลเฉลยจากทั้งสองกรณีนำมารวมกัน (เมื่อแปลงลาปลาซแล้ว)

$$\begin{align} Y(s) &= Y_h(s) + Y_p(s) \\ &= \frac{ y(0)}{(s+1)} + \frac{y(0) + y'(0) }{ (s+1)^2 } \\ & \;\; + \frac{ 2/5 }{(s+1)^2} + \frac{ 4/25 }{(s+1)} - \frac{ (4s + 6)/25 }{(s^2 + 4)} \\ \end{align}$$

จากเงื่อนไขเริ่มต้น $y(0)=−1, y'(0)=1$ แล้วแทนค่าลงในสมการสำหรับ $Y(s)$ จะได้

$$\begin{align} Y(s) &= Y_h(s) + Y_p(s) \\ &= \frac{(-1)}{(s+1)} + \frac{(-1+1)}{ (s+1)^2 } \\ & \;\; + \frac{ 2/5 }{(s+1)^2} + \frac{ 4/25 }{(s+1)} - \frac{ (4s + 6)/25 }{(s^2 + 4)} \\ &= \frac{ 2/5 }{ (s+1)^2 } - \frac{ 21/25 }{ (s+1) } - \frac{ (4s + 6)/25 }{(s^2 + 4)} \\ \end{align}$$

การหาผลการแปลงลาปลาซผกผันของ $Y(s)$ จะได้ผลลัพธ์เป็นฟังก์ชัน $y(t)$ ดังนี้ $$y(t) = \frac{2}{5} t e^{-t} - \frac{21}{25} e^{-t} - \frac{4}{25} cos(2t) - \frac{3}{25}sin(2t)$$

ถ้าให้ $cos(B) = 4/5$ และ $sin(B) = 3/5$ จะได้ $$sin^2(B) + cos^2(B) = 1 \\ tan(B) = \frac{sin(B)}{cos(B)} = \frac{3}{4} \\ B = tan^{-1}( \frac{3}{4} ) \\ $$ จากสูตร $$cos(A-B) = cos(A)cos(B) + sin(A)sin(B) \\\ $$ ดังนั้น $y(t)$ สามารถเขียนให้อยู่ในรูปใหม่ได้ดังนี้ $$\begin{align} y(t) &= \frac{2}{5} t e^{-t} - \frac{21}{25} e^{-t} - \frac{1}{5}\Big[ cos(2t) (\frac{4}{5}) + sin(2t)(\frac{3}{5}) \Big] \\ &= \frac{2}{5} t e^{-t} - \frac{21}{25} e^{-t} - \frac{1}{5} cos( 2t - atan(\frac{3}{4}) ) \end{align}$$

 

ตัวอย่างการเขียนโค้ด SymPy เพื่อหาผลเฉลยของสมการเชิงอนุพันธ์ โดยใช้คำสั่ง dsolve() มีดังนี้

from sympy import *
from IPython.display import display, Math
# define symbols
t = Symbol('t', real=True)
y = symbols('y', cls=Function)
# define D1y and D2y as the first and second derivative of y(t)
D1y = y(t).diff(t)
D2y = y(t).diff(t,t)
# write the ODE
ode = Eq( D2y + 2*D1y + y(t), sin(2*t) )
display( ode )
# define the initial conditions for the ODE
ics = {y(0): -1, D1y.subs(t, 0): 1} 
# solve the ODE without initial conditions
sol = dsolve( ode, y(t) )
# solve the ODE with initial conditions (IVP)
sol_ivp = dsolve( ode, y(t), ics=ics )
# display solutions 
display( sol )
display( sol_ivp )

เอาต์พุตที่ได้จากการรันโค้ดตัวอย่าง มีดังนี้

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง

 

หากต้องการเขียนโค้ดโดยใช้วิธีการแปลงลาปลาซ เพื่อหาผลเฉลยของสมการเชิงอนุพันธ์สามัญ (ODE) ก็มีแนวทางดังนี้

• สร้างฟังก์ชันที่เป็นสัญลักษณ์ y(t) และหาอนุพันธ์ของฟังก์ชันดังกล่าว แล้วนำไปแปลงลาปลาซ โดยใช้คำสั่ง laplace_transform() และนำผลลัพธ์ที่ได้ไปใช้เป็นพจน์ที่อยู่ทางด้านซ้ายของ ODE
• ฟังก์ชันที่อยู่ทางขวามือของสมการ และในตัวอย่างนี้คือ $sin(2t)$ ให้นำไปแปลงลาปลาซ แล้วนำไปใช้เป็นพจน์ที่อยู่ทางขวามือของ ODE
• ใช้สัญลักษณ์ Y(s) แทนที่ผลการแปลงลาปลาซของฟังก์ชัน y(t)
• สร้างเงื่อนไขเริ่มต้น (Initial Condition) แล้วนำไปแทนที่ลงในสมการ ODE หลังจากที่ได้มีการแปลงลาปลาซแล้ว
• จัดรูปสมการที่ทำให้ Y(s) ย้ายมาอยู่ด้านซ้ายมือของสมการ
• หากต้องการทำขั้นตอน Partial Fraction Expansion ให้ใช้คำสั่ง apart() ดำเนินการกับ Y(s)
• หาผลเฉลยโดยการแปลงลาปลาซผกผันของ Y(s)

from sympy import *
from IPython.display import display
# define symbols for independent variables
t = Symbol('t', real=True)
s = Symbol('s', complex=True)
# define symbols for functions
y = Function('y')(t)
Y = Function('Y')(s)
# define lambda function for Laplace transform
LT = lambda y: laplace_transform(y,t,s)
# find the Laplace Transform of y(t) and its derivatives
lt_y   = LT( y )
lt_D1y = LT( y.diff(t,1) )
lt_D2y = LT( y.diff(t,2) )
# write thhe Laplace tranform of the ODE
lt_ode = Eq( lt_D2y + 2*lt_D1y + lt_y, LT( sin(2*t) )[0] )
# substitute the Laplace transform of y(t) with Y(s)
Ys = LaplaceTransform(y,t,s)
lt_ode = lt_ode.subs( Ys, Y )
display( lt_ode )
# substitute the initital conditions into the ODE
ics = {y.subs(t,0): -1, y.diff(t).subs(t,0): 1}
lt_ode = lt_ode.subs( ics )
display( lt_ode ) 
# solve the ODE using Laplace transform
Y = solve( lt_ode, Y )[0]
# display Y(s)
display( Y )
# perform partial fraction expansion of Y(s)
Y = Y.apart( s ) 
display( Y )
# perform inverse Laplace transform
sol = inverse_laplace_transform( Y,s,t )
# simplify the solution y(t)
sol = simplify( expand( sol ) )
display( collect(expand( sol ), Heaviside(t)) )

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง (ต่อจากส่วนแรก)

 

---

▷ การวิเคราะห์วงจร RC ด้วยการแปลงลาปลาซ#

วงจรตัวอย่างต่อไปนี้ประกอบด้วยตัวต้านทาน $R_1$ และตัวเก็บประจุ $C_1$ ที่นำมาต่ออนุกรมกัน และมีแหล่งจ่าย $v_S(t)$ เป็นสัญญาณอินพุตของระบบและให้ $v_C(t)$ ซึ่งเป็นแรงดันตกคร่อมที่ตัวเก็บประจุ เป็นสัญญาณเอาต์พุตของระบบ รูปแบบของวงจร RC ในลักษณะนี้ ทำหน้าที่เป็นตัวกรองแบบพาสซีฟที่ยอมให้ความถี่ต่ำผ่านได้ดี (Low-Pass Passive RC Filter)

รูป: วงจร RC

ถ้าวิเคราะห์วงจรด้วยทฤษฎีทางไฟฟ้า KVL (Mesh Analysis) จะได้สมการต่อไปนี้

$$v_R(t) + v_C(t) = v_S(t) \\ \begin{align} i(t) &= i_C(t) = i_R(t) \\ v_R(t) &= R \, i_R(t) \\ i_C(t) &= C \frac{d}{dt} v_C(t) \\ \end{align}$$

และนำไปสู่สมการ ODE (อันดับหนึ่ง) $$RC \frac{d}{dt} v_C(t) + v_C(t) = v_S(t) \\ $$

ถ้าหาผลการแปลงลาปลาซ จะเขียนรูปสมการได้ดังนี้

$$\begin{align} \mathcal{L}\{ v_C(t) \} &= V_C(s) \\ \mathcal{L}\{ v'_C(t) \} &= s V_C(s) - v_C(0) \\ \mathcal{L}\{ v_S(t) \} &= V_S(s) \\ \end{align} \\ RC\, [s\,V_C(s) - v_c(0)] + V_C(s) = V_S(s) \\ (1 + sRC) V_C(s) = RC v_c(0) + V_S(s) \\ $$

ถ้าให้เงื่อนไขเริ่มต้นของวงจรเป็น $0$ (Zero Initial Condition) ซึ่งหมายถึง $v_C(0) = 0$ ดังนั้นจะได้ $V_C(s)$ $$\begin{align} V_C(s) &= \frac{1}{(1+sRC)} V_S(s) \\ H(s) &= \frac{V_{out}}{V_{in}} = \frac{V_C(s)}{V_S(s)} \\ &= \frac{1}{1+sRC} = \frac{ \frac{1}{RC} }{s + \frac{1}{RC}} \\ \end{align}$$

$H(s)$ เรียกว่า ฟังก์ชันถ่ายโอน (Transfer Function) ของวงจร และเป็นอัตราส่วนระหว่างสัญญาณเอาต์พุตต่อสัญญาณอินพุต และ $h(t)$ คือ ผลตอบสนองอิมพัลส์ของระบบ (Impulse Response) ซึ่งได้จากการแปลงลาปลาซผกผันของฟังก์ชัน $H(s)$

$$\begin{align} h(t) &= \mathcal{L}^{-1}\{ H(s) \} \\ &= \frac{1}{RC} e^{-t/RC} \\ &= \frac{1}{\tau} e^{-t/\tau}, \;\; t \geq 0, \; \tau = RC \\ \end{align}$$

ถัดไปมาดูตัวอย่างการเขียนโค้ด Python สำหรับการวิเคราะห์วงจรด้วยการแปลงลาปลาซ หารูปแบบของฟังก์ชันถ่ายโอน $H(s)$ และฟังก์ชันที่เป็นผลตอบสนองอิมพัลส์ $h(t)$ ตามลำดับ

from sympy import *
from IPython.display import display, Math
# define symbols for independent variables
t = Symbol('t', real=True)
s = Symbol('s') 
# define symbols for RC constant
R,C,tau = symbols('R C tau', positive=True)
# define symbols for time-dependent functions
v_C = Function('v_C')(t)
v_S = Function('v_S')(t)
# define symbols for s-dependent functions
V_C = Function('V_C')(s)
V_S = Function('V_S')(s)
LT = lambda y: laplace_transform(y,t,s)
lt_ode = Eq( R*C*LT(v_C.diff(t)) + LT(v_C), LT(v_S) )
LT_V_C = LaplaceTransform(v_C,t,s)
LT_V_S = LaplaceTransform(v_S,t,s)
# substitute symbols
lt_ode = lt_ode.subs( {LT_V_C: V_C, LT_V_S: V_S, R*C: tau} )
# assume v_C(0) = 0 (zero initial condition)
lt_ode = lt_ode.subs( {v_C.subs(t,0): 0} )
# after applying Laplace transform to the ODE, display the final result
display( lt_ode )
# move V_C(s) to the left side of the equation
V_C = solve( lt_ode, V_C )[0]
h = Function('h')(t)
# get the transfer function H(s) = V_C(s)/V_S(s)
H = V_C/V_S
# get the impulse response h(t)
h = inverse_laplace_transform( H, s, t)
# display the transfer function
display( H )
# display the impulse response
display( h )

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง (ต่อจากส่วนแรก)

ในกรณีที่สัญญาณอินพุตเป็นรูปคลื่นไซน์ (Sinusoidal Wave) ที่มีความถี่เชิงมุม $\omega = 2\pi f$ สัญญาณเอาต์พุตจะได้เป็นรูปคลื่นไซน์เช่นกัน เมื่อเข้าสู่สภาวะคงตัว (Steady State)

ถ้าแทนตัวแปร $s$ ด้วยสัญลักษณ์ $j\omega$ (ความถี่เชิงซ้อน) ฟังก์ชัน $H(s)$ จะเป็นฟังก์ชันที่ให้ค่าเป็นเลขเชิงซ้อน ใช้สัญลักษณ์เป็น $H(j\omega)$ ซึ่งขึ้นอยู่ตัวแปรอิสระ $\omega$ และสามารถเขียนให้อยู่ของขนาด (Magnitude) และมุม (Phase) ได้ตามรูปแบบของออยเลอร์ (Euler's Exponential Form) ดังนี้

$$H(j\omega) = |H(j\omega)|\, e^{ j\angle{ H(j\omega) } } \\ $$

ดังนั้นสำหรับวงจร RC ในตัวอย่างนี้จะได้ฟังก์ชันถ่ายโอน

$$\begin{align} H(j\omega) &= \frac{1}{1 + j\omega RC} = \frac{1 - j\omega RC}{1 + (\omega RC)^2} \\ |H(j\omega)| &= \frac{1}{|1 + j\omega RC|} = \frac{1}{\sqrt{1 + (\omega RC)^2}} \\ \angle(H(j\omega)| &= -tan^{-1}(\omega RC) \\ \end{align}$$

 

---

▷ การวิเคราะห์ผลตอบสนองชั่วขณะสำหรับวงจร RC#

ถ้ากำหนดให้แหล่งจ่ายแรงดันไฟฟ้า $v_S(t)$ เป็นอินพุตของวงจร RC และเป็นสัญญาณรูปคลื่นไซน์

$$v_S(t) = V_S\, sin(\omega_0 t), \; V_S > 0 \\ \omega_0 = 2\pi f_0$$

ให้ $v_C(t)$ เป็นเอาต์พุตของวงจรดังกล่าว หรือผลเฉลยของสมการเชิงอนุพันธ์ที่เกี่ยวข้องกับวงจร ลองมาดูตัวอย่างการเขียนโค้ด Python เพื่อวิเคราะห์วงจร และหาผลตอบสนองชั่วขณะในเชิงเวลาของวงจร (Transient Response) โดยกำหนดให้ใช้ค่าตัวเลข $f_0 = 50\,Hz$, $V_S = 5\,V$, $R = 1k\,\Omega$ และ $C = 10\,uF$ สำหรับการแสดงรูปกราฟ $v_C(t)$ และ $v_S(t)$

from sympy import *
from IPython.display import display, Math
# define symbols for independent variables
t = Symbol('t', real=True)
s = Symbol('s') 
# define symbols for RC constant
R,C,tau = symbols('R C tau', positive=True )
V_S,omega_0 = symbols( 'V_S omega_0', positive=True )
# define symbols for time-dependent functions
v_C = Function('v_C')(t)
v_S = V_S*sin(omega_0*t)
# define symbols for s-dependent function
V_C = Function('V_C')(s)
LT = lambda y: laplace_transform(y,t,s)
# apply the Laplace transform to the ODE
lt_ode = Eq( R*C*LT(v_C.diff(t)) + LT(v_C), LT(v_S)[0] )
# display the Laplace tranform of the ODE
display(lt_ode)
# substitute symbols
LT_V_C = LaplaceTransform(v_C,t,s)
LT_V_S = LaplaceTransform(v_S,t,s)
lt_ode = lt_ode.subs( {LT_V_C: V_C, LT_V_S: V_S, R*C: tau} )
# move V_C(s) to the left side of the equation
V_C = solve( lt_ode, V_C )[0]
# display V_C(s) which is the Laplace transform of v_C(t)
display( factor(V_C) )
# assume v_C(0) = 0 (zero initial condition)
V_C = V_C.subs( {v_C.subs(t,0): 0} )
# perform the inverse Laplace transform
v_C = inverse_laplace_transform( V_C, s, t )
# display v_C(t) which is the solution of the ODE
display( collect( expand(v_C), 
    {sin(omega_0*t), exp(-t/tau), Heaviside(t)} ) )

import numpy as np
import matplotlib.pyplot as plt
# substitute numeric values for the circuit parameters
params = {tau: (1e3)*(10e-6), omega_0: 2*pi*50, V_S: 5 }
f1 = lambdify( t, v_S.subs( params ), "numpy" )
f2 = lambdify( t, v_C.subs( params ), "numpy" )
t = np.linspace(0, 0.1, 200)
plt.figure(figsize=(10,6))
plt.plot( t, f1(t), t, f2(t) )
plt.ylabel('$v_C(t)$')
plt.xlabel('t')
plt.grid()
plt.legend(['$v_S(t)$', '$v_C(t)$'])
plt.show()

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง

รูป: เอาต์พุตเมื่อรันโค้ดตัวอย่าง (ต่อจากส่วนแรก)

รูป: ตัวอย่างรูปกราฟสำหรับ $v_S(t)$ และ $v_C(t)$ ในช่วงเวลา $0 .. 0.1\,sec$

 

ถ้าเปลี่ยน $v_S(t)$ จากสัญญาณรูปคลื่นไซน์ ให้เป็นแบบขั้นบันได เช่น

$$v_S(t) = V_S \cdot u(t) = \begin{cases} V_S & \; t \geq 0 \\ 0 & \; t < 0\\ \end{cases}$$ โดยที่ $u(t)$ หมายถึง ฟังก์ชันเฮฟวีไซด์ (Heaviside Function)

หากจะศึกษาผลตอบสนองต่อสัญญาณอินพุตที่เป็นฟังก์ชันแบบขั้นบันได (Step Function) ก็สามารถเขียนโค้ดเพื่อหาผลตอบสนองชั่วขณะสำหรับ $v_C(t)$ ได้ดังนี้ (ในตัวอย่างนี้ ได้กำหนดเงื่อนไขเริ่มต้น $v_C(0) = 0$)

from sympy import *
from IPython.display import display, Math
# define symbols for independent variables
t = Symbol('t', real=True)
s = Symbol('s') 
# define symbols for RC constant
R,C,tau = symbols('R C tau', positive=True )
V_S = symbols( 'V_S', positive=True )
# define symbols for time-dependent functions
v_C = Function('v_C')(t)
v_S = V_S*Heaviside(t) # <---- step function
# define symbols for s-dependent function
V_C = Function('V_C')(s)
LT = lambda y: laplace_transform(y,t,s)
# apply the Laplace transform to the ODE
lt_ode = Eq( R*C*LT(v_C.diff(t)) + LT(v_C), LT(v_S)[0] )
# display the Laplace tranform of the ODE
display(lt_ode)
# substitute symbols
LT_V_C = LaplaceTransform(v_C,t,s)
LT_V_S = LaplaceTransform(v_S,t,s)
lt_ode = lt_ode.subs( {LT_V_C: V_C, LT_V_S: V_S, R*C: tau} )
# move V_C(s) to the left side of the equation
V_C = solve( lt_ode, V_C )[0]
# display V_C(s) which is the Laplace transform of v_C(t)
display( factor(V_C) )
# assume v_C(0) = 0 (zero initial condition)
V_C = V_C.subs( {v_C.subs(t,0): 0} )
# perform the inverse Laplace transform
v_C = inverse_laplace_transform( V_C, s, t )
# display v_C(t) which is the solution of the ODE
display( collect( expand(v_C), 
    {exp(-t/tau), Heaviside(t)} ) )

import numpy as np
import matplotlib.pyplot as plt
# substitute numeric values for the circuit parameters
params = {tau: (1e3)*(10e-6), V_S: 5 }
f1 = lambdify( t, v_S.subs( params ), "numpy" )
f2 = lambdify( t, v_C.subs( params ), "numpy" )
t = np.linspace(-0.02, 0.1, 200)
plt.figure(figsize=(10,6))
plt.plot( t, f1(t), t, f2(t), linewidth=2.0 )
plt.ylabel('$v_C(t)$')
plt.xlabel('t')
plt.grid()
plt.legend(['$v_S(t)$', '$v_C(t)$'])
plt.show()

รูป: เอาต์พุตที่ได้จากการรันโค้ดตัวอย่าง

รูป: ตัวอย่างรูปกราฟสำหรับ $v_S(t)$ และ $v_C(t)$

 

---

▷ กล่าวสรุป#

บทความนี้ได้นำเสนอตัวอย่างการแปลงลาปลาซสำหรับฟังก์ชันคณิตศาสตร์พื้นฐาน การนำไปใช้เพื่อหาผลเฉลยของสมการเชิงอนุพันธ์สามัญเชิงเส้นที่มีสัมประสิทธิ์คงตัว และตัวอย่างการเขียนโค้ดและใช้คำสั่งของ Python / SymPy เพื่อหาผลเฉลยด้วยวิธีการประมวลเชิงสัญลักษณ์ เปรียบเทียบผลลัพธ์ที่ได้ในแต่ละวิธี

 

---

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Created: 2022-04-28 | Last Updated: 2022-04-29