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

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

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

---

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

การวิเคราะห์วงจรไฟฟ้าพื้นฐานที่ประกอบด้วยตัวต้านทาน ตัวเก็บประจุ ตัวเหนี่ยวนำไฟฟ้า หรือเรียกว่า วงจรประเภท 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(t) = 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)$ ที่มีลักษณะเป็นขั้นบันไดหนึ่งหน่วย (Unit-step Function) หรือเรียกว่า ฟังก์ชันเฮฟวีไซด์ (Heaviside Function) มีค่าแบ่งเป็นช่วงดังนี้

$$u(t) = \begin{cases} 1, & \; t \geq 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)$

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

 

---

▷ การเขียนโค้ดด้วย MATLAB เพื่อหาผลการแปลงลาปลาซด้านเดียว#

จากตัวอย่างฟังก์ชันพื้นฐาน ลองมาดูการใช้คำสั่ง laplace() ของ MATLAB - Symbolic Math Toolbox เพื่อหาผลการแปลงลาปลาซด้านเดียว โดยมีการประกาศให้ $t$ และ $s$ เป็นสัญลักษณ์สำหรับตัวแปรอิสระของฟังก์ชัน

clearvars
syms t s a c
syms omega_0 t_0 positive
% Laplace transform of the Dirac-delta function
laplace( dirac(t),t,s )
% Laplace transform of a constant function
laplace( c,t,s )
% Laplace transform of the Heaviside function
laplace( heaviside(t),t,s )
% Laplace transform of exponential function
laplace( exp(-a*t),t,s )
% Laplace transform of a sine function
laplace( sin(omega_0*t),t,s )
% Laplace transform of a cosine function
laplace( cos(omega_0*t),t,s )
% Laplace transform of a pulse function
laplace( heaviside(t) - heaviside(t-t_0),t,s )
% Laplace transform of a ramp function
laplace( t*heaviside(t),t,s )
% Laplace transform of an exponentially damped sine
laplace( exp(-a*t)*sin(omega_0*t),t,s)
% Laplace transform of an exponentially damped cosine
laplace( exp(-a*t)*cos(omega_0*t),t,s )

รูป: ตัวอย่างการเขียนโค้ดเพื่อหาผลการแปลงลาปลาซโดยใช้ MATLAB (R2022a)

หากต้องการหาผลการแปลงลาปลาซผกผัน ก็สามารถใช้คำสั่ง ilaplace() ของ MATLAB ในการคำนวณได้

รูป: ตัวอย่างการใช้คำสั่ง ilaplace() เพื่อหาผลการแปลงลาปลาซผกผัน

 

---

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

เราจะพบว่า ผลการแปลงลาปลาซมักอยู่ในรูปของเศษส่วนย่อย (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) \\ \end{align}$$

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

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

 

MATLAB ก็มีคำสั่ง partfrac() และ residue() สำหรับการแยกเป็นเศษส่วนย่อย ลองมาดูตัวอย่างการใช้คำสั่ง partfrac() เพื่อแยกฟังก์ชันเศษส่วนออกเป็นเศษส่วนย่อยแล้วหาผลการแปลงลาปลาซผกผันของแต่ละพจน์ย่อย

clearvars
syms t s f(t) F(s)
F(s) = 5*(s-1)/((s^2 + 4)*(s+1))
% find the inverse Laplace transform of F(s)
f(t) = ilaplace( F(s),s,t )
%  perform partial fraction expansion of F(s)
P = partfrac( F(s), s )
% find inverse Laplace transform of each term of P
f(t) = 0;
for p = children(P)
    f(t) = f(t) + ilaplace( p, s, t );
end
% show the result
disp( f(t) )

รูป: ตัวอย่างการใช้คำสั่ง partfrac()

 

---

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

ถัดไปเป็นตัวอย่างของสมการ 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}$$

 

ตัวอย่างการเขียนโค้ด MATLAB มีดังนี้

