Engineering Mathematics 1: Ordinary Differential Equations (ODE)

**Engineering mathematics** is an unavoidable subject for any engineering student. In particular, Ordinary Differential Equations (ODEs) form the mathematical foundation of virtually every engineering discipline — from circuit analysis and control system design to mechanical vibration analysis. This article systematically explains ODEs from the basics to Laplace transforms, centered on engineering application examples.

1. Fundamentals of Differential Equations

What is a Differential Equation?

A **differential equation** is an equation that expresses a relationship between an unknown function and its derivatives. For example:

$$\frac{dy}{dx} = 2x, \quad y'' + 3y' + 2y = 0, \quad \frac{\partial u}{\partial t} = k\frac{\partial^2 u}{\partial x^2}$$

Differential equations are a core tool for mathematically modeling natural phenomena and engineering systems.

Classification: ODE vs PDE

**Ordinary Differential Equation (ODE)**: Contains derivatives with respect to only one independent variable.

$$\frac{d^2x}{dt^2} + 2\frac{dx}{dt} + x = 0$$

**Partial Differential Equation (PDE)**: Contains partial derivatives with respect to two or more independent variables.

$$\frac{\partial u}{\partial t} = c^2\frac{\partial^2 u}{\partial x^2}$$

This article focuses on ODEs.

Order and Degree

The **order** is the highest derivative appearing in the equation.

- First-order ODE: $y' + y = x$ (highest derivative is $y'$)

- Second-order ODE: $y'' + 2y' + y = 0$ (highest derivative is $y''$)

The **degree** is the exponent of the highest-order derivative. In most engineering problems, the degree is 1.

Linear vs Nonlinear

**Linear ODE**: The unknown function and its derivatives appear only as first-degree terms.

$$a_n(x)y^{(n)} + a_{n-1}(x)y^{(n-1)} + \cdots + a_1(x)y' + a_0(x)y = f(x)$$

**Nonlinear ODE**: Contains nonlinear terms such as $y^2$, $y \cdot y'$, $\sin(y)$, etc.

$$y'' + \sin(y) = 0 \quad \text{(pendulum equation)}$$

General Solution, Particular Solution, Singular Solution

- **General solution**: A complete solution containing arbitrary constants

- **Particular solution**: A solution obtained by applying initial conditions to determine the constants

- **Singular solution**: A special solution that cannot be obtained from the general solution

ODEs in Engineering

Most engineering systems are described by ODEs.

**RC charging circuit**:

$$RC\frac{dv}{dt} + v = V_s$$

**Spring-mass-damper system**:

$$m\ddot{x} + c\dot{x} + kx = F(t)$$

**Population growth model**:

$$\frac{dP}{dt} = kP$$

**RL circuit**:

$$L\frac{di}{dt} + Ri = E(t)$$

All of these are ODEs. Learning how to solve them allows you to fully analyze the dynamic behavior of engineering systems.

2. First-Order Ordinary Differential Equations

2.1 Separation of Variables

**Form**: $\frac{dy}{dx} = f(x)g(y)$

**Solution procedure**:

$$\frac{dy}{g(y)} = f(x)dx$$

Integrate both sides:

$$\int \frac{dy}{g(y)} = \int f(x)dx + C$$

**Example: RC Charging Circuit**

Circuit equation:

$$\frac{dv}{dt} = \frac{V_s - v}{RC}$$

Separation of variables:

$$\frac{dv}{V_s - v} = \frac{dt}{RC}$$

Integration:

$$-\ln|V_s - v| = \frac{t}{RC} + C_1$$

Applying initial condition $v(0) = 0$:

$$V_s - v = V_s e^{-t/(RC)}$$

$$\boxed{v(t) = V_s\left(1 - e^{-t/(RC)}\right)}$$

This result represents the **transient response** of an RC circuit. The voltage converges to $V_s$ at a rate determined by the time constant $\tau = RC$.

**Physical interpretation**:

- At $t = \tau$: $v \approx 0.632 V_s$ (63.2% charged)

- At $t = 5\tau$: $v \approx 0.993 V_s$ (almost fully charged)

2.2 Integrating Factor Method

**Standard form of a first-order linear ODE**:

$$\frac{dy}{dx} + P(x)y = Q(x)$$

**Integrating factor**: $\mu(x) = e^{\int P(x)dx}$

**Derivation**:

Multiply both sides by $\mu(x)$:

$$\mu(x)\frac{dy}{dx} + \mu(x)P(x)y = \mu(x)Q(x)$$

The left side, by the product rule:

$$\frac{d}{dx}[\mu(x)y] = \mu(x)Q(x)$$

Integrating:

$$\mu(x)y = \int \mu(x)Q(x)dx + C$$

$$\boxed{y = \frac{1}{\mu(x)}\left[\int \mu(x)Q(x)dx + C\right]}$$

**Example: RL Circuit**

$$L\frac{di}{dt} + Ri = E$$

Converting to standard form:

$$\frac{di}{dt} + \frac{R}{L}i = \frac{E}{L}$$

Integrating factor: $\mu(t) = e^{(R/L)t}$

$$\frac{d}{dt}\left[e^{(R/L)t} \cdot i\right] = \frac{E}{L}e^{(R/L)t}$$

Integrating:

$$e^{(R/L)t} \cdot i = \frac{E}{R}e^{(R/L)t} + C$$

Applying initial condition $i(0) = 0$:

$$\boxed{i(t) = \frac{E}{R}\left(1 - e^{-(R/L)t}\right)}$$

The current converges to steady state $E/R$ with time constant $\tau = L/R$.

2.3 Exact ODEs

**Form**: $M(x,y)dx + N(x,y)dy = 0$

**Exactness condition**:

$$\frac{\partial M}{\partial y} = \frac{\partial N}{\partial x}$$

When this condition is satisfied, there exists a potential function $F$ such that $F(x,y) = C$.

**Solution procedure**:

1. Integrate $\frac{\partial F}{\partial x} = M$ with respect to $x$ to find $F$

2. Use the condition $\frac{\partial F}{\partial y} = N$ to determine the integration constant

3. General solution: $F(x,y) = C$

**Example**:

$(2xy + y^2)dx + (x^2 + 2xy)dy = 0$

Checking exactness:

$$\frac{\partial M}{\partial y} = 2x + 2y = \frac{\partial N}{\partial x}$$

Since it is exact:

$$F = \int M\,dx = \int (2xy + y^2)dx = x^2y + xy^2 + g(y)$$

Since $\frac{\partial F}{\partial y} = x^2 + 2xy + g'(y) = N = x^2 + 2xy$, we have $g'(y) = 0$, $g(y) = C$

General solution: $x^2y + xy^2 = C$

2.4 Bernoulli Equation

**Form**:

$$\frac{dy}{dx} + P(x)y = Q(x)y^n \quad (n \neq 0, 1)$$

**Substitution**: Setting $v = y^{1-n}$ transforms this into a linear ODE.

Since $v' = (1-n)y^{-n}y'$, dividing the original equation by $y^n$:

$$y^{-n}\frac{dy}{dx} + P(x)y^{1-n} = Q(x)$$

$$\frac{1}{1-n}\frac{dv}{dx} + P(x)v = Q(x)$$

$$\frac{dv}{dx} + (1-n)P(x)v = (1-n)Q(x)$$

This is a first-order linear ODE in $v$, solvable by the integrating factor method.

**Application**: The logistic growth model $\frac{dP}{dt} = rP\left(1 - \frac{P}{K}\right)$ is a Bernoulli equation with $n=2$.

3. Second-Order Linear Ordinary Differential Equations

3.1 Homogeneous Equations

**Standard form**:

$$ay'' + by' + cy = 0$$

**Characteristic equation** derivation: Substitute $y = e^{rx}$.

$$a r^2 e^{rx} + b r e^{rx} + c e^{rx} = 0$$

$$\boxed{ar^2 + br + c = 0}$$

Three cases arise depending on the sign of the discriminant $\Delta = b^2 - 4ac$.

**Case 1: Two distinct real roots** ($\Delta > 0$)

For $r_1 \neq r_2$:

$$y = C_1 e^{r_1 x} + C_2 e^{r_2 x}$$

**Case 2: Repeated real root** ($\Delta = 0$)

For $r_1 = r_2 = r = -b/(2a)$:

$$y = (C_1 + C_2 x)e^{rx}$$

The term $C_2 x e^{rx}$ must be added so that the **Wronskian** is nonzero, constructing an independent second solution.

**Case 3: Complex conjugate roots** ($\Delta < 0$)

For $r = \alpha \pm i\beta$ (where $\alpha = -b/(2a)$, $\beta = \sqrt{4ac - b^2}/(2a)$):

$$y = e^{\alpha x}(C_1\cos\beta x + C_2\sin\beta x)$$

This form is converted to real solutions using Euler's formula $e^{i\beta x} = \cos\beta x + i\sin\beta x$.

3.2 Complete Analysis of RLC Circuits

**Circuit equation** (in terms of charge $q(t)$):

$$L\frac{d^2q}{dt^2} + R\frac{dq}{dt} + \frac{q}{C} = 0$$

In standard form:

$$q'' + \frac{R}{L}q' + \frac{1}{LC}q = 0$$

**Characteristic equation**:

$$r^2 + \frac{R}{L}r + \frac{1}{LC} = 0$$

$$r = -\frac{R}{2L} \pm \sqrt{\left(\frac{R}{2L}\right)^2 - \frac{1}{LC}}$$

Defining the damping coefficient $\alpha = R/(2L)$ and resonant frequency $\omega_0 = 1/\sqrt{LC}$:

$$r = -\alpha \pm \sqrt{\alpha^2 - \omega_0^2}$$

**Overdamped**: $\alpha > \omega_0$ ($R > 2\sqrt{L/C}$)

$$q(t) = C_1 e^{r_1 t} + C_2 e^{r_2 t} \quad (r_1, r_2 < 0)$$

A sum of two exponential terms — decays monotonically without oscillation.

**Critically Damped**: $\alpha = \omega_0$ ($R = 2\sqrt{L/C}$)

$$q(t) = (C_1 + C_2 t)e^{-\alpha t}$$

Returns to equilibrium most rapidly without oscillation. Very important in control system design.

**Underdamped**: $\alpha < \omega_0$ ($R < 2\sqrt{L/C}$)

Defining the damped oscillation frequency $\omega_d = \sqrt{\omega_0^2 - \alpha^2}$:

$$q(t) = e^{-\alpha t}(C_1\cos\omega_d t + C_2\sin\omega_d t)$$

Oscillates while decaying. Commonly seen in communication circuits and vibration systems.

**Summary of physical interpretation**:

| Condition | Name | Response Characteristics |

| ----------------- | ----------------- | --------------------------------- |

| $R > 2\sqrt{L/C}$ | Overdamped | Slow decay without oscillation |

| $R = 2\sqrt{L/C}$ | Critically damped | Fastest decay without oscillation |

| $R < 2\sqrt{L/C}$ | Underdamped | Oscillates and decays |

| $R = 0$ | Undamped | Continuous oscillation |

3.3 Non-homogeneous Equations

**Form**:

$$ay'' + by' + cy = f(x)$$

**Complete solution**: $y = y_h + y_p$

- $y_h$: General solution of the homogeneous equation (complementary function)

- $y_p$: Particular solution

Method of Undetermined Coefficients

Guess the form of $y_p$ based on the form of $f(x)$.

| Form of $f(x)$ | Guessed form of $y_p$ |

| ---------------------------------- | --------------------------------------- |

| $ke^{ax}$ | $Ae^{ax}$ |

| $kx^n$ | $A_n x^n + \cdots + A_1 x + A_0$ |

| $k\cos\omega x$ or $k\sin\omega x$ | $A\cos\omega x + B\sin\omega x$ |

| $ke^{ax}\cos\omega x$ | $e^{ax}(A\cos\omega x + B\sin\omega x)$ |

**Note**: If the guessed form of $y_p$ overlaps with a term in $y_h$, multiply by $x$.

**Example**: $y'' + 4y = 3\cos 2x$

Homogeneous solution: $y_h = C_1\cos 2x + C_2\sin 2x$

Since the guessed form of $y_p$ overlaps with $y_h$, modify:

$$y_p = x(A\cos 2x + B\sin 2x)$$

After substitution: $A = 0$, $B = 3/4$

$$y_p = \frac{3}{4}x\sin 2x$$

This is **resonance**, where the amplitude increases proportionally with time.

**Forced vibration response**:

$$m\ddot{x} + c\dot{x} + kx = F_0\cos\omega t$$

Steady-state particular solution:

$$x_p(t) = \frac{F_0/k}{\sqrt{(1-r^2)^2 + (2\zeta r)^2}}\cos(\omega t - \phi)$$

Where frequency ratio $r = \omega/\omega_n$ and damping ratio $\zeta = c/(2\sqrt{mk})$.

The amplitude is maximum at $r = 1$ (resonance condition).

Variation of Parameters

Used for general $f(x)$ where the method of undetermined coefficients cannot be applied.

Given $y_h = C_1 y_1 + C_2 y_2$:

$$y_p = -y_1 \int \frac{y_2 f}{W}dx + y_2 \int \frac{y_1 f}{W}dx$$

Wronskian: $W = y_1 y_2' - y_1' y_2$

This formula applies to any continuous function $f(x)$.

4. Systems of Differential Equations

4.1 Matrix Representation

A system of $n$ first-order ODEs:

$$\frac{d\mathbf{X}}{dt} = A\mathbf{X}$$

Where $\mathbf{X} = [x_1, x_2, \ldots, x_n]^T$ and $A$ is an $n \times n$ coefficient matrix.

4.2 Eigenvalue/Eigenvector Solution

Form of solution: $\mathbf{X} = \mathbf{v}e^{\lambda t}$

Substituting: $\lambda \mathbf{v} = A\mathbf{v}$

So $\lambda$ is an **eigenvalue** of $A$, and $\mathbf{v}$ is the corresponding **eigenvector**.

Characteristic equation: $\det(A - \lambda I) = 0$

**Example**: Two connected tanks system

$$\frac{d}{dt}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} -2 & 1 \\ 1 & -2 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}$$

Characteristic equation: $(\lambda+2)^2 - 1 = 0$

$$\lambda_1 = -1, \quad \lambda_2 = -3$$

Eigenvectors: $\mathbf{v}_1 = [1, 1]^T$, $\mathbf{v}_2 = [1, -1]^T$

General solution:

$$\mathbf{X} = C_1 \begin{pmatrix} 1 \\ 1 \end{pmatrix} e^{-t} + C_2 \begin{pmatrix} 1 \\ -1 \end{pmatrix} e^{-3t}$$

4.3 Application to Circuit Network Analysis

Circuit with two loops (Kirchhoff's voltage law):

$$L_1\frac{di_1}{dt} + R_1 i_1 - M\frac{di_2}{dt} = E_1(t)$$

$$L_2\frac{di_2}{dt} + R_2 i_2 - M\frac{di_1}{dt} = E_2(t)$$

Where $M$ is the mutual inductance. This can be written in matrix form and solved using the eigenvalue method or Laplace transform.

5. Laplace Transform

5.1 Definition and Convergence

**Definition**:

$$\mathcal{L}\{f(t)\} = F(s) = \int_0^\infty f(t)e^{-st}dt$$

$s = \sigma + j\omega$ is a complex variable. For convergence, $\text{Re}(s) > \sigma_c$ (half-plane of convergence) must be satisfied.

**Why use the Laplace transform?**

- Differential operations are converted to algebraic operations

- Initial conditions are automatically incorporated

- Complex ODEs can be solved as simple algebraic equations

- It forms the basis of the transfer function concept

5.2 Key Transform Pairs

| $f(t)$ | $F(s) = \mathcal{L}\{f(t)\}$ |

| ----------------------- | ---------------------------------- |

| $1$ (unit step) | $\dfrac{1}{s}$ |

| $t$ | $\dfrac{1}{s^2}$ |

| $t^n$ | $\dfrac{n!}{s^{n+1}}$ |

| $e^{at}$ | $\dfrac{1}{s-a}$ |

| $\sin\omega t$ | $\dfrac{\omega}{s^2+\omega^2}$ |

| $\cos\omega t$ | $\dfrac{s}{s^2+\omega^2}$ |

| $e^{at}\sin\omega t$ | $\dfrac{\omega}{(s-a)^2+\omega^2}$ |

| $e^{at}\cos\omega t$ | $\dfrac{s-a}{(s-a)^2+\omega^2}$ |

| $\delta(t)$ (impulse) | $1$ |

| $u(t-a)$ (delayed step) | $\dfrac{e^{-as}}{s}$ |

**Derivation example**: $\mathcal{L}\{e^{at}\}$

$$\mathcal{L}\{e^{at}\} = \int_0^\infty e^{at} e^{-st} dt = \int_0^\infty e^{-(s-a)t} dt = \left[\frac{e^{-(s-a)t}}{-(s-a)}\right]_0^\infty = \frac{1}{s-a}$$

Converges when $s > a$.

5.3 Key Properties

**Linearity**:

$$\mathcal{L}\{af(t) + bg(t)\} = aF(s) + bG(s)$$

**Time shift (First shifting theorem)**:

$$\mathcal{L}\{e^{at}f(t)\} = F(s-a)$$

**Frequency shift (Second shifting theorem)**:

$$\mathcal{L}\{f(t-a)u(t-a)\} = e^{-as}F(s)$$

**Differentiation transform** (most important property):

$$\mathcal{L}\{f'(t)\} = sF(s) - f(0)$$

$$\mathcal{L}\{f''(t)\} = s^2F(s) - sf(0) - f'(0)$$

$$\mathcal{L}\{f^{(n)}(t)\} = s^nF(s) - s^{n-1}f(0) - s^{n-2}f'(0) - \cdots - f^{(n-1)}(0)$$

Initial conditions are automatically included in the s-domain.

**Integration transform**:

$$\mathcal{L}\left\{\int_0^t f(\tau)d\tau\right\} = \frac{F(s)}{s}$$

**Convolution theorem**:

$$\mathcal{L}\{(f*g)(t)\} = F(s)G(s)$$

Where $(f*g)(t) = \int_0^t f(\tau)g(t-\tau)d\tau$

This theorem mathematically proves that the output of a system is the convolution of the impulse response and the input signal.

**Final value theorem**:

$$\lim_{t\to\infty} f(t) = \lim_{s\to 0} sF(s)$$

Allows determination of the steady-state value without solving in the t-domain.

**Initial value theorem**:

$$\lim_{t\to 0^+} f(t) = \lim_{s\to\infty} sF(s)$$

5.4 Inverse Laplace Transform

**Partial Fraction Decomposition**:

For $F(s) = P(s)/Q(s)$, factor the denominator and decompose into a sum of table transform pairs.

**Case 1: Distinct real poles**

$$F(s) = \frac{N(s)}{(s-p_1)(s-p_2)\cdots(s-p_n)} = \frac{A_1}{s-p_1} + \frac{A_2}{s-p_2} + \cdots$$

$$A_k = \lim_{s\to p_k}(s-p_k)F(s)$$

**Case 2: Repeated poles** ($s = p$ is an $m$-fold root)

$$F(s) = \cdots + \frac{A_m}{(s-p)^m} + \frac{A_{m-1}}{(s-p)^{m-1}} + \cdots + \frac{A_1}{s-p}$$

$$A_k = \frac{1}{(m-k)!}\left[\frac{d^{m-k}}{ds^{m-k}}(s-p)^m F(s)\right]_{s=p}$$

**Case 3: Complex conjugate poles**

$$\frac{As + B}{s^2 + 2\alpha s + (\alpha^2 + \beta^2)} \Rightarrow e^{-\alpha t}(C_1\cos\beta t + C_2\sin\beta t)$$

**Heaviside Expansion Theorem**:

When the denominator is a product of distinct linear factors:

$$F(s) = \sum_{k=1}^n \frac{N(p_k)}{Q'(p_k)} \cdot \frac{1}{s - p_k}$$

5.5 Solving ODEs with the Laplace Transform

**Example: RLC Series Circuit Transient Response**

$$L\frac{d^2i}{dt^2} + R\frac{di}{dt} + \frac{i}{C} = E_0\delta(t)$$

Initial conditions: $i(0) = 0$, $i'(0) = 0$ (no charge on capacitor)

**Solution**:

Applying the Laplace transform to both sides:

$$I(s) = \frac{E_0/L}{s^2 + (R/L)s + 1/(LC)}$$

Defining $\alpha = R/(2L)$, $\omega_0 = 1/\sqrt{LC}$, $\omega_d = \sqrt{\omega_0^2 - \alpha^2}$ (underdamped case):

$$I(s) = \frac{E_0/L}{(s+\alpha)^2 + \omega_d^2}$$

Inverse Laplace transform (from table: $\mathcal{L}\{e^{-\alpha t}\sin\omega_d t\} = \omega_d/[(s+\alpha)^2+\omega_d^2]$):

$$i(t) = \frac{E_0}{L\omega_d}e^{-\alpha t}\sin\omega_d t$$

This is the **impulse response** of the RLC circuit, which is also its **Green's function**.

6. Numerical ODE Solutions with Python

6.1 Using scipy.integrate.solve_ivp

from scipy.integrate import solve_ivp

RLC circuit simulation (2nd order ODE -> system of 1st order ODEs)

L * q'' + R * q' + q/C = Vs (step input)

State variables: y[0] = q (charge), y[1] = i = dq/dt (current)

def rlc_circuit(t, y, R, L, C, Vs):

q, dq_dt = y

q'' = (Vs - R*i - q/C) / L

d2q_dt2 = (Vs - R * dq_dt - q / C) / L

return [dq_dt, d2q_dt2]

Circuit parameters

R = 100 # Resistance (Ohm)

L = 0.1 # Inductance (H)

C = 1e-6 # Capacitance (F)

Vs = 10.0 # Supply voltage (V)

Calculate resonant frequency and damping ratio

omega_0 = 1.0 / np.sqrt(L * C)

alpha = R / (2 * L)

zeta = alpha / omega_0

print(f"Resonant frequency: {omega_0/(2*np.pi):.1f} Hz")

print(f"Damping ratio zeta: {zeta:.3f}")

if zeta > 1:

print("-> Overdamped")

elif zeta == 1:

print("-> Critically damped")

else:

print("-> Underdamped")

Numerical solution

y0 = [0.0, 0.0] # Initial conditions: q(0) = 0, i(0) = 0

t_span = (0, 0.01)

t_eval = np.linspace(0, 0.01, 2000)

sol = solve_ivp(

rlc_circuit, t_span, y0,

t_eval=t_eval,

args=(R, L, C, Vs),

method='RK45',

rtol=1e-8,

atol=1e-10

)

Voltage calculation: v_C = q/C

v_C = sol.y[0] / C

i = sol.y[1]

Visualization

fig, axes = plt.subplots(2, 1, figsize=(10, 8))

axes[0].plot(sol.t * 1e3, v_C, 'b-', linewidth=2, label='Capacitor voltage')

axes[0].axhline(y=Vs, color='r', linestyle='--', label=f'Steady state {Vs}V')

axes[0].set_xlabel('Time (ms)')

axes[0].set_ylabel('Voltage (V)')

axes[0].set_title('RLC Series Circuit - Step Response')

axes[0].legend()

axes[0].grid(True, alpha=0.3)

axes[1].plot(sol.t * 1e3, i * 1e3, 'g-', linewidth=2, label='Current')

axes[1].set_xlabel('Time (ms)')

axes[1].set_ylabel('Current (mA)')

axes[1].set_title('RLC Circuit Current')

axes[1].legend()

axes[1].grid(True, alpha=0.3)

plt.tight_layout()

plt.savefig('rlc_response.png', dpi=150)

plt.show()

6.2 Analytical Solutions with SymPy

Symbol definitions

t, s = sp.symbols('t s')

R_sym, L_sym, C_sym = sp.symbols('R L C', positive=True)

Vs_sym = sp.Symbol('Vs')

RC circuit analytical solution (solving ODE with sympy)

v = sp.Function('v')

tau = sp.Symbol('tau', positive=True)

dv/dt = (Vs - v) / (R*C)

ode = sp.Eq(v(t).diff(t), (Vs_sym - v(t)) / (R_sym * C_sym))

Initial condition v(0) = 0

solution = sp.dsolve(ode, v(t), ics={v(0): 0})

print("RC circuit solution:")

print(solution)

sp.pprint(solution)

Laplace transform calculation

f = sp.exp(-t) * sp.sin(2*t)

F = sp.laplace_transform(f, t, s, noconds=True)

print(f"\nL{{e^(-t) sin(2t)}} = {F}")

Inverse Laplace transform

F_inv = sp.inverse_laplace_transform(

1 / (s**2 + 2*s + 5), s, t

)

print(f"L^(-1){{1/(s^2+2s+5)}} = {F_inv}")

6.3 Comparing Different Integrators

Euler method vs RK4 comparison

Example: y' = -y, y(0) = 1, exact solution y = e^(-t)

def f(t, y):

return -y

def euler_method(f, t0, y0, h, n):

t = np.zeros(n+1)

y = np.zeros(n+1)

t[0], y[0] = t0, y0

for i in range(n):

y[i+1] = y[i] + h * f(t[i], y[i])

t[i+1] = t[i] + h

return t, y

def rk4_method(f, t0, y0, h, n):

t = np.zeros(n+1)

y = np.zeros(n+1)

t[0], y[0] = t0, y0

for i in range(n):

k1 = h * f(t[i], y[i])

k2 = h * f(t[i] + h/2, y[i] + k1/2)

k3 = h * f(t[i] + h/2, y[i] + k2/2)

k4 = h * f(t[i] + h, y[i] + k3)

y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6

t[i+1] = t[i] + h

return t, y

h = 0.5

n = 10

t0, y0 = 0.0, 1.0

t_euler, y_euler = euler_method(f, t0, y0, h, n)

t_rk4, y_rk4 = rk4_method(f, t0, y0, h, n)

t_exact = np.linspace(0, n*h, 200)

y_exact = np.exp(-t_exact)

plt.figure(figsize=(10, 5))

plt.plot(t_exact, y_exact, 'k-', linewidth=2, label='Exact solution')

plt.plot(t_euler, y_euler, 'r--o', label=f'Euler method (h={h})')

plt.plot(t_rk4, y_rk4, 'b--s', label=f'RK4 (h={h})')

plt.xlabel('t')

plt.ylabel("y(t)")

plt.title("Euler vs RK4 Comparison: y' = -y")

plt.legend()

plt.grid(True, alpha=0.3)

plt.tight_layout()

plt.show()

Error comparison

print(f"Euler method error at t=5: {abs(y_euler[-1] - np.exp(-5)):.6f}")

print(f"RK4 error at t=5: {abs(y_rk4[-1] - np.exp(-5)):.6f}")

6.4 Automating Laplace Transforms with Python

s, t = sp.symbols('s t', real=True)

a, omega = sp.symbols('a omega', positive=True)

Auto-generate Laplace transform table

functions = {

'1': sp.Integer(1),

't': t,

't^2': t**2,

't^n': t**3,

'e^(at)': sp.exp(a*t),

'sin(wt)': sp.sin(omega*t),

'cos(wt)': sp.cos(omega*t),

't*e^(at)': t * sp.exp(a*t),

'e^(at)*sin(wt)': sp.exp(a*t) * sp.sin(omega*t),

}

print("Laplace Transform Table:")

print("-" * 50)

for name, func in functions.items():

F = sp.laplace_transform(func, t, s, noconds=True)

print(f"L{{{name}}} = {sp.simplify(F)}")

7. Comprehensive Engineering Application Examples

7.1 Transient Analysis of RLC Circuit with Power Source

**Problem**: A series RLC circuit with $R = 50\,\Omega$, $L = 0.2\,\text{H}$, $C = 100\,\mu\text{F}$ is connected to a $V_s = 100\,\text{V}$ DC source at $t=0$. Initial conditions are $v_C(0) = 0$, $i(0) = 0$.

**Solving with the Laplace Transform**:

Kirchhoff's voltage law:

$$L\frac{di}{dt} + Ri + \frac{1}{C}\int i\,dt = V_s u(t)$$

Substituting $i = C dv_C/dt$ and writing in terms of charge $q$:

$$L\ddot{q} + R\dot{q} + \frac{q}{C} = V_s$$

Laplace transform:

$$\left(Ls^2 + Rs + \frac{1}{C}\right)Q(s) = \frac{V_s}{s}$$

$$Q(s) = \frac{V_s}{s\left(Ls^2 + Rs + \frac{1}{C}\right)} = \frac{V_s/L}{s\left(s^2 + \frac{R}{L}s + \frac{1}{LC}\right)}$$

Substituting values: $\alpha = R/(2L) = 125\,\text{s}^{-1}$, $\omega_0 = 1/\sqrt{LC} \approx 223.6\,\text{s}^{-1}$

Since $\omega_0 > \alpha$, this is underdamped.

$$\omega_d = \sqrt{\omega_0^2 - \alpha^2} \approx 185.0\,\text{s}^{-1}$$

After partial fraction decomposition and inverse Laplace transform:

$$v_C(t) = V_s\left[1 - e^{-\alpha t}\left(\cos\omega_d t + \frac{\alpha}{\omega_d}\sin\omega_d t\right)\right]$$

$$= 100\left[1 - e^{-125t}\left(\cos 185t + 0.676\sin 185t\right)\right]\,\text{V}$$

7.2 Resonance Analysis of Spring-Mass-Damper System

**Problem**: $m = 1\,\text{kg}$, $k = 100\,\mathrm{N/m}$, $c = 4\,\mathrm{N \cdot s/m}$, $F(t) = 10\cos\omega t\,\mathrm{N}$

Natural frequency: $\omega_n = \sqrt{k/m} = 10\,\text{rad/s}$

Damping ratio: $\zeta = c/(2\sqrt{mk}) = 0.2$ (underdamped)

Steady-state response amplitude:

$$X = \frac{F_0/k}{\sqrt{(1-r^2)^2 + (2\zeta r)^2}}$$

At $r = \omega/\omega_n = 1$ (resonance):

$$X_{resonance} = \frac{F_0/k}{2\zeta} = \frac{10/100}{2 \times 0.2} = 0.25\,\text{m}$$

This is **2.5 times** the static deflection $F_0/k = 0.1\,\text{m}$ without resonance. This is the danger of resonance.

Summary and Next Steps

Topics covered in this article:

1. **Differential equation basics**: Classification, order, linearity, types of solutions

2. **First-order ODE methods**: Separation of variables, integrating factor, exact ODE, Bernoulli equation

3. **Second-order linear ODE**: Characteristic equation, three cases, complete RLC circuit analysis

4. **Non-homogeneous equations**: Undetermined coefficients, variation of parameters, forced vibration

5. **Systems of ODEs**: Matrix representation, eigenvalue/eigenvector solution

6. **Laplace transform**: Definition, transform pairs, properties, inverse transform, ODE solving

7. **Python implementation**: Numerical and symbolic solutions with scipy and sympy

The next article will cover **Fourier series/transform and partial differential equations (PDEs)** — the core mathematical tools for signal processing and electromagnetics.

References

- Kreyszig, E. "Advanced Engineering Mathematics", 10th Edition, Wiley

- Boyce, W. & DiPrima, R. "Elementary Differential Equations", 11th Edition, Wiley

- Simmons, G. "Differential Equations with Applications and Historical Notes", 3rd Edition

- [SciPy ODE Solver Documentation](https://docs.scipy.org/doc/scipy/reference/integrate.html)

- [SymPy Laplace Transform Documentation](https://docs.sympy.org/latest/modules/integrals/integrals.html)

Quiz

**Answer**: The integrating factor method is used to solve first-order linear ODEs of the standard form dy/dx + P(x)y = Q(x). The integrating factor is mu(x) = exp(integral of P(x) dx).

**Explanation**: Multiplying both sides by mu(x) transforms the left side into the derivative of mu(x)\*y, allowing direct integration. This is a systematic technique applicable to any first-order linear ODE.

**Answer**: The three cases depend on the damping coefficient alpha = R/(2L) and resonant frequency omega_0 = 1/sqrt(LC):

- Overdamped: alpha greater than omega_0, i.e., R greater than 2\*sqrt(L/C) — slow decay without oscillation

- Critically damped: alpha equals omega_0, i.e., R equals 2\*sqrt(L/C) — fastest decay without oscillation

- Underdamped: alpha less than omega_0, i.e., R less than 2\*sqrt(L/C) — oscillates and decays

**Explanation**: The characteristic equation r^2 + (R/L)r + 1/(LC) = 0 has roots r = -alpha plus/minus sqrt(alpha^2 - omega_0^2). The sign of the discriminant determines whether roots are real or complex.

**Answer**: The convolution theorem states that the Laplace transform of the convolution of two functions equals the product of their individual Laplace transforms: L\{(f*g)(t)\} = F(s)*G(s).

**Explanation**: This theorem is fundamental to engineering because it proves that the output of a linear time-invariant system equals the convolution of the input signal with the system's impulse response. In the s-domain, this becomes simple multiplication, making system analysis much easier.

**Answer**: Resonance occurs when the frequency ratio r = omega/omega_n equals 1, meaning the driving frequency equals the natural frequency of the system. Under resonance conditions, the steady-state amplitude is X = (F_0/k) / (2\*zeta), which can be many times larger than the static deflection.

**Explanation**: In the method of undetermined coefficients for y'' + 4y = 3\*cos(2x), the forcing frequency matches the natural frequency, causing the particular solution to include a factor of x (secular growth). In real systems with damping, amplitude does not grow to infinity but reaches a finite maximum at resonance. Avoiding resonance is critical in structural and mechanical engineering.

**Answer**: For a system dX/dt = AX, we assume solutions of the form X = v*exp(lambda*t). Substituting gives lambda*v = A*v, so lambda are eigenvalues and v are eigenvectors of matrix A. The general solution is a linear combination of all eigenvalue-eigenvector pairs: X = sum of C*k * v*k * exp(lambda_k \* t).

**Explanation**: Eigenvalues determine the stability and time constants of each mode. If all eigenvalues have negative real parts, the system is stable. The eigenvectors define the directions of each natural mode of the system.