Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section you will be able to:

Explain what the Laplace transform does to a function of $t$ and why it turns differential equations into algebraic ones.
Derive the rules $\mathcal{L}[y'] = sY - y(0)$ and $\mathcal{L}[y''] = s^2 Y - sy(0) - y'(0)$ straight from the definition.
Solve an initial-value problem for $ay'' + by' + cy = f(t)$ by applying the transform, solving an algebraic equation in $Y(s)$ , and inverting back to $y(t)$ .
Interpret the poles of $Y(s)$ in the s-plane as stability/oscillation indicators of the time-domain response.
Handle discontinuous forcings (step pulses, impulses) using the Heaviside and Dirac symbols — exactly the case where undetermined coefficients breaks down.
Verify a Laplace-derived solution numerically using SciPy/PyTorch's ODE solvers to machine precision.

The Big Picture: Why Transform At All?

"A differential equation is hard. An algebraic equation is easy. The Laplace transform turns the first into the second."

Up to now, solving $ay'' + by' + cy = f(t)$ has required two separate skills: find the homogeneous solution from the characteristic roots, then guess a particular solution that matches the forcing. The recipe works, but it breaks the moment $f(t)$ is a non-elementary expression — a rectangular pulse, an instantaneous impulse, a piecewise input from a controller. Every real engineering signal looks like that.

The Laplace transform offers a different deal. Pay one fixed cost up front — set up an integral that turns the time function $f(t)$ into a function $F(s)$ of a new variable $s$ — and in exchange, all derivatives become multiplications by $s$ , all integrals become divisions, all initial conditions get baked in automatically, and discontinuous forcings become tidy algebraic expressions.

The contract

The Laplace transform converts:

Calculus operations on $y(t)$ → algebra on $Y(s)$ .
Initial conditions → constants inside the same algebra (no separate step).
Piecewise / discontinuous forcing → smooth exponential factors $e^{-as}$ in the s-domain.

Control engineering

Transfer functions $H(s) = Y(s)/U(s)$ are the lingua franca of feedback control. They are literally Laplace transforms of impulse responses.

Signal processing

Filter design uses pole-zero plots in the s-plane. Every audio equaliser and image deblurring algorithm relies on this picture.

Electrical circuits

Replace $L$ , $R$ , $C$ with impedances $sL$ , $R$ , $1/(sC)$ and Kirchhoff's laws become ordinary Ohm-style algebra.

Defining the Laplace Transform

The Laplace transform of a function $f(t)$ defined for $t \ge 0$ is the function $F(s)$ given by

F(s) = \mathcal{L}\{f(t)\} = \int_{0}^{\infty} e^{-st}\,f(t)\,dt.

The variable $s$ is in general complex; the integral converges for $\operatorname{Re}(s)$ larger than some threshold determined by how fast $f$ grows. Read the formula slowly:

Multiply $f(t)$ by an exponentially decaying weight $e^{-st}$ . Large $s$ kills the tail, small $s$ probes it.
Integrate the product from $t = 0$ to $\infty$ . The output is a single number for each chosen $s$ .
Sweep $s$ over all admissible values; you now have a brand-new function $F(s)$ that encodes everything about $f(t)$ in a different variable.

A mental image

Think of $e^{-st}$ as a spotlight. Different $s$ values shine the light with different decay rates. The integral measures "how much of $f$ overlaps with this exponential spotlight." Big $s$ sees only the early part of $f$ ; small $s$ averages over the long tail.

Three transforms by hand

Three pairs you will use over and over. Each one is one integral.

1. The constant 1.

\mathcal{L}\{1\} = \int_0^{\infty} e^{-st}\,dt = \left[-\frac{e^{-st}}{s}\right]_0^{\infty} = \frac{1}{s}\quad (\operatorname{Re}(s) > 0).

2. The exponential $e^{at}$ .

\mathcal{L}\{e^{at}\} = \int_0^{\infty} e^{-(s-a)t}\,dt = \frac{1}{s - a}\quad (\operatorname{Re}(s) > a).

3. The sinusoid $\sin(\omega t)$ . Use

\sin(\omega t) = (e^{i\omega t} - e^{-i\omega t})/(2i)

and apply rule 2 twice:

\mathcal{L}\{\sin(\omega t)\} = \frac{1}{2i}\left(\frac{1}{s - i\omega} - \frac{1}{s + i\omega}\right) = \frac{\omega}{s^2 + \omega^2}.

The exponential is the eigenfunction

Look at pair 2: the Laplace transform of $e^{at}$ is a single shifted reciprocal $1/(s-a)$ . This is the deeper reason everything works — the exponential $e^{-st}$ inside the integral "sees" only the corresponding component of $f$ . Sinusoids and polynomials are just combinations of exponentials, so each gets a clean rational image in $s$ .

Linearity: The Foundational Property

Because the integral is linear, so is the Laplace transform:

\mathcal{L}\{\alpha\,f(t) + \beta\,g(t)\} = \alpha\,F(s) + \beta\,G(s).

This is the property that lets us transform an ODE term by term. Add and multiply by constants on the time side, and the same operations appear on the s side. Combined with the derivative rule below, every linear ODE with constant coefficients reduces to an algebraic problem.

The Magic Rule: Derivatives Become Multiplication by $s$

This is the rule that makes the whole machine work. We will derive it from the definition using a single integration by parts.

First derivative

Start from the definition with $f$ replaced by $f'$ :

\mathcal{L}\{f'(t)\} = \int_0^{\infty} e^{-st}\,f'(t)\,dt.

Integrate by parts with $u = e^{-st}$ , $dv = f'(t)\,dt$ . Then $du = -s\,e^{-st}\,dt$ and $v = f(t)$ . So

\int_0^{\infty} e^{-st} f'(t)\,dt = \bigl[e^{-st} f(t)\bigr]_0^{\infty} + s\int_0^{\infty} e^{-st} f(t)\,dt.

Assume $f$ does not grow faster than some exponential, so the boundary term at $t = \infty$ vanishes for $\operatorname{Re}(s)$ large enough. At $t = 0$ the boundary term is $-f(0)$ . We get the clean rule:

\boxed{\;\mathcal{L}\{f'(t)\} = s\,F(s) - f(0).\;}

Second derivative

Apply the first-derivative rule twice — that is, treat $f''$ as the derivative of $f'$ :

\mathcal{L}\{f''(t)\} = s\,\mathcal{L}\{f'(t)\} - f'(0) = s\bigl(s F(s) - f(0)\bigr) - f'(0).

Simplify:

\boxed{\;\mathcal{L}\{f''(t)\} = s^2 F(s) - s\,f(0) - f'(0).\;}

The pattern keeps going

By induction, $\mathcal{L}\{f^{(n)}\} = s^n F(s) - s^{n-1} f(0) - s^{n-2} f'(0) - \cdots - f^{(n-1)}(0)$ . Each differentiation in $t$ pulls out one power of $s$ and pays the price of one initial-condition term. The transform literally records how much initial information was needed to integrate the derivative away.

The whole trick in one breath

When you Laplace-transform $a y'' + b y' + c y = f(t)$ , the LHS becomes

a\bigl(s^2 Y - s y(0) - y'(0)\bigr) + b\bigl(s Y - y(0)\bigr) + c\,Y \;=\; (as^2 + bs + c)\,Y - \bigl(\text{IC terms}\bigr).

The coefficient of $Y$ on the s side is exactly the characteristic polynomial $as^2 + bs + c$ . The IC terms appear as a known constant on the right. Solve for $Y$ , invert. Done.

A Small Table You Will Use Forever

You only need a handful of pairs to handle every constant-coefficient 2nd-order ODE you will meet in physics and engineering. Memorise the first eight; derive the rest on the fly.

f(t) (time domain)	F(s) (s domain)	Comment
1	1/s	The DC level. Re(s) > 0.
t^n (n integer ≥ 0)	n! / s^(n+1)	Polynomial inputs.
e^(at)	1 / (s − a)	The exponential eigenfunction.
sin(ωt)	ω / (s² + ω²)	Pure oscillation.
cos(ωt)	s / (s² + ω²)	90° phase-shifted partner.
e^(at) sin(ωt)	ω / ((s − a)² + ω²)	Damped/growing sinusoid.
e^(at) cos(ωt)	(s − a) / ((s − a)² + ω²)	Same with cosine.
u(t − a)	e^(−a s) / s	Heaviside step at t = a.
δ(t − a)	e^(−a s)	Dirac impulse at t = a.
f'(t)	s F(s) − f(0)	Derivative rule — the engine of everything.
f''(t)	s² F(s) − s f(0) − f'(0)	Direct corollary; this is what we use for second-order ODEs.
∫₀ᵗ f(τ) dτ	F(s) / s	Integration ↔ division by s (the symmetric partner of the derivative rule).
u(t − a) f(t − a)	e^(−a s) F(s)	Second shifting theorem — vital for piecewise forcing.

Why the table is short

Every entry follows from $\mathcal{L}\{e^{at}\} = 1/(s-a)$ plus linearity. Sinusoids are sums of complex exponentials; polynomials come from differentiating $1/(s-a)$ with respect to $a$ ; damped sinusoids are products and shifts. If you ever forget an entry, derive it in 30 seconds with one of these tricks.

The Five-Step Recipe

Every constant-coefficient initial-value problem now follows the same five steps:

Transform. Apply $\mathcal{L}$ to both sides of the ODE. Replace $y(t)$ with $Y(s)$ , derivatives with their $s$ -rules.
Plug in ICs. Substitute $y(0)$ and $y'(0)$ . Now the equation is purely algebraic in $Y$ .
Solve. Isolate $Y(s) = (\text{stuff})/(\text{characteristic polynomial})$ .
Partial fractions. Break $Y(s)$ into pieces that match table entries.
Invert. Read each piece off the table backwards to get $y(t)$ .

The diagram in your head

Time domain (hard) → Laplace → s domain (easy algebra) → Inverse Laplace → Time domain (answer).

Worked Example (Step-by-Step)

Solve $y'' + 4 y = \sin(t)$ subject to $y(0) = 1,\; y'(0) = 0$ .

Click to expand the full hand calculation

Step 1 — Apply $\mathcal{L}$ to both sides.

\mathcal{L}\{y''\} + 4\,\mathcal{L}\{y\} = \mathcal{L}\{\sin t\}.

From the table: $\mathcal{L}\{y''\} = s^2 Y - s y(0) - y'(0)$ and $\mathcal{L}\{\sin t\} = 1/(s^2 + 1)$ .

Step 2 — Insert the ICs.

$y(0) = 1$ , $y'(0) = 0$ , so $\mathcal{L}\{y''\} = s^2 Y - s$ . The transformed equation is

(s^2 + 4)\,Y - s = \frac{1}{s^2 + 1}.

Step 3 — Solve for $Y(s)$ .

Y(s) = \frac{s}{s^2 + 4} + \frac{1}{(s^2 + 1)(s^2 + 4)}.

First term: response to the initial condition. Second term: response to the forcing. Notice that the algebra cleanly separates these two physical effects.

Step 4 — Partial fractions on the forcing piece.

Write $\frac{1}{(s^2 + 1)(s^2 + 4)} = \frac{A}{s^2 + 1} + \frac{B}{s^2 + 4}$ . Clear denominators: $1 = A(s^2 + 4) + B(s^2 + 1)$ . Set $s^2 = -1$ : $1 = 3A$ → $A = 1/3$ . Set $s^2 = -4$ : $1 = -3B$ → $B = -1/3$ .

Y(s) = \frac{s}{s^2 + 4} + \frac{1/3}{s^2 + 1} - \frac{1/3}{s^2 + 4}.

Step 5 — Inverse-Laplace term by term.

$\mathcal{L}^{-1}\{s/(s^2+4)\} = \cos(2t)$ . $\mathcal{L}^{-1}\{1/(s^2+1)\} = \sin(t)$ , scaled by $1/3$ . $\mathcal{L}^{-1}\{1/(s^2+4)\} = (1/2)\sin(2t)$ , scaled by $-1/3$ gives $-(1/6)\sin(2t)$ .

\boxed{\;y(t) = \cos(2t) + \tfrac{1}{3}\sin(t) - \tfrac{1}{6}\sin(2t).\;}

Step 6 — Sanity-check the ICs.

$y(0) = 1 + 0 - 0 = 1$ ✓. $y'(t) = -2\sin(2t) + \tfrac{1}{3}\cos(t) - \tfrac{1}{3}\cos(2t)$ gives $y'(0) = 0 + \tfrac{1}{3} - \tfrac{1}{3} = 0$ ✓.

Step 7 — Plug back into the ODE.

$y''(t) = -4\cos(2t) - \tfrac{1}{3}\sin(t) + \tfrac{2}{3}\sin(2t)$ . Then $y'' + 4 y = (-4 + 4)\cos(2t) + (-\tfrac{1}{3} + \tfrac{4}{3})\sin(t) + (\tfrac{2}{3} - \tfrac{2}{3})\sin(2t) = \sin(t)$ ✓. The Laplace machine produced exactly the solution.

Step 8 — Physical reading.

The system $y'' + 4y = f$ is an undamped oscillator with natural frequency $\omega_0 = 2$ . The forcing has frequency $\omega = 1 \neq \omega_0$ , so there is no resonance and the response stays bounded — exactly the mixture of cos(2t), sin(t), sin(2t) we found. If we had used $\sin(2t)$ as forcing instead, the inverse-Laplace would produce a $t\cos(2t)$ term that grows without bound: classical resonance, revealed by a repeated pole at $\pm 2i$ .

Interactive Laplace Explorer

Drive a generic damped second-order system $y'' + 2\zeta\omega_0\,y' + \omega_0^2\,y = f(t)$ with several canonical forcings and watch $Y(s)$ , its poles, and $y(t)$ update in real time. Click the forcing buttons to pick step, impulse, sinusoid, or free response (initial conditions only).

Loading Laplace explorer…

Three experiments worth doing

Resonance. Set forcing to $A\sin(\omega t)$ with $\zeta \approx 0.05$ , then drag $\omega$ toward $\omega_0$ . The two pole markers (✕ for the system, ○ for the forcing) get close to each other and the time-domain response amplitude spikes.
Stability boundary. Drop $\zeta$ to 0. The characteristic poles sit on the imaginary axis — perpetual oscillation.
Step settling. Pick the step input preset. Increase $\zeta$ and watch the response transition from ringing (underdamped) to a smooth s-curve (overdamped).

Drag the Poles: s-Plane ↔ t-Plane

The point of all this algebra is that the location of the poles of $Y(s)$ in the s-plane is the entire story of $y(t)$ . Real part < 0 means decay; real part > 0 means growth; imaginary part means oscillation; pole on the imaginary axis means undamped oscillation; two equal real poles mean critical damping with a $t\,e^{rt}$ signature. Move the poles, and the time response follows.

Loading s-plane dragger…

What you are really doing

The impulse response $h(t)$ is the inverse Laplace of $1/(s^2 + 2\zeta\omega_0 s + \omega_0^2)$ — a function with two poles. By dragging those poles you are literally manipulating $Y(s)$ and watching $y(t)$ . This is the same picture every control engineer uses to design feedback loops: place the closed-loop poles wherever you want the response to live.

Partial Fractions: How the Inverse Transform Really Works

After solving for $Y(s)$ , you almost always have a proper rational function $Y(s) = N(s)/D(s)$ . The inverse-Laplace algorithm is then: factor the denominator, write $Y$ as a sum of simpler rationals, and invert each piece off the table. Two flavours cover almost everything.

Distinct simple poles

If $D(s) = (s - r_1)(s - r_2)\cdots(s - r_n)$ with all $r_k$ different, then

Y(s) = \sum_{k=1}^{n} \frac{A_k}{s - r_k},\qquad A_k = \lim_{s \to r_k}(s - r_k)\,Y(s).

The residue formula on the right is the fastest way to compute $A_k$ by hand. Each term inverts to $A_k\,e^{r_k t}$ .

Complex conjugate pair

When the denominator contains an irreducible quadratic $(s - \alpha)^2 + \beta^2$ , do not factor into complex pieces — keep it as a real quadratic and write the numerator as $A(s - \alpha) + B\beta$ :

\frac{A(s - \alpha) + B\,\beta}{(s - \alpha)^2 + \beta^2} \;\longleftrightarrow\; e^{\alpha t}\bigl(A\cos\beta t + B\sin\beta t\bigr).

Repeated poles

If a pole $r$ has multiplicity $m$ , the partial-fraction expansion includes terms up to $1/(s - r)^m$ . The inverse transform of $1/(s - r)^k$ is $t^{k-1}\,e^{rt}/(k-1)!$ . That extra $t$ factor is exactly the same $t\,e^{rt}$ that appeared in section 22.03 — the Laplace transform makes its origin algebraically transparent.

Why partial fractions is unavoidable

The Laplace transform turns calculus into algebra, but it does NOT remove the need to factor a denominator. The hard work has just moved from "invent a particular solution" to "break a rational function into table entries." The good news: there is a deterministic mechanical procedure (partial fractions) — no guessing, ever.

Why Bother: Discontinuous Forcing

Now we earn the cost we paid up front. Suppose the input to our system is a rectangular pulse:

f(t) = \begin{cases} 1, & 1 \le t < 3 \\ 0, & \text{otherwise} \end{cases} \;=\; u(t-1) - u(t-3).

Where $u$ is the Heaviside step function: zero before its argument, one after. The method of undetermined coefficients chokes here — what would you guess for a particular solution that is on between $t = 1$ and $t = 3$ ? Laplace doesn't care:

\mathcal{L}\{u(t - a)\} = \frac{e^{-a s}}{s},\qquad \mathcal{L}\{u(t-a) g(t-a)\} = e^{-a s} G(s).

The first identity says "a switch turning on at $t = a$ appears as an exponential factor $e^{-as}$ in s." The second (the second shifting theorem) says "a time-shifted function gets multiplied by the same exponential." Together they are the entire vocabulary for piecewise inputs.

The Dirac impulse

The instantaneous "kick" $\delta(t - a)$ is the derivative of $u(t - a)$ . By the derivative rule (formally):

\mathcal{L}\{\delta(t - a)\} = e^{-a s}.

Impulses become just exponentials in s — no $1/s$ factor. Geometrically that makes sense: an impulse has zero width, so there is no integration happening — only the exponential timestamp survives.

Two physical scenarios you can finally solve

Hammer strike on a piano string. The string obeys a second-order ODE; the hammer is well-modelled as $A\,\delta(t - t_0)$ . Laplace gives you the post-strike vibration in one line.
Voltage switch onto an RLC circuit. The supply is $V\,u(t)$ ; the circuit equation is second-order. Laplace converts the discontinuity into a clean algebraic factor.

Symbolic Computation in Python (SymPy)

SymPy makes the entire pipeline literal: transform, substitute, solve, invert. The script below redoes the worked example end-to-end so you can confirm every intermediate quantity matches what you computed by hand.

Solve y'' + 4y = sin(t) symbolically with SymPy

🐍laplace_solve.py

Explanation(12)

Code(41)

1Import SymPy

SymPy is the Python computer-algebra system. It carries symbols exactly, so 'sin(t)' stays as the function sin(t) rather than a floating-point number — exactly what we need for Laplace algebra.

4Declare s and t

We mark s and t as positive real variables. The Laplace transform's integral runs over t ≥ 0, and the convergence region is Re(s) > a for some real a. Telling SymPy these facts up front gives much cleaner intermediate expressions.

5Declare the unknown function

y is a generic function of t. We have not committed to a formula for it — that is exactly what we are about to solve for.

8Encode the ODE

sp.Eq(LHS, RHS) is SymPy's way of writing 'LHS = RHS'. Notice that y(t).diff(t, 2) literally means d²y/dt² — the second derivative.

9Initial conditions as a dict

Two initial conditions because the ODE is second-order. We will plug these in by hand at the algebraic step — SymPy's symbolic Laplace transform leaves the IC terms abstract.

13Apply Laplace to the LHS

sp.laplace_transform is the workhorse. Setting noconds=True drops the convergence-domain metadata and returns just F(s). The LHS still contains a LaplaceTransform(y(t), t, s) placeholder because SymPy does not know what y is yet.

15Apply Laplace to the RHS

L[sin(t)] = 1/(s²+1). This is the single most-used pair in the entire table. SymPy returns it instantly.

20Introduce the symbol Y

We rename the abstract LaplaceTransform(y(t),t,s) to a plain symbol Y. From this point on the equation is ALGEBRAIC — no derivatives, no integrals. This is the whole reason the transform is worth doing.

26Bake the initial conditions into the LHS

Using L[y''] = s²Y − s y(0) − y'(0). With y(0)=1, y'(0)=0 the LHS becomes s²Y − s + 4Y. The IC information is now an additive constant on the LHS — it does not require a separate step later.

30Solve for Y(s)

Just one-variable algebra. SymPy returns Y as a rational function of s. After simplification you should see Y(s) = s/(s²+4) + 1/((s²+1)(s²+4)) — the first term is the response to the initial condition, the second is the response to the forcing.

35Inverse Laplace

sp.inverse_laplace_transform reads Y(s) and matches each term against its time-domain partner. For our example it returns y(t) = cos(2t) + (1/3) sin(t) − (1/6) sin(2t). The whole ODE — derivatives, ICs, forcing — has been crushed into one line of inverse transform.

39Sanity checks

Always verify the recovered solution against the original ICs. y(0) should equal 1 and y'(0) should equal 0. If either fails, the most common cause is a sign mistake when collecting the y(0), y'(0) terms inside L[y''].

29 lines without explanation

1import sympy as sp
2
3# Symbols.  s is the Laplace variable; t is time.
4t, s = sp.symbols("t s", positive=True, real=True)
5y = sp.Function("y")
6
7# The ODE: y'' + 4 y = sin(t),  y(0) = 1,  y'(0) = 0
8ode = sp.Eq(y(t).diff(t, 2) + 4 * y(t), sp.sin(t))
9ics = {y(0): 1, y(t).diff(t).subs(t, 0): 0}
10
11# Step 1.  Apply the Laplace transform to BOTH sides.
12# laplace_transform returns a 3-tuple (F(s), a, condition);
13# we want just F(s) so we use noconds=True.
14LHS = sp.laplace_transform(y(t).diff(t, 2) + 4 * y(t), t, s, noconds=True)
15RHS = sp.laplace_transform(sp.sin(t), t, s, noconds=True)
16print("L[y'' + 4y] =", LHS)        # contains LaplaceTransform(y(t), t, s)
17print("L[sin(t)]   =", RHS)        # 1/(s**2 + 1)
18
19# Step 2.  Substitute the symbol Y for L[y(t)], y(0)=1, y'(0)=0.
20Y = sp.symbols("Y")
21LHS_sub = LHS.rewrite(sp.LaplaceTransform).subs(
22    sp.LaplaceTransform(y(t), t, s), Y
23)
24# Manual substitution of the IC terms hidden inside L[y''(t)]:
25# L[y''] = s^2 Y - s y(0) - y'(0) = s^2 Y - s
26algebraic_eq = sp.Eq(s**2 * Y - s + 4 * Y, RHS)
27print("Algebraic equation:", algebraic_eq)
28
29# Step 3.  Solve for Y(s).
30Y_of_s = sp.solve(algebraic_eq, Y)[0]
31Y_of_s = sp.simplify(Y_of_s)
32print("Y(s) =", Y_of_s)
33
34# Step 4.  Inverse Laplace transform back to time domain.
35y_of_t = sp.inverse_laplace_transform(Y_of_s, s, t)
36y_of_t = sp.simplify(y_of_t)
37print("y(t) =", y_of_t)
38
39# Sanity checks: y(0) and y'(0).
40print("y(0)  =", sp.simplify(y_of_t.subs(t, 0)))
41print("y'(0) =", sp.simplify(sp.diff(y_of_t, t).subs(t, 0)))

Run the script. The output ends with $y(t) = \cos(2t) + \sin(t)/3 - \sin(2t)/6$ , exactly matching the by-hand answer. The $y(0)$ and $y'(0)$ checks return 1 and 0 — the initial conditions are encoded correctly through the algebra.

Numerical Cross-Check (SciPy / PyTorch)

Symbolic answers can lie if you typed something wrong. Numerical integration is a fully independent check: a different algorithm (Runge-Kutta) attacking the original ODE in the time domain, with no awareness of the Laplace machinery. Agreement to $10^{-9}$ or better is your assurance that the symbolic derivation is right.

Cross-check with SciPy's adaptive Runge-Kutta solver

🐍rk_verify.py

Explanation(6)

Code(27)

1Imports

NumPy for vectorized math, SciPy's solve_ivp for the adaptive Runge-Kutta integrator. No symbolic engine here — this is pure floating-point validation that the Laplace answer is correct.

4Convert to first-order system

Every classical numerical solver wants y'(t) = f(t, y), where y is a vector. So we set z₁=y, z₂=y'. The ODE y''=sin(t)−4y rewrites as the 2-D system (z₁', z₂') = (z₂, sin(t) − 4 z₁).

9Right-hand-side function

The body of rhs is one line of vector arithmetic. NumPy evaluates sin(t) and the linear combination element-wise — no Python loop is needed.

13Run the IVP solver

solve_ivp picks an adaptive step-size 5(4) Runge-Kutta method by default. dense_output=True lets us evaluate the spline approximation at any t in [0, 10] later. Tight tolerances rtol=1e-10, atol=1e-12 ensure that 'numerical error' is well below 'symbolic error' so the comparison is meaningful.

17Analytic formula

The closed form from the Laplace derivation. Vectorized over t thanks to NumPy.

21Compare at 5 samples

Both solvers should agree to about 10 digits. If they do, you have an independent confirmation that the Laplace algebra produced the correct answer. Printing the absolute error column lets you spot any systematic drift.

21 lines without explanation

1import numpy as np
2from scipy.integrate import solve_ivp
3
4# Re-cast the 2nd-order ODE as a 2-D first-order system:
5#   z1 = y,   z2 = y'
6#   z1' = z2
7#   z2' = sin(t) - 4 z1
8def rhs(t, z):
9    z1, z2 = z
10    return [z2, np.sin(t) - 4.0 * z1]
11
12# Solve from t=0 to t=10, dense output so we can sample anywhere.
13sol = solve_ivp(rhs, (0.0, 10.0), y0=[1.0, 0.0], dense_output=True,
14                rtol=1e-10, atol=1e-12)
15
16# Analytic solution from the Laplace derivation.
17def y_analytic(t):
18    return np.cos(2 * t) + (1.0 / 3.0) * np.sin(t) - (1.0 / 6.0) * np.sin(2 * t)
19
20# Compare at 5 sample times.
21ts = np.linspace(0, 10, 5)
22y_num = sol.sol(ts)[0]          # row 0 is y
23y_exact = y_analytic(ts)
24err = np.abs(y_num - y_exact)
25for ti, yn, ye, e in zip(ts, y_num, y_exact, err):
26    print(f"t = {ti:5.2f}   numeric = {yn:+.10f}   "
27          f"exact = {ye:+.10f}   |err| = {e:.2e}")

Expected output: at every sample time, the numeric and exact columns agree to ten or eleven decimal places. The same $solve\_ivp$ pattern works in PyTorch via $\texttt{torchdiffeq}$ ; the right-hand-side function becomes a callable returning a tensor, the integration is differentiable end-to-end, and you can place the entire ODE inside a neural network (a neural ODE) that learns the characteristic polynomial — which, by everything in this section, is literally the denominator of $Y(s)$ .

Bonus: Heaviside forcing with SymPy

The script below tackles the case that breaks undetermined coefficients: a rectangular pulse drives an undamped oscillator. SymPy returns the answer as a piecewise expression in two Heaviside functions — exactly the shape of the physical response.

Rectangular pulse forcing — the case Laplace was made for

🐍laplace_pulse.py

Explanation(6)

Code(25)

1SymPy import

We reuse the same symbolic stack. The new ingredient is Heaviside, which gives an exact symbolic representation of the unit-step function.

9Build the rectangular pulse

Subtracting two shifted step functions gives a rectangular pulse: turns on at t=1, turns off at t=3, height 1 in between, zero outside. This is the kind of forcing that breaks undetermined-coefficients but Laplace eats for breakfast.

13Laplace of the pulse

L[u(t − a)] = e^(−a s)/s, so L[u(t−1) − u(t−3)] = (e^(−s) − e^(−3 s))/s. The two exponentials in s are the algebraic signatures of the on/off times in t.

17Algebraic equation

Zero initial conditions kill the s y(0) and y'(0) terms, so L[y''] reduces to s²Y. The whole equation collapses to (s² + 4) Y = (e^(−s) − e^(−3s)) / s.

19Solve symbolically

Y(s) = (e^(−s) − e^(−3s)) / (s (s² + 4)). The two exponential factors are time-shift instructions that the inverse Laplace will interpret.

22Inverse Laplace with shifts

SymPy returns y(t) as a piecewise expression involving u(t−1) and u(t−3): nothing happens for t<1, the system starts ringing at t=1, then receives a second 'kick' at t=3 that adjusts the ringing. This is exactly the kind of answer that takes pages by hand and one line of inverse-Laplace.

19 lines without explanation

1import sympy as sp
2
3t, s = sp.symbols("t s", positive=True, real=True)
4y = sp.Function("y")
5
6# Forcing: a "pulse" — on between t=1 and t=3, zero elsewhere.
7#   f(t) = u(t - 1) - u(t - 3)
8# where u is the unit-step / Heaviside function.
9u = sp.Heaviside
10f = u(t - 1) - u(t - 3)
11
12# ODE:  y'' + 4 y = f(t),  y(0) = y'(0) = 0
13LHS = sp.laplace_transform(y(t).diff(t, 2) + 4 * y(t), t, s, noconds=True)
14RHS = sp.laplace_transform(f, t, s, noconds=True)
15print("L[f(t)]  =", sp.simplify(RHS))   # (e^{-s} - e^{-3s}) / s
16
17# With zero ICs:  L[y''] = s^2 Y.
18Y = sp.symbols("Y")
19algebraic = sp.Eq(s**2 * Y + 4 * Y, RHS)
20Y_sol = sp.solve(algebraic, Y)[0]
21print("Y(s) =", sp.simplify(Y_sol))
22
23# Inverse — note the exponentials e^{-as} become time shifts u(t - a) f(t - a).
24y_sol = sp.inverse_laplace_transform(Y_sol, s, t)
25print("y(t) =", sp.simplify(y_sol))

Common Pitfalls

Forgetting to plug in the ICs

The most common mistake: write $\mathcal{L}\{y''\} = s^2 Y$ and continue. That drops the $-s y(0) - y'(0)$ term. Your answer will satisfy the ODE but the wrong initial conditions.

Mis-shifting in the second shifting theorem

$\mathcal{L}\{f(t - a)\,u(t - a)\} = e^{-as}F(s)$ requires the function and the step to both be shifted by the same $a$ . A common slip is writing $\mathcal{L}\{f(t)\,u(t-a)\}$ with the unshifted $f$ ; that requires extra work and is NOT just $e^{-as}F(s)$ .

Partial fractions with complex poles

You can split a real quadratic $(s - \alpha)^2 + \beta^2$ into two complex linear factors, but it is almost always faster to keep it together and use the $e^{\alpha t}\cos\beta t$ / $e^{\alpha t}\sin\beta t$ table entries directly.

The integral does not always converge

For $f(t) = e^{t^2}$ there is no $s$ for which the defining integral converges. Functions that grow faster than any exponential do not have a Laplace transform. In practice this is rare in engineering — most physical signals are at worst exponentially bounded — but it is the correct mathematical caveat.

Summary

The Laplace transform is the single most powerful tool for solving linear ODEs with constant coefficients. Its full leverage comes from a small, structured set of rules:

Definition: $F(s) = \int_0^{\infty} e^{-st} f(t)\,dt$ .
Linearity: $\mathcal{L}\{\alpha f + \beta g\} = \alpha F + \beta G$ .
Derivative rule: $\mathcal{L}\{f''\} = s^2 F - s f(0) - f'(0)$ .
Shifting: $\mathcal{L}\{u(t-a) f(t-a)\} = e^{-a s} F(s)$ .
Inverse via partial fractions and a small table.

The s-plane picture

Pole location	Time-domain behaviour	Physical regime
Re(s) < 0, Im(s) = 0	Pure decay C e^(rt)	Overdamped — fast non-oscillating return
Re(s) = 0, Im(s) = 0 (repeated)	Linear in t times decay	Critically damped
Re(s) < 0, Im(s) ≠ 0	Damped sinusoid	Underdamped
Re(s) = 0, Im(s) ≠ 0	Pure sinusoid	Undamped oscillator
Re(s) > 0	Exponential blow-up	Unstable

The slogan

"Calculus in t is algebra in s. The poles of Y(s) are the entire fate of y(t)."

Coming up next: Chapter 23 opens the door to systems of differential equations — coupled equations where Laplace becomes a matrix identity and the poles of $Y(s)$ become the eigenvalues of the system matrix. Everything you learned here generalises one dimension at a time.