clearvars
syms s t y(t) Y(s)
% derivatives of y(t)
D1y = diff(y(t),t,1);
D2y = diff(y(t),t,2);
%  formulate the ODE
ode = D2y + 2*D1y + y(t) == sin(2*t);
% Solution 1: solve the ODE using the dsolve() command
sol = dsolve( ode, [y(0)==-1, subs(D1y,t,0)==1] )
% Laplace transform of the ODE
lt_ode = laplace( ode, t, s );
lt_ode = subs( lt_ode, laplace(y,t,s), Y(s) );
collect(lt_ode,Y(s))
lt_ode = subs( lt_ode, [y(0), subs(D1y,t,0)], [-1,1]);
collect(lt_ode,Y(s))
Y(s) = isolate( lt_ode, Y(s) )
Y(s) = partfrac( Y(s) )
% Solution 2: perform the inverse Laplace transform of Y(s)
y(t) = ilaplace( rhs(Y(s)),s, t )
% compare the solutions 
isAlways( y(t) == sol ) 

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

รูป: ตัวอย่างการทำคำสั่ง MATLAB (ส่วนที่หนึ่ง)

รูป: ตัวอย่างการทำคำสั่ง MATLAB (ส่วนที่สอง)

 

---

▷ การวิเคราะห์วงจร 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}$$

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

clearvars
syms s t v_C(t) v_S(t) V_C(s) V_S(s) H(s) h(t)
syms R C tau positive 
% the first-order ODE for the RC filter circuit
ode = R*C*diff( v_C(t), t) + v_C(t) == v_S(t);
% take the Laplace transform to the ODE
lt_ode = laplace( ode,t,s );
% substitute Laplace transforms of v_C(t) and v_S(t) with V_C(s) and V_S(s)
lt_ode = subs( lt_ode, laplace(v_C(t),t,s), V_C(s) );
lt_ode = subs( lt_ode, laplace(v_S(t),t,s), V_S(s) );
% rewrite the equation so that V_C(s) appears on the left side 
V_C(s) = isolate( lt_ode, V_C(s) );
% substitute RC with tau and V_C(0) with 0 (zero initial condition)
V_C(s) = subs( V_C(s), [R*C, v_C(0)], [tau, 0] )
% find the transfer function H(s)
H(s) = rhs(V_C(s))/ V_S(s)
% find the impulse response h(t)
h(t) = ilaplace( H(s),s,t )

รูป: ตัวอย่างเอาต์พุตเมื่อทำคำสั่ง MATLAB

ในกรณีที่สัญญาณอินพุตเป็นรูปคลื่นไซน์ (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}$$

 

---

▷ การแสดงรูปกราฟ Bode Plot สำหรับวงจร RC#

เราสามารถใช้คำสั่งของ MATLAB เช่น tf() เพื่อสร้างฟังก์ชันถ่ายโอนที่อยู่ในรูปเศษส่วนพหุนาม ของตัวแปร $s$ โดยจะต้องระบุค่าสัมประสิทธิ์ของเศษส่วนพหุนามสำหรับ $H(s) = N(s)/D(s)$ แล้วจึงใช้คำสั่ง bodeplot() แสดงรูปกราฟของฟังก์ชันถ่ายโอน (เรียกว่า Bode Plot) โดยแบ่งเป็นขนาด (Magnitude Plot) และเฟส (Phase Plot) ที่ขึ้นอยู่กับค่าความถี่

clearvars;
RC = (1e+3)*(10e-6); % use R=1kOhm and C=10uF
% define the transfer function H(s)
% H(s) = N(s)/D(s) = [1/RC] / [s + 1/RC]
H = tf(1/RC,[1 1/RC])
% set Bode plot options
opts = bodeoptions;
opts.FreqUnits = 'Hz';
opts.FreqScale = 'log';
opts.MagUnits  = 'dB';
opts.MagScale  = 'linear';
opts.PhaseUnits = 'deg';
% set the frequency range for Bode plot
freq_range = {0.1 10^6};
% show Bode plot 
bodeplot( H, freq_range, opts ); grid on;

รูป: ตัวอย่างคำสั่ง MATLAB สำหรับการแสดงรูปกราฟ Bode Plot

รูป: รูปกราฟ Bode Plot สำหรับฟังก์ชันถ่ายโอนของวงจร RC (R=1kΩ, C=10uF)

จากรูปกราฟจะเห็นได้ว่า เมื่อความถี่ $f$ เพิ่มขึ้น จะทำให้ฟังก์ชันถ่ายโอน $H(j\omega)$ มีขนาดลดลง ซึ่งหมายความว่า อัตราส่วนระหว่างเอาต์พุตต่ออินพุต จะมีขนาดลดลงเมื่อความถี่เพิ่มขึ้น

ถ้ามีความถี่ใกล้ $0$ จะทำให้ $|H(j\omega)|$ มีค่าใกล้เคียง $1$ (หรือ $0 \,dB$) แต่จะลดลงเมื่อความถี่เพิ่มขึ้น ดังนั้นวงจร RC ในรูปแบบนี้จึงทำหน้าที่เป็นตัวกรองให้ความถี่ต่ำผ่านได้ดี

เฟส $\angle H(j\omega)$ จะเริ่มลดลงจาก $0^{\circ}$ เมื่อความถี่เพิ่มขึ้น ซึ่งหมายความว่า สัญญาณเอาต์พุตจะมีความต่างเฟสเพิ่มมาขึ้นเมื่อเปรียบเทียบกับสัญญาณอินพุตถ้าความถี่เพิ่มมากขึ้น

ถ้าไม่ใช้คำสั่ง tf() และ bodeplot() ก็มีตัวอย่างการเขียนโค้ด MATLAB ดังนี้

clearvars; clf;
syms s f omega
syms R C positive
% define the transfer function H(s) of the RC circuit
H = 1/(1+s*R*C)
% substitute s with j*omega 
H = subs( H, s, 1i*omega ) 
% substitute omega with 2*pi*f
% substitute R and C with numeric values
H = subs( H, [omega,R,C], [2*pi*f, 1e+3, 10e-6] );
% convert symbolic functions to numeric functions
H_mag   = matlabFunction( 20*log10( abs(H) ) );
H_phase = matlabFunction( (180/pi)*angle( H ) );
% define the frequency range (log-scale): 10^0 .. 10^4
f = logspace( 0, 4 );
% plot the magnitude of H(j*omega)
semilogx( f, H_mag(f) ), grid on, 
xlabel('Hz'), ylabel('Magnitude (dB)');
% plot the phase of H(j*omega)
semilogx( f, H_phase(f) ), grid on,
xlabel('Hz'), ylabel('Phase (deg)');

รูป: ตัวอย่างคำสั่ง MATLAB สำหรับการแสดงรูปกราฟ Bode Plot

รูป: รูปกราฟสำหรับขนาดของ $H(j\omega)$ สำหรับความถี่ในช่วง $1..10^4$ (log-scale)

รูป: รูปกราฟสำหรับเฟสของ $H(j\omega)$ สำหรับความถี่ในช่วง $1..10^4$ Hz (log-scale)

 

สำหรับวงจร RC ที่เป็นวงจรตัวกรองให้ความถี่ต่ำผ่านได้ดี (Low-pass RC Filter) ความถี่ $\omega=\omega_c$ ที่ทำให้ $|H(j\omega_c)| = \frac{1}{\sqrt{2}}$ หรือ $20\cdot log_{10}(|H(j\omega_c)|) = -3\,dB\;$ เรียกว่า Cutoff Frequency หรือ Corner Frequency และสามารถคำนวณได้ดังนี้

$$|H(j\omega_c)|^2 = \frac{1}{2} = \frac{1}{1 + (\omega_c RC)^2} \\ \Rightarrow \quad \omega_c = 2\pi f_c = \frac{1}{RC} \\ f_c = \frac{1}{2\pi RC} \\ $$

ถ้าให้ $R=1k \Omega$ และ $C=10uF$ จะได้ $f_c = \frac{1}{2\pi (0.01)} = 15.915 \,Hz$ สำหรับวงจรตัวอย่าง

 

---

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

clearvars
syms s t v_C(t) v_S(t) V_C(s)
syms R C V_S omega_0 tau positive 
% The input is a sinusoidal voltage source.
v_S(t) = V_S*sin(omega_0*t);
% define the ODE for the RC filter
ode = R*C*diff( v_C(t), t) + v_C(t) == v_S(t);
% perform LT of the ODE
lt_ode = laplace( ode,t,s );
% substitute a symbol for V_C(s)
lt_ode = subs( lt_ode, laplace(v_C(t),t,s), V_C(s) );
% rewrite the equation so that V_C(s) appears on the left-hand side
sol = isolate( lt_ode, V_C(s) );
% substitute R*C with tau (the time constant symbol)
sol = subs( sol, R*C, tau );
% show V_C(s)
disp( sol )
% assume zero initial condition: v_C(0) = 0
sol = subs( sol, v_C(0), 0 );
% find the inverse Laplace  transform to get v_C(t)
sol = ilaplace( sol, s, t);
% get the output signal v_C(t)
v_C(t) = rhs(sol)
% substitute tau with a specific value (R=1k, C=10u)
v_C(t) = subs( v_C(t), tau, 1e3 * 10e-6 );
% substitute omega_0 with 2*pi*50Hz and V_S with 5V
v_C(t) = subs( v_C(t), [omega_0, V_S], [2*pi*50, 5] );
v_S(t) = subs( v_S(t), [omega_0, V_S], [2*pi*50, 5] );
fplot( [v_S(t), v_C(t)], [0,0.1] ), 
grid on, xlabel('t'), legend( ['v_S(t)';'v_C(t)'] )

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

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

 

ถ้าเปลี่ยน $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) ก็สามารถเขียนโค้ด MATLAB เพื่อหาผลตอบสนองชั่วขณะสำหรับ $v_C(t)$ ได้ดังนี้ (ในตัวอย่างนี้ ได้กำหนดเงื่อนไขเริ่มต้น $v_C(0) = 0$)

clearvars; clf;
syms s t v_C(t) v_S(t) V_C(s)
syms R C V_S tau positive
% The input is a DC source.
v_S(t) = V_S*heaviside(t);
% define the ODE for the RC filter
ode = R*C*diff( v_C(t), t) + v_C(t) == v_S(t);
% perform LT of the ODE
lt_ode = laplace( ode,t,s );
% substitute a symbol for V_C(s)
lt_ode = subs( lt_ode, laplace(v_C(t),t,s), V_C(s) );
% rewrite the equation so that V_C(s) appears on the left-hand side
sol = isolate( lt_ode, V_C(s) );
% substitute R*C with tau (the time constant symbol)
sol = subs( sol, R*C, tau );
% show V_C(s)
disp( sol )
% assume zero initial condition: v_C(0) = 0
sol = subs( sol, v_C(0), 0 );
% find the inverse Laplace transform to get v_C(t)
sol = ilaplace( sol, s, t);
% get the output signal v_C(t), for t>=0
v_C(t) = rhs(sol)*heaviside(t)
% substitute tau with a specific value (R=1k, C=10u)
v_C(t) = subs( v_C(t), tau, 1e3 * 10e-6 );
% substitute V_S with 5V
v_C(t) = subs( v_C(t), V_S, 5 );
v_S(t) = subs( v_S(t), V_S, 5 );
fplot( [v_S(t), v_C(t)], 'linewidth', 1.25 ), 
grid on, xlabel('t'), 
xlim([-0.01 0.1]), ylim([-0.2 5.2]), 
legend( ['v_S(t)';'v_C(t)'] )

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

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

 

---

▷ กล่าวสรุป#

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

 

---

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

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