Variation of Parameters for Second-Order Linear ODEs
import numpy as np
import sympy as sym
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from IPython.display import Math, display
mpl.rcParams['figure.dpi'] = 150
mpl.rcParams['axes.spines.top'] = False
mpl.rcParams['axes.spines.right'] = FalseIn this section we will:
State the general theory for second-order linear ODEs with variable coefficients and explain when the method of undetermined coefficients fails.
Develop the method of variation of parameters as a general technique for finding particular solutions — applicable to any second-order linear ODE, including those with variable coefficients.
Work through a suite of detailed examples, including equations with trigonometric, exponential, and rational forcing functions, as well as an Euler equation with variable-coefficient homogeneous part.
Use SymPy and Python to verify solutions and visualize particular solutions alongside their homogeneous counterparts.
This material corresponds to Section 2.4.2 of (Logan 2015).
In Notes 4 we introduced variation of parameters as one of two methods for finding particular solutions of constant-coefficient equations, alongside undetermined coefficients. Here we develop the method more fully: we derive the formula from scratch, prove when it applies, work through four detailed examples (including the Euler equation, which has variable coefficients), and clarify exactly when to prefer it over undetermined coefficients. If you used variation of parameters in Notes 4, everything here will be familiar in structure — the new content is the generality and the examples.
A second-order linear ODE with variable coefficients has the form \[ p(t)\,x'' + q(t)\,x' + r(t)\,x = f(t), \tag{VC} \] where the coefficient functions \(p(t)\), \(q(t)\), \(r(t)\) and the forcing function \(f(t)\) are continuous on some interval \(I\), with \(p(t) \neq 0\) on \(I\). Dividing by \(p(t)\) puts the equation in standard form: \[ x'' + P(t)\,x' + Q(t)\,x = g(t), \tag{SF} \] where \(P = q/p\), \(Q = r/p\), and \(g = f/p\).
The Superposition Principle and the general solution structure \(x = x_h + x_p\) hold for (SF) exactly as for the constant-coefficient case. What changes is the method we use to find \(x_h\) and \(x_p\), since the exponential ansatz \(e^{\lambda t}\) no longer works when the coefficients vary.
The method of undetermined coefficients requires the right-hand side \(g(t)\) to be a polynomial, exponential, sine, cosine, or product thereof. For variable-coefficient equations the homogeneous solution \(x_h\) no longer consists of these simple forms, so there is no longer a predictable guess for \(x_p\). Moreover, even for constant-coefficient equations, forcing functions such as \(\sec t\), \(\tan t\), \(\ln t\), \(1/t\), or \(e^t/t\) lie outside the undetermined-coefficients table and cannot be handled by guessing.
Variation of parameters (developed by Lagrange) works for any continuous \(f(t)\).
| Forcing \(g(t)\) | Why UC fails |
|---|---|
| \(\sec t\), \(\tan t\), \(\csc t\) | Not in the guess table; derivatives generate new function families |
| \(\ln t\), \(\sqrt{t}\), \(1/t\) | Not closed under differentiation — infinitely many new terms appear |
| \(e^t/t\) | Product of an exponential and \(1/t\); \(1/t\) grows under differentiation |
| Variable coefficients | No exponential \(x_h\) to compare against; guessing is not predictable |
Variation of parameters handles all of these by working directly with the homogeneous solutions, whatever form they take.
Suppose we know two linearly independent solutions \(x_1(t)\) and \(x_2(t)\) of the homogeneous equation \[ x'' + P(t)\,x' + Q(t)\,x = 0. \tag{H} \] Their general solution is \(x_h = C_1 x_1 + C_2 x_2\). In variation of parameters, we allow the constants to vary: seek a particular solution of the form \[ x_p(t) = u_1(t)\,x_1(t) + u_2(t)\,x_2(t) \] where \(u_1(t)\) and \(u_2(t)\) are functions to be determined.
Computing \(x_p'\): \[ x_p' = u_1'x_1 + u_1 x_1' + u_2'x_2 + u_2 x_2'. \] We impose the side condition \[ u_1'x_1 + u_2'x_2 = 0 \tag{SC} \] to eliminate second derivatives of \(u_1\), \(u_2\) when we differentiate again. This gives \(x_p' = u_1 x_1' + u_2 x_2'\) and \[ x_p'' = u_1'x_1' + u_1 x_1'' + u_2'x_2' + u_2 x_2''. \] Substituting into (SF) and using the fact that \(x_1\), \(x_2\) satisfy (H): \[ u_1'x_1' + u_2'x_2' = g(t). \tag{D} \] Equations (SC) and (D) form the \(2\times 2\) linear system for \(u_1'\) and \(u_2'\): \[ \begin{pmatrix} x_1 & x_2 \\ x_1' & x_2' \end{pmatrix} \begin{pmatrix} u_1' \\ u_2' \end{pmatrix} = \begin{pmatrix} 0 \\ g(t) \end{pmatrix}. \] The coefficient matrix is the Wronskian matrix whose determinant \[ W(x_1, x_2) = x_1 x_2' - x_2 x_1' \neq 0 \] is non-zero since \(x_1\), \(x_2\) are linearly independent. Solving by Cramer’s rule: \[ u_1' = -\frac{x_2\,g}{W}, \qquad u_2' = \frac{x_1\,g}{W}. \] Integrating: \[ u_1 = -\int\frac{x_2\,g}{W}\,dt, \qquad u_2 = \int\frac{x_1\,g}{W}\,dt. \]
Given two linearly independent solutions \(x_1\), \(x_2\) of (H) with Wronskian \(W = x_1 x_2' - x_2 x_1'\), a particular solution of \(x'' + Px' + Qx = g\) is \[ \boxed{x_p(t) = -x_1(t)\int\frac{x_2(t)\,g(t)}{W(t)}\,dt \;+\; x_2(t)\int\frac{x_1(t)\,g(t)}{W(t)}\,dt.} \] The general solution is \(x(t) = C_1 x_1(t) + C_2 x_2(t) + x_p(t)\).
The Wronskian of any two solutions of (H) satisfies Abel’s identity: \[ W(t) = W(t_0)\,\exp\!\left(-\int_{t_0}^t P(s)\,ds\right). \] For constant-coefficient equations (\(P\) constant) this gives \(W(t) = W(t_0)\,e^{-P(t-t_0)}\), and in particular \(W \neq 0\) everywhere if it is non-zero at one point. More generally, since the exponential factor is always nonzero, Abel’s identity shows that \(W\) either vanishes identically on \(I\) or is nonzero everywhere on \(I\) — in particular, if \(W(t_0)\neq 0\) for any single \(t_0\in I\), then \(W(t)\neq 0\) for all \(t\in I\). This justifies checking the Wronskian at just one convenient point.
The forcing function \(\sec t\) is not in the undetermined-coefficients table — we must use variation of parameters.
Step 1. Already in standard form; \(g(t) = \sec t\).
Step 2. Homogeneous equation \(x'' + x = 0\) has eigenvalues \(\lambda = \pm i\), giving \[x_1 = \cos t, \qquad x_2 = \sin t.\]
Step 3. Wronskian: \[W = \cos t\cdot\cos t - \sin t\cdot(-\sin t) = \cos^2 t + \sin^2 t = 1.\]
Step 4. \[u_1' = -\frac{\sin t\cdot\sec t}{1} = -\tan t, \qquad u_2' = \frac{\cos t\cdot\sec t}{1} = 1.\]
Step 5. \[u_1 = \ln|\cos t|, \qquad u_2 = t.\]
Step 6. \[x_p = \cos t\cdot\ln|\cos t| + t\sin t.\]
General solution: \[\boxed{x(t) = C_1\cos t + C_2\sin t + \cos t\ln|\cos t| + t\sin t.}\]
The solution exists on the interval \((-\pi/2,\,\pi/2)\) containing \(t=0\), and more generally on any open interval between consecutive singularities of \(\sec t\) (i.e., between adjacent odd multiples of \(\pi/2\)). The \(t\sin t\) term in \(x_p\) causes the amplitude to grow linearly — behavior analogous to resonance, even though \(\omega = 1\) and the forcing \(\sec t\) is not a pure sinusoid. Undetermined coefficients never produces this kind of growth.
t_plot = np.linspace(0, 1.4, 400)
# Particular solution (no constants)
xp1 = np.cos(t_plot)*np.log(np.cos(t_plot)) + t_plot*np.sin(t_plot)
# IVP: xp(0)=0, xp'(0)=0, so x(0)=C1=1, x'(0)=C2=0
xh1 = np.cos(t_plot) # C1=1, C2=0
x1_full = xh1 + xp1
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(t_plot, x1_full, color='steelblue', lw=2.5, label='Full solution $x(t)$')
ax.plot(t_plot, xh1, color='darkorange', lw=1.8, ls='--', label='Homogeneous $x_h=\\cos t$')
ax.plot(t_plot, xp1, color='seagreen', lw=1.8, ls='-.',
label='Particular $x_p=\\cos t\\ln|\\cos t|+t\\sin t$')
# Numerical check
def ode1(t, y):
return [y[1], -y[0] + 1/np.cos(t)]
sol1 = solve_ivp(ode1, (0, 1.4), [1.0, 0.0], dense_output=True, max_step=0.005)
t_dots = np.linspace(0, 1.4, 18)
ax.plot(t_dots, sol1.sol(t_dots)[0], 'ro', markersize=5, label='Numerical check')
ax.axvline(np.pi/2, color='gray', ls=':', lw=1.2, label=r'Singularity $t=\pi/2$')
ax.axhline(0, color='k', lw=0.5)
ax.set_xlabel('$t$'); ax.set_ylabel('$x(t)$')
ax.set_title(r"$x''+x=\sec t$, $x(0)=1$, $x'(0)=0$")
ax.legend(fontsize=8.5); ax.set_ylim(-0.5, 2.5)
plt.tight_layout(); plt.show()
This equation has a repeated eigenvalue — the homogeneous solutions are \(x_1 = e^t\) and \(x_2 = te^t\), and the forcing \(e^t/t\) is outside the undetermined-coefficients table (the \(1/t\) factor makes it inadmissible).
Step 1. Already in standard form; \(g(t) = e^t/t\).
Step 2. Characteristic equation \((\lambda-1)^2 = 0\) gives \(\lambda = 1\) (repeated): \[x_1 = e^t, \qquad x_2 = te^t.\]
Step 3. Wronskian: \[W = e^t(e^t + te^t) - te^t\cdot e^t = e^{2t} + te^{2t} - te^{2t} = e^{2t}.\]
Step 4. \[u_1' = -\frac{te^t\cdot(e^t/t)}{e^{2t}} = -\frac{e^{2t}}{e^{2t}} = -1, \qquad u_2' = \frac{e^t\cdot(e^t/t)}{e^{2t}} = \frac{1}{t}.\]
Step 5. \[u_1 = -t, \qquad u_2 = \ln t \quad (t > 0).\]
Step 6. \[x_p = e^t\cdot(-t) + te^t\cdot\ln t = te^t(\ln t - 1).\]
General solution: \[\boxed{x(t) = C_1 e^t + C_2 te^t + te^t(\ln t - 1), \quad t > 0.}\]
Note that \(te^t(\ln t - 1) = te^t\ln t - te^t\), and the \(-te^t\) part could be absorbed into the \(C_2 te^t\) term by redefining \(C_2\). This is valid because particular solutions are not unique — any particular solution differs from any other by a solution of the homogeneous equation, so absorbing constant-multiple terms from \(x_p\) into \(x_h\) yields the same general solution. In practice we keep the particular solution as found and let the constants adjust when ICs are applied.
# SymPy verification
t_sym = sym.Symbol('t', positive=True)
x_sym = sym.Function('x')
ode2_sym = sym.Eq(x_sym(t_sym).diff(t_sym,2) - 2*x_sym(t_sym).diff(t_sym) + x_sym(t_sym),
sym.exp(t_sym)/t_sym)
sol2_sym = sym.dsolve(ode2_sym, x_sym(t_sym))
print("SymPy solution:")
display(Math(sym.latex(sol2_sym)))SymPy solution:\(\displaystyle x{\left(t \right)} = \left(C_{1} + t \left(C_{2} + \log{\left(t \right)}\right)\right) e^{t}\)
t_plot2 = np.linspace(1.0, 2.8, 400)
# IVP: x(1)=0, x'(1)=0
# Solving: C1 + C2 - 1 = 0 and C1 + 2*C2 + 1 = 0 => C1=2, C2=-1
C1_2, C2_2 = 2.0, -1.0
xp2 = t_plot2 * np.exp(t_plot2) * (np.log(t_plot2) - 1)
xh2 = C1_2*np.exp(t_plot2) + C2_2*t_plot2*np.exp(t_plot2) # 2e^t - te^t
x2_full = xh2 + xp2
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(t_plot2, x2_full, color='steelblue', lw=2.5, label='Full solution $x(t)$')
ax.plot(t_plot2, xh2, color='darkorange', lw=1.8, ls='--', label='$x_h = 2e^t - te^t$ (with $C_1=2$, $C_2=-1$)')
ax.plot(t_plot2, xp2, color='seagreen', lw=1.8, ls='-.', label='$x_p = te^t(\\ln t-1)$')
def ode2(t, y): return [y[1], 2*y[1] - y[0] + np.exp(t)/t]
sol2n = solve_ivp(ode2, (1.0, 2.8), [0.0, 0.0], dense_output=True, max_step=0.005)
t_dots2 = np.linspace(1.0, 2.8, 18)
ax.plot(t_dots2, sol2n.sol(t_dots2)[0], 'ro', markersize=5, label='Numerical check')
ax.axhline(0, color='k', lw=0.5)
ax.set_xlabel('$t$'); ax.set_ylabel('$x(t)$')
ax.set_title(r"$x''-2x'+x=e^t/t$, $x(1)=0$, $x'(1)=0$")
ax.legend(fontsize=8.5)
plt.tight_layout(); plt.show()
This is a second-order linear ODE with variable coefficients — the coefficients \(-2t\) and \(2\) are functions of \(t\), not constants. The homogeneous part is an Euler–Cauchy equation, whose solutions are powers of \(t\).
Step 1. Divide by \(t^2\) (standard form, \(t \neq 0\)): \[x'' - \frac{2}{t}\,x' + \frac{2}{t^2}\,x = t\sin t. \qquad g(t) = t\sin t.\]
Step 2. For the homogeneous Euler equation \(t^2 x'' - 2tx' + 2x = 0\), try \(x = t^r\). The motivation is that the Euler equation has a scaling symmetry: replacing \(t\) by \(ct\) leaves the equation’s structure unchanged, and \(t^r\) is the family of functions preserved under this scaling — just as \(e^{\lambda t}\) is preserved under differentiation for constant-coefficient equations. Substituting: \[r(r-1)t^r - 2rt^r + 2t^r = t^r[r^2 - 3r + 2] = t^r(r-1)(r-2) = 0.\] So \(r = 1\) and \(r = 2\), giving: \[x_1 = t, \qquad x_2 = t^2.\]
Step 3. Wronskian: \[W = t\cdot 2t - t^2\cdot 1 = 2t^2 - t^2 = t^2.\]
Step 4. \[u_1' = -\frac{t^2\cdot t\sin t}{t^2} = -t\sin t, \qquad u_2' = \frac{t\cdot t\sin t}{t^2} = \sin t.\]
Step 5. Using integration by parts for \(u_1\): \[u_1 = -\int t\sin t\,dt = t\cos t - \sin t, \qquad u_2 = \int\sin t\,dt = -\cos t.\]
Step 6. \[x_p = t(t\cos t - \sin t) + t^2(-\cos t) = t^2\cos t - t\sin t - t^2\cos t = -t\sin t.\]
General solution: \[\boxed{x(t) = C_1 t + C_2 t^2 - t\sin t.}\]
# SymPy verification
t_sym = sym.Symbol('t', positive=True)
x_sym = sym.Function('x')
ode3_sym = sym.Eq(t_sym**2*x_sym(t_sym).diff(t_sym,2)
- 2*t_sym*x_sym(t_sym).diff(t_sym)
+ 2*x_sym(t_sym),
t_sym**3*sym.sin(t_sym))
sol3_sym = sym.dsolve(ode3_sym, x_sym(t_sym))
print("SymPy solution:")
display(Math(sym.latex(sol3_sym)))SymPy solution:\(\displaystyle x{\left(t \right)} = t \left(C_{1} + C_{2} t - \sin{\left(t \right)}\right)\)
t_plot3 = np.linspace(np.pi, 3*np.pi, 500)
# IVP: x(pi) = pi^2, x'(pi) = 2*pi
# x = C1*t + C2*t^2 - t*sin(t)
# x(pi): C1*pi + C2*pi^2 - pi*0 = pi*(C1 + C2*pi) = pi^2 => C1 + C2*pi = pi => C1 = pi - C2*pi
# x'(t) = C1 + 2*C2*t - sin(t) - t*cos(t)
# x'(pi): C1 + 2*C2*pi - 0 - pi*(-1) = C1 + 2*C2*pi + pi = 2*pi
# => C1 + 2*C2*pi = pi; with C1 = pi - C2*pi: pi - C2*pi + 2*C2*pi = pi + C2*pi = pi => C2=0, C1=pi
xp3 = -t_plot3 * np.sin(t_plot3)
xh3 = np.pi * t_plot3 # C1=pi, C2=0
x3_full = xh3 + xp3
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(t_plot3, x3_full, color='steelblue', lw=2.5, label='Full solution $x(t)$')
ax.plot(t_plot3, xh3, color='darkorange', lw=1.8, ls='--', label='$x_h = \\pi t$ (with $C_1=\\pi$, $C_2=0$)')
ax.plot(t_plot3, xp3, color='seagreen', lw=1.8, ls='-.', label='$x_p = -t\\sin t$')
def ode3(t, y):
return [y[1], (2*t*y[1] - 2*y[0] + t**3*np.sin(t))/t**2]
x0_3 = np.pi**2
xp0_3 = 2*np.pi
sol3n = solve_ivp(ode3, (np.pi, 3*np.pi), [x0_3, xp0_3], dense_output=True, max_step=0.01)
t_dots3 = np.linspace(np.pi, 3*np.pi, 20)
ax.plot(t_dots3, sol3n.sol(t_dots3)[0], 'ro', markersize=5, label='Numerical check')
ax.axhline(0, color='k', lw=0.5)
ax.set_xlabel('$t$'); ax.set_ylabel('$x(t)$')
ax.set_title(r"Euler eq: $t^2x''-2tx'+2x=t^3\sin t$, $x(\pi)=\pi^2$, $x'(\pi)=2\pi$")
ax.legend(fontsize=8.5)
plt.tight_layout(); plt.show()
Another trigonometric forcing function outside the undetermined-coefficients table.
Step 1. Already in standard form; \(g(t) = \csc t\).
Step 2. Same homogeneous equation as Example 1: \[x_1 = \cos t, \qquad x_2 = \sin t, \qquad W = 1.\]
Step 3. \[u_1' = -\frac{\sin t\cdot\csc t}{1} = -1, \qquad u_2' = \frac{\cos t\cdot\csc t}{1} = \cot t = \frac{\cos t}{\sin t}.\]
Step 4. \[u_1 = -t, \qquad u_2 = \int\cot t\,dt = \ln|\sin t|.\]
Step 5. \[x_p = \cos t\cdot(-t) + \sin t\cdot\ln|\sin t| = -t\cos t + \sin t\ln|\sin t|.\]
General solution: \[\boxed{x(t) = C_1\cos t + C_2\sin t - t\cos t + \sin t\ln|\sin t|.}\]
Compare Examples 1 and 4: both use \(x_1 = \cos t\), \(x_2 = \sin t\), and \(W = 1\). The only difference is the forcing — \(\sec t\) versus \(\csc t\) — yet the particular solutions are qualitatively different (\(\cos t\ln|\cos t| + t\sin t\) vs. \(-t\cos t + \sin t\ln|\sin t|\)). The method of variation of parameters is identical in both cases; only the integrals change. This illustrates how the method adapts automatically to the forcing, without any guessing or modification rule.
t_plot4 = np.linspace(0.1, np.pi - 0.05, 400) # between singularities
xp4 = -t_plot4*np.cos(t_plot4) + np.sin(t_plot4)*np.log(np.sin(t_plot4))
# IVP: x(pi/2)=pi/2, x'(pi/2)=0
# x = C1*cos(t) + C2*sin(t) + xp
# xp(pi/2) = -pi/2*0 + 1*ln(1) = 0
# xp'(t) = -cos(t) + t*sin(t) + cos(t)*ln(sin(t)) + sin(t)*cos(t)/sin(t)
# = -cos(t) + t*sin(t) + cos(t)*ln(sin(t)) + cos(t)
# = t*sin(t) + cos(t)*ln(sin(t))
# xp'(pi/2) = pi/2*1 + 0*ln(1) = pi/2
# x(pi/2): C1*0 + C2*1 + 0 = pi/2 => C2 = pi/2
# x'(pi/2): -C1*0 + C2*0 + pi/2 = ... wait
# x'(t) = -C1*sin(t) + C2*cos(t) + xp'(t)
# x'(pi/2) = -C1*1 + C2*0 + pi/2 = 0 => C1 = pi/2
C1_4, C2_4 = np.pi/2, np.pi/2
xh4 = C1_4*np.cos(t_plot4) + C2_4*np.sin(t_plot4)
x4_full = xh4 + xp4
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(t_plot4, x4_full, color='steelblue', lw=2.5, label='Full solution $x(t)$')
ax.plot(t_plot4, xh4, color='darkorange', lw=1.8, ls='--',
label=f'$x_h = \\frac{{\\pi}}{{2}}\\cos t + \\frac{{\\pi}}{{2}}\\sin t$')
ax.plot(t_plot4, xp4, color='seagreen', lw=1.8, ls='-.', label='$x_p=-t\\cos t+\\sin t\\ln|\\sin t|$')
def ode4(t, y):
return [y[1], -y[0] + 1/np.sin(t)]
sol4n = solve_ivp(ode4, (np.pi/2, np.pi-0.05), [np.pi/2, 0.0], dense_output=True, max_step=0.005)
t_dots4 = np.linspace(np.pi/2, np.pi-0.08, 10)
ax.plot(t_dots4, sol4n.sol(t_dots4)[0], 'ro', markersize=5, label='Numerical check')
ax.axhline(0, color='k', lw=0.5)
ax.axvline(np.pi, color='gray', ls=':', lw=1.2, label=r'Singularity at $t=\pi$')
ax.set_xlabel('$t$'); ax.set_ylabel('$x(t)$')
ax.set_title(r"$x''+x=\csc t$, $x(\pi/2)=\pi/2$, $x'(\pi/2)=0$, on $(0,\pi)$")
ax.legend(fontsize=8); ax.set_ylim(-1.5, 3)
plt.tight_layout(); plt.show()
t_sym = sym.Symbol('t', positive=True)
x_sym = sym.Function('x')
examples = {
r"x''+x=\sec t":
sym.Eq(x_sym(t_sym).diff(t_sym,2) + x_sym(t_sym), sym.sec(t_sym)),
r"x''-2x'+x=e^t/t":
sym.Eq(x_sym(t_sym).diff(t_sym,2) - 2*x_sym(t_sym).diff(t_sym) + x_sym(t_sym),
sym.exp(t_sym)/t_sym),
r"t^2x''-2tx'+2x=t^3\sin t":
sym.Eq(t_sym**2*x_sym(t_sym).diff(t_sym,2) - 2*t_sym*x_sym(t_sym).diff(t_sym)
+ 2*x_sym(t_sym), t_sym**3*sym.sin(t_sym)),
r"x''+x=\csc t":
sym.Eq(x_sym(t_sym).diff(t_sym,2) + x_sym(t_sym), sym.csc(t_sym)),
}
for label, ode in examples.items():
print(f"${label}$")
sol = sym.dsolve(ode, x_sym(t_sym))
display(Math(r"x(t)=" + sym.latex(sol.rhs)))
print()$x''+x=\sec t$\(\displaystyle x(t)=\left(C_{1} + t\right) \sin{\left(t \right)} + \left(C_{2} + \log{\left(\cos{\left(t \right)} \right)}\right) \cos{\left(t \right)}\)
$x''-2x'+x=e^t/t$\(\displaystyle x(t)=\left(C_{1} + t \left(C_{2} + \log{\left(t \right)}\right)\right) e^{t}\)
$t^2x''-2tx'+2x=t^3\sin t$\(\displaystyle x(t)=t \left(C_{1} + C_{2} t - \sin{\left(t \right)}\right)\)
$x''+x=\csc t$\(\displaystyle x(t)=\left(C_{1} - t\right) \cos{\left(t \right)} + \left(C_{2} + \log{\left(\sin{\left(t \right)} \right)}\right) \sin{\left(t \right)}\)
| Example | Homog. solutions | \(W\) | Forcing | Key feature |
|---|---|---|---|---|
| \(x''+x=\sec t\) | \(\cos t\), \(\sin t\) | \(1\) | Outside UC table | \(x_p\) grows like \(t\sin t\) |
| \(x''-2x'+x=e^t/t\) | \(e^t\), \(te^t\) | \(e^{2t}\) | Repeated eigenvalue + \(1/t\) | \(x_p\) grows like \(te^t\ln t\) |
| \(t^2x''-2tx'+2x=t^3\sin t\) | \(t\), \(t^2\) | \(t^2\) | Variable coefficients | Power-law homog. solutions |
| \(x''+x=\csc t\) | \(\cos t\), \(\sin t\) | \(1\) | Outside UC table | Log term in \(x_p\) |
| Method | Use when |
|---|---|
| Undetermined coefficients | \(f(t)\) is a polynomial, exponential, \(\sin\), \(\cos\), or their products; and the ODE has constant coefficients |
| Variation of parameters | \(f(t)\) is any continuous function (\(\sec t\), \(\tan t\), \(\ln t\), \(1/t\), etc.) or the ODE has variable coefficients |
Always try undetermined coefficients first when applicable — it is faster and avoids potentially difficult integrals. Switch to variation of parameters when the forcing falls outside the UC table, or when the coefficients of the ODE depend on \(t\).
Variation of parameters applies to any second-order linear ODE in standard form — constant or variable coefficients, any continuous \(g(t)\).
The side condition \(u_1'x_1+u_2'x_2=0\) is what makes the system solvable and keeps the algebra manageable.
The Wronskian \(W=x_1 x_2'-x_2 x_1'\) is the key quantity — it must be non-zero on the interval (guaranteed by linear independence of \(x_1\), \(x_2\)).
Forcing functions outside the UC table (\(\sec t\), \(\csc t\), \(1/t\), \(\ln t\), etc.) lead to non-elementary integrals — but SymPy handles them automatically.
For Euler equations, the homogeneous solutions are powers \(t^r\) rather than exponentials; variation of parameters proceeds identically once \(x_1\), \(x_2\) are found.
These notes are also viewable as a slide deck presentation: Open slides in full screen
With the method of variation of parameters, we now have a complete toolkit for second-order linear ODEs — constant or variable coefficients, any continuous forcing. The next major step is to reframe everything in matrix language: a second-order ODE \(x'' + px' + qx = 0\) is equivalent to the first-order system \(\mathbf{y}' = A\mathbf{y}\), and the homogeneous solutions \(x_1\), \(x_2\) become columns of the fundamental matrix. The Wronskian becomes the determinant of that matrix, Abel’s identity becomes Liouville’s theorem, and variation of parameters extends naturally to systems of any size. This unification is the subject of Chapter 3 of Logan (2015).
Next: Introduction to Laplace transforms — Logan §3.1.
import sys
print("Python version:", sys.version)
print('\n'.join(f'{m.__name__}=={m.__version__}' for m in globals().values() if getattr(m, '__version__', None)))Python version: 3.14.4 | packaged by conda-forge | (main, Apr 8 2026, 02:33:53) [Clang 20.1.8 ]
numpy==2.4.3
sympy==1.14.0
matplotlib==3.10.8