Variation of Parameters

Once you can solve the free problem $L[y]=0$ , variation of parameters hands you the forced problem $L[y]=f$ almost for nothing. The trick is disarmingly simple and, in hindsight, inevitable. The homogeneous general solution is a linear combination $c_1y_1+\cdots+c_ny_n$ of a fundamental system with constant coefficients. Freeze those constants and the trajectory can only trace free motions; let them vary — promote $c_i$ to functions $u_i(x)$ — and the same combination acquires exactly enough freedom to bend, at every instant, into whatever shape an arbitrary forcing $f(x)$ demands. Lagrange saw that the extra freedom is more than you need: $n$ unknown functions are pinned down by the single equation $L[y_p]=f$ , so you may impose $n-1$ auxiliary constraints for free. Choose them to annihilate the higher derivatives of the $u_i$ before they can appear, and substitution into $L$ collapses from a tangle of second derivatives to one clean algebraic line. What is left is a linear algebraic system for the derivatives $u_i'$ whose coefficient matrix is the Wronskian of the fundamental system — and because linear independence guarantees $W\neq0$ , that system is always solvable by Cramer's rule. Integrate the resulting $u_i'$ and you have a particular solution in closed form: an integral of the forcing against a kernel built purely from homogeneous data. That kernel is the Green's function in embryo, and the whole construction is the bridge from elementary ODE technique to boundary-value problems, Sturm–Liouville theory, and the impulse-response calculus of driven physical systems. Where undetermined coefficients stalls — on $\sec x$ , $\tan x$ , $1/x$ , a measured drive — variation of parameters simply works, at the honest price of an integral.

Key takeaways

Variation of parameters builds a particular solution of an inhomogeneous linear ODE by letting the constants of the homogeneous general solution become unknown functions $u_i(x)$ .
For $y''+py'+qy=f$ with fundamental pair $y_1,y_2$ , write $y_p=u_1y_1+u_2y_2$ and impose the gauge constraint $u_1'y_1+u_2'y_2=0$ to reduce the order.
Cramer's rule then gives the master formulae $u_1'=-y_2f/W$ and $u_2'=y_1f/W$ , with $W$ the Wronskian — always nonzero for an independent pair, so the method never fails to be solvable.
It handles any continuous forcing ( $\sec x$ , $\tan x$ , $1/x$ , tabulated data), unlike undetermined coefficients which needs exponential–polynomial drives.
The integrated answer is a convolution of $f$ against a kernel of homogeneous solutions — the Green's function of the problem, and the impulse (Duhamel) response of a driven oscillator.
Key idea: promote the integration constants to functions, impose one fewer constraint than the order, and the inhomogeneous equation collapses to a linear system for the $u_i'$ whose coefficient matrix is the Wronskian.

The core insight

Promote the integration constants of the homogeneous solution to functions, then impose just enough auxiliary constraints — one fewer than the order — so that the inhomogeneous equation reduces to a solvable linear system for the derivatives $u_i'$ , whose coefficient matrix is the Wronskian.

Where this lives today

The method is quiet infrastructure, present whenever a system is driven. In classical and quantum mechanics the Duhamel integral $x_p(t)=\frac{1}{m\omega_0}\int_0^t\sin\omega_0(t-\tau)F(\tau)\,\mathrm{d}\tau$ is the exact response of an oscillator to any force history — the workhorse behind vibration analysis, control theory, and the impulse-response formalism of linear systems. In electromagnetism and acoustics the same kernel, upgraded to satisfy boundary conditions, is the Green's function whose numerical evaluation underlies boundary-element solvers. In quantum scattering the radial Schrödinger equation is solved by treating the free spherical waves as the fundamental system and the potential as forcing, and variation of parameters reproduces the Lippmann–Schwinger integral equation that is the foundation of scattering theory and the Born series. Where the homogeneous solutions are special functions — Bessel, Legendre, hypergeometric — the tabulated Wronskians (for instance $W[J_\nu,Y_\nu]=2/(\pi x)$ ) turn the method into closed-form integral representations of inhomogeneous special-function equations. As a computational primitive it is the analytic backbone of Sturm–Liouville eigenfunction expansions and of the Green's-function machinery in every graduate methods course.

The physical picture

Picture the fundamental solutions $y_1,y_2$ as a moving frame in solution space. With constant coefficients, $c_1y_1+c_2y_2$ sweeps out every possible free motion of the system and nothing more. To follow a forced motion you let the coefficients drift, $u_1(x),u_2(x)$ , tuning them continuously so the trajectory tracks the driving force at every instant. But two drifting functions carry one redundant degree of freedom, because they must satisfy only the single equation $L[y_p]=f$ . Lagrange spends that redundancy on a deliberate gauge choice — the constraint $u_1'y_1+u_2'y_2=0$ — picked precisely so that the second derivatives of the $u_i$ never appear. What survives after substitution is a $2\times2$ linear system whose determinant is the Wronskian, sitting in the denominator of the answer. The first figure traces that logic as a pipeline: forcing in, particular solution out, Wronskian at the pivot.

Figure

Variation of parameters as a pipeline: promote the constants to functions, impose the gauge constraint, and the operator equation collapses to a 2×2 linear system whose determinant is the Wronskian W = y₁y₂' − y₁'y₂, solved by Cramer's rule and then integrated

The pipeline is abstract; the payoff is concrete. Feed the method a forcing that undetermined coefficients cannot touch — say $\sec x$ driving the unit oscillator $y''+y=\sec x$ — and out comes $y_p=\cos x\,\ln|\cos x|+x\sin x$ . The bounded piece $\cos x\,\ln|\cos x|$ is the routine response; the term $x\sin x$ is the fingerprint of resonance, a secular growth that appears because $\sec x$ overlaps the natural frequency $\omega=1$ . The second figure plots that solution, with real numbers, so you can watch the amplitude climb.

Figure

The worked answer for y'' + y = sec x is yₚ = cos x·ln(cos x) + x sin x. It climbs from 0 to yₚ(1.0) = 0.509 and yₚ(1.5) = 1.309, tracking the linear envelope y = x: the x sin x piece is resonant secular growth from driving at the natural frequency ω = 1

The core equations

Everything begins from the equation in normalized form — leading coefficient one — because the master formulae are derived by substituting into the operator $y''+py'+qy$ with a unit coefficient on $y''$ . Any deviation from this normalization is the single most common source of error in the whole method.

y''+p(x)\,y'+q(x)\,y=f(x) \qquad (1)

where: $p(x),q(x)$ — continuous coefficient functions, with the leading coefficient already divided out; $f(x)$ — the continuous forcing or source term; $y$ — the unknown function of $x$ . Reading it: the left side is the linear operator $L$ ; the right side is what drives it away from free motion. The construction rests on knowing a fundamental system $y_1,y_2$ of the homogeneous equation. The ansatz then replaces each constant by a function, and this is the one genuinely creative act — everything after it is forced.

y_p(x)=u_1(x)\,y_1(x)+u_2(x)\,y_2(x) \qquad (2)

where: $y_1,y_2$ — a fundamental (linearly independent) pair solving $y''+py'+qy=0$ ; $u_1,u_2$ — the unknown functions that have replaced the integration constants. Reading it: two unknown functions where before there were two numbers. Two unknown functions constrained by one equation leave a degree of freedom to spend. We spend it on a gauge choice designed to make the algebra collapse — set the "coefficient-velocity" contribution to $y_p'$ to zero.

u_1'\,y_1+u_2'\,y_2=0 \qquad (3)

where: $u_i'$ — the derivative of the $i$ -th coefficient function. Reading it: the part of $y_p'$ coming from the changing coefficients is switched off, so $y_p'$ looks formally like a homogeneous combination. With the constraint imposed, differentiating once more and substituting $y_p$ into $L$ makes every term multiplied by an undifferentiated $u_i$ vanish (because each $y_i$ solves the homogeneous equation), leaving a single surviving line.

u_1'\,y_1'+u_2'\,y_2'=f \qquad (4)

where: $f$ — the forcing that the surviving first-derivative terms must reproduce. Reading it: this is the second algebraic equation, partnering the constraint (3) to pin down both $u_i'$ . Equations (3) and (4) form a $2\times2$ linear system in $u_1',u_2'$ . Its coefficient matrix is exactly the Wronskian matrix, and its determinant is the Wronskian — the quantity that decides solvability.

W(x)=\begin{vmatrix} y_1 & y_2\\ y_1' & y_2'\end{vmatrix}=y_1 y_2'-y_2 y_1' \qquad (5)

where: $W(x)$ — the Wronskian of the fundamental pair. Reading it: by Abel's identity $W$ is either everywhere zero or nowhere zero on the interval, and for an independent pair it is nowhere zero — so the system is always invertible. Cramer's rule on the system delivers the master formulae, the two equations worth memorizing.

u_1'=-\frac{y_2\,f}{W},\qquad u_2'=\frac{y_1\,f}{W} \qquad (6)

where: $u_1',u_2'$ — the derivatives obtained by inverting the $2\times2$ system. Reading it: the forcing $f$ is weighted by the other solution and divided by the Wronskian; the sign flips between the two. Integrating (6) reconstructs the coefficient functions and hence the particular solution in indefinite-integral form.

y_p(x)=-y_1(x)\!\int\!\frac{y_2\,f}{W}\,\mathrm{d}x+y_2(x)\!\int\!\frac{y_1\,f}{W}\,\mathrm{d}x \qquad (7)

where: $\int$ — indefinite integral, whose constants of integration merely add homogeneous pieces $c_1y_1+c_2y_2$ and are absorbed into the general solution. Reading it: two quadratures finish the job. Choosing definite limits from a base point $x_0$ fuses the two integrals into a single convolution against a kernel built only from homogeneous data — the Green's-function form.

y_p(x)=\int_{x_0}^{x}\frac{y_1(t)\,y_2(x)-y_1(x)\,y_2(t)}{W(t)}\,f(t)\,\mathrm{d}t \qquad (8)

where: $x_0$ — base point, chosen so that $y_p(x_0)=y_p'(x_0)=0$ ; $t$ — the integration (source) variable. Reading it: the response at $x$ is the accumulated effect of the source at every earlier point $t$ , propagated by the kernel $G(x,t)$ . The whole scheme lifts verbatim to order $n$ : impose $n-1$ constraints that zero the lower-derivative combinations, and one final equation carries the forcing, all solved by Cramer's rule on the full $n\times n$ Wronskian.

\sum_{k=1}^{n}u_k'\,y_k^{(j)}=0\ \ (j=0,\dots,n-2),\qquad \sum_{k=1}^{n}u_k'\,y_k^{(n-1)}=f \qquad (9)

where: $y_k$ — the $n$ fundamental solutions; $y_k^{(j)}$ — the $j$ -th derivative of $y_k$ ; $u_k'$ — solved via Cramer's rule with the full $n\times n$ Wronskian in the denominator. Reading it: the second-order case is just $n=2$ .

Derivation AFrom the ansatz to the Wronskian system

The goal is to show, with no gaps, why the single operator equation $L[y_p]=f$ is equivalent to the tidy pair of algebraic equations (3) and (4), and hence why the master formulae (6) follow by Cramer's rule. We start from the ansatz $y_p=u_1y_1+u_2y_2$ , differentiate twice, deploy the gauge constraint at exactly the right moment to kill the second derivatives of the $u_i$ , and lean on the fact that each $y_i$ annihilates the operator. The result is that two functions constrained by one equation are made determinate by one freely chosen extra equation.

Step 1. Differentiate the ansatz once: $y_p'=(u_1'y_1+u_2'y_2)+(u_1y_1'+u_2y_2')$ . The first bracket carries the derivatives of the unknown coefficients; the second looks like an ordinary homogeneous derivative. Why: isolating the $u_i'$ terms is what lets us cancel them by a choice, rather than being forced to carry $u_i''$ later.

Step 2. Impose the gauge constraint $u_1'y_1+u_2'y_2=0$ , so that $y_p'=u_1y_1'+u_2y_2'$ . We are entitled to one extra condition because two functions answer to only one equation; this particular one is chosen so that $u_i''$ can never appear. Why: without it, the next differentiation would produce $u_1''y_1+u_2''y_2$ , and the problem would be harder than the one we started with.

Step 3. Differentiate again: $y_p''=(u_1'y_1'+u_2'y_2')+(u_1y_1''+u_2y_2'')$ . Only first derivatives of the $u_i$ survive, exactly as the gauge choice intended. Why: the constraint has done its work — the second derivatives of the coefficients are gone for good.

Step 4. Substitute $y_p,y_p',y_p''$ into $L[y]=y''+py'+qy$ and regroup by $u_1$ and $u_2$ : $L[y_p]=u_1(y_1''+py_1'+qy_1)+u_2(y_2''+py_2'+qy_2)+(u_1'y_1'+u_2'y_2')$ . Why: collecting terms exposes that the bracketed factors are just $L[y_1]$ and $L[y_2]$ .

Step 5. Each $y_i$ solves the homogeneous equation, so $y_i''+py_i'+qy_i=0$ and both bracketed groups vanish, leaving $L[y_p]=u_1'y_1'+u_2'y_2'$ . Setting this equal to $f$ gives the second equation (4). Why: this is the whole magic — the operator, applied to the ansatz, distils to a single first-order combination of the $u_i'$ .

Step 6. Equations (3) and (4) are a linear system $\big(\begin{smallmatrix}y_1&y_2\\ y_1'&y_2'\end{smallmatrix}\big)\big(\begin{smallmatrix}u_1'\\ u_2'\end{smallmatrix}\big)=\big(\begin{smallmatrix}0\\ f\end{smallmatrix}\big)$ with determinant $W$ . Cramer's rule gives $u_1'=-y_2f/W$ and $u_2'=y_1f/W$ . Why: $W\neq0$ for an independent pair, so the inversion is legitimate everywhere on the interval — the method cannot break down.

Go deeper · Derivation B · 1 minThe kernel and its built-in zero initial conditions

Derivation BThe kernel and its built-in zero initial conditions

The aim is to turn the indefinite-integral answer (7) into a definite integral against a single kernel $G(x,t)$ , and to prove that this particular choice of limits makes the particular solution satisfy homogeneous initial conditions $y_p(x_0)=y_p'(x_0)=0$ . This is what promotes variation of parameters from a bookkeeping recipe into the Green's-function representation used throughout boundary-value and Sturm–Liouville theory. We fix a base point $x_0$ , combine the two integrals, and differentiate under the integral sign with the Leibniz rule.

Step 1. Replace the indefinite integrals in (7) with definite ones from $x_0$ to $x$ : $y_p=-y_1(x)\!\int_{x_0}^{x}\frac{y_2(t)f(t)}{W(t)}\mathrm{d}t+y_2(x)\!\int_{x_0}^{x}\frac{y_1(t)f(t)}{W(t)}\mathrm{d}t$ . Why: definite limits fix the arbitrary constants, selecting one specific particular solution rather than a family.

Step 2. The factors $y_1(x),y_2(x)$ are constants with respect to $t$ , so pull them inside and combine: $y_p=\int_{x_0}^{x}\frac{y_1(t)y_2(x)-y_1(x)y_2(t)}{W(t)}f(t)\,\mathrm{d}t\equiv\int_{x_0}^{x}G(x,t)f(t)\,\mathrm{d}t$ . Why: writing the integrand as one kernel $G(x,t)$ makes the causal structure explicit — the source at $t$ acts through $G$ .

Step 3. Evaluate at $x=x_0$ : the integration range collapses, $y_p(x_0)=\int_{x_0}^{x_0}(\cdots)\,\mathrm{d}t=0$ . Why: the first homogeneous initial condition is automatic, with no extra algebra.

Step 4. Differentiate with the Leibniz rule: $y_p'(x)=G(x,x)f(x)+\int_{x_0}^{x}\partial_xG(x,t)f(t)\,\mathrm{d}t$ , and note $G(x,x)=\frac{y_1(x)y_2(x)-y_1(x)y_2(x)}{W(x)}=0$ . At $x=x_0$ the remaining integral is again empty, so $y_p'(x_0)=0$ . Why: the diagonal value $G(x,x)=0$ is precisely the property that keeps the boundary term from spoiling the derivative — it is the signature of a second-order Green's function.

Go deeper · Derivation C · 1 minThe Duhamel impulse response of a driven oscillator

Derivation CThe Duhamel impulse response of a driven oscillator

The purpose here is to specialize the kernel to the constant-coefficient oscillator and recover the impulse-response (Duhamel) integral that physicists use for arbitrarily forced oscillators — step forces, ramps, measured or noisy drives — where undetermined coefficients is powerless. We take $y_1=\cos\omega_0 t$ , $y_2=\sin\omega_0 t$ , compute the Wronskian, and read off the convolution kernel, showing it is the free response to a unit impulse.

Step 1. For $\ddot x+\omega_0^2x=F(t)/m$ the homogeneous solutions are $y_1=\cos\omega_0 t$ and $y_2=\sin\omega_0 t$ , with forcing $f=F/m$ . Why: these are the free oscillations, the natural building blocks the method demands.

Step 2. Compute the Wronskian: $W=\cos\omega_0 t\,(\omega_0\cos\omega_0 t)-\sin\omega_0 t\,(-\omega_0\sin\omega_0 t)=\omega_0$ . It is constant, as Abel's identity predicts for an equation with no first-derivative term. Why: a constant $W$ makes the kernel a pure function of $t-\tau$ , which is the hallmark of a time-invariant system.

Step 3. Insert into the kernel (8): $G(t,\tau)=\frac{\cos\omega_0\tau\,\sin\omega_0 t-\cos\omega_0 t\,\sin\omega_0\tau}{\omega_0}=\frac{\sin\omega_0(t-\tau)}{\omega_0}$ using the sine subtraction identity. Why: the kernel depends only on the elapsed time $t-\tau$ — it is the memory of the oscillator.

Step 4. Assemble the particular solution: $x_p(t)=\frac{1}{m\omega_0}\int_0^t\sin\omega_0(t-\tau)\,F(\tau)\,\mathrm{d}\tau$ . This is the Duhamel integral; $G$ is literally the response to a delta-function impulse $F(\tau)=m\,\delta(\tau)$ . Why: it recasts "solve the ODE" as "convolve the force with the impulse response", the working definition of a Green's function in physics. *The kernel $\sin\omega_0(t-\tau)/\omega_0$ satisfies $\ddot G+\omega_0^2 G=\delta(t-\tau)$ with $G=0$ and $\partial_tG=1$ at $t=\tau$ — the defining conditions of the impulse-response Green's function.*

Limiting cases

Variation of parameters is the general case from which several familiar methods drop out as limits, and into which physical intuition can be read directly.

Case	Result	Meaning
Only one solution known	$y=u(x)\,y_1(x)$ , one constraint	The method degenerates to reduction of order; a single unknown $u$ and a first-order equation for $u'$ reconstruct $y_2$ en route.
Exponential–polynomial forcing	$f(x)=e^{\alpha x}P_n(x)$	Undetermined coefficients is faster; variation of parameters still works but demands integrals you could have guessed the form of.
Constant coefficients, sinusoidal drive	$y''+\omega_0^2y=\cos\omega t$	Reproduces the standard resonance curve; at $\omega\to\omega_0$ the integral exposes the secular $t\sin\omega_0 t$ growth automatically.
Boundary conditions imposed	$y_p(x)=\int_a^b G(x,t)\,f(t)\,\mathrm{d}t$	Choosing $y_1,y_2$ to satisfy the left and right conditions turns the kernel into the Green's function of the two-point boundary-value problem.
Time-invariant system	$G(t,\tau)=G(t-\tau)$	A constant Wronskian makes the kernel a convolution — the Duhamel impulse response of a driven physical system.

Common trapForgetting to normalize the leading coefficient

The master formulae $u_1'=-y_2f/W$ , $u_2'=y_1f/W$ assume the equation is in the normalized form $y''+py'+qy=f$ , i.e. the coefficient of $y''$ is exactly $1$ . If you begin from $a(x)\,y''+b(x)\,y'+c(x)\,y=g(x)$ , you must first divide through by $a(x)$ , so the forcing that enters the formulae is $f=g/a$ , not $g$ . Skipping this is the single most frequent mistake, and it is silent — the algebra still runs, but every answer is off by a factor of $a(x)$ . In the Cauchy–Euler example below, the forcing is $x^3$ before normalization and $f=x$ after dividing by $x^2$ ; using $x^3$ would give a wrong $y_p$ by a factor of $x^2$ . Always normalize first, then read off $f$ .

Worked examples

1 · The secant-driven oscillator (fully worked, with check)

Solve $y''+y=\sec x$ — a forcing undetermined coefficients cannot handle, since $\sec x$ is not of exponential–polynomial type.

Step 1. The homogeneous equation $y''+y=0$ has fundamental pair $y_1=\cos x$ , $y_2=\sin x$ . Wronskian: $W=\cos x\,(\cos x)-\sin x\,(-\sin x)=\cos^2x+\sin^2x=\mathbf{1}$ . The equation is already normalized, so $f=\sec x$ .

Step 2. Master formula for $u_1'$ : $u_1'=-\dfrac{y_2f}{W}=-\sin x\,\sec x=-\tan x$ . Integrate: $u_1=\int(-\tan x)\,\mathrm{d}x=\ln|\cos x|$ , using $\int\tan x\,\mathrm{d}x=-\ln|\cos x|$ .

Step 3. Master formula for $u_2'$ : $u_2'=\dfrac{y_1f}{W}=\cos x\,\sec x=1$ , which integrates to $u_2=\int1\,\mathrm{d}x=x$ . The neat collapse to a constant is a small gift of this particular forcing.

Step 5. Check by substitution. $y_p'=-\sin x\,\ln|\cos x|+x\cos x$ (the $\pm\sin x$ terms cancel), and $y_p''=-\cos x\,\ln|\cos x|+\dfrac{\sin^2x}{\cos x}+\cos x-x\sin x$ . Then $y_p''+y_p=\dfrac{\sin^2x}{\cos x}+\cos x=\dfrac{\sin^2x+\cos^2x}{\cos x}=\dfrac{1}{\cos x}=\sec x$ . **Check: $y_p''+y_p=\sec x$ exactly.* The term $x\sin x$ is the secular signature of driving at the natural frequency $\omega=1$ ; it is the unbounded growth plotted in the second figure.

2 · Cauchy–Euler equation, and the normalization trap

Solve $x^2y''-2xy'+2y=x^3$ for $x>0$ .

Step 1. The homogeneous Cauchy–Euler equation is solved by $y=x^r$ : substituting gives $r(r-1)-2r+2=r^2-3r+2=(r-1)(r-2)=0$ , so $r=1,2$ and $y_1=x$ , $y_2=x^2$ .

Step 2. Normalize first. Divide by $x^2$ : $y''-\dfrac{2}{x}y'+\dfrac{2}{x^2}y=x$ , so the forcing seen by the formulae is $\mathbf{f=x}$ — not $x^3$ . This is the trap from the warning box.

Step 3. Wronskian: $W=\begin{vmatrix}x&x^2\\1&2x\end{vmatrix}=2x^2-x^2=x^2$ , nonzero for $x\neq0$ , confirming independence.

Step 4. Master formulae: $u_1'=-\dfrac{y_2f}{W}=-\dfrac{x^2\cdot x}{x^2}=-x$ and $u_2'=\dfrac{y_1f}{W}=\dfrac{x\cdot x}{x^2}=1$ . Integrate: $u_1=-\tfrac12x^2$ , $u_2=x$ .

Step 5. Assemble: $y_p=u_1y_1+u_2y_2=-\tfrac12x^2\cdot x+x\cdot x^2=-\tfrac12x^3+x^3=\tfrac12x^3$ . The particular solution is $y_p=\tfrac12x^3$ .

Step 6. Check against the original (un-normalized) equation: $x^2(3x)-2x(\tfrac32x^2)+2(\tfrac12x^3)=3x^3-3x^3+x^3=x^3$ . Verified. The check uses $x^3$ on the right because we back-substitute into the original equation, while the formula used $f=x$ — the factor of $x^2$ is exactly the normalization.

Physics of the portrait

Variation of parameters is where three centuries of mechanics and analysis converge: a computational trick born in the perturbation theory of the planets, sharpened into a general method for linear equations, and eventually recognized as the classical shadow of the Green's function.

Origins · 1 minJoseph-Louis Lagrange and how this was discovered

Joseph-Louis Lagrange

1736–1813 · variation of parameters and analytical mechanics

Joseph-Louis Lagrange gave the method its name and its rigour. Working on the slow, cumulative perturbations of planetary orbits — problems where the "constants" of a Keplerian ellipse drift over centuries under the tug of other planets — he developed the variation of arbitrary constants (la variation des constantes arbitraires), promoting orbital elements to slowly varying functions of time. Euler had already varied constants in specific perturbation calculations from the 1740s, but it was Lagrange who, in memoirs of the 1770s and then decisively in his celestial-mechanics work around 1808–1810, turned the device into a systematic algorithm and extended it to general linear differential equations. The auxiliary constraint that makes the method collapse so cleanly is precisely his — a deliberate choice among the redundant freedoms, the same instinct for economical structure that runs through his Mécanique analytique of 1788, where he banished geometric diagrams in favour of pure analysis. The coefficient matrix of his linear system would later be named the Wronskian, after Józef Hoene-Wroński, and the kernel his integrals produce would be recognized as George Green's function. What Lagrange bequeathed is the enduring idea that the solution of the forced problem is encoded entirely in the free problem — you need only let its constants breathe.

History

1748

— Euler varies orbital constants in planetary perturbation calculations

1774

— Lagrange develops variation of arbitrary constants for celestial mechanics

1788

— Lagrange publishes Mécanique analytique, recasting mechanics as pure analysis

1808

— Lagrange systematizes variation of constants for general linear equations

1812

— Hoene-Wroński introduces the determinant of a solution system

1882

— Muir names that determinant the Wronskian

Going further

Texts. Boyce & DiPrima, Elementary Differential Equations, §3.6 (variation of parameters, second order) · Arfken & Weber, Mathematical Methods for Physicists, 7th ed., Ch. 7 (ODEs) and Ch. 10 (Green's functions) · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed., §15.2 · Stone & Goldbart, Mathematics for Physics, Ch. 1–2 (ODE methods and Green's functions) · Hassani, Mathematical Physics, Ch. 14 · Byron & Fuller, Mathematics of Classical and Quantum Physics, Ch. 2. The master formulae $u_1'=-y_2f/W$ , $u_2'=y_1f/W$ follow from Cramer's rule on the Wronskian system. The secant example $y''+y=\sec x\Rightarrow y_p=\cos x\ln|\cos x|+x\sin x$ is confirmed by direct substitution ( $y_p''+y_p=\sec x$ ). The Cauchy–Euler example $x^2y''-2xy'+2y=x^3\Rightarrow y_p=\tfrac12x^3$ is confirmed by back-substitution ( $3x^3-3x^3+x^3=x^3$ ). The Duhamel kernel $\sin\omega_0(t-\tau)/\omega_0$ and the Bessel Wronskian $W[J_\nu,Y_\nu]=2/(\pi x)$ are standard results.

Prerequisites. Linear Second-Order ODEs · The Wronskian and Linear Independence · Reduction of Order · First-Order Ordinary Differential Equations. Leads to. Two-Point Boundary-Value Problems · Sturm–Liouville Theory · Existence and Uniqueness Theorems for ODEs.

Ask the tutor: Bring me any inhomogeneous linear ODE and I will walk you through normalizing it, building the Wronskian, and turning the master formulae into a checked particular solution.

When should I use variation of parameters instead of undetermined coefficients?

Use undetermined coefficients when the forcing is of exponential–polynomial–sinusoidal type ( $e^{\alpha x}$ , $x^n$ , $\cos\omega x$ , and products thereof), because guessing a trial form and matching coefficients is faster and avoids integration. Reach for variation of parameters the moment the forcing falls outside that family — $\sec x$ , $\tan x$ , $1/x$ , $\ln x$ , or tabulated/measured data — or when the coefficients are variable rather than constant. Variation of parameters is strictly more general: it works for any continuous forcing once the homogeneous solutions are known, at the cost of possibly non-elementary integrals. In short, undetermined coefficients is the specialist shortcut; variation of parameters is the universal tool.

Why can I choose the auxiliary constraint so freely?

Because the ansatz $y_p=u_1y_1+u_2y_2$ introduces two unknown functions but the original ODE imposes only one equation on them. That mismatch leaves exactly one degree of freedom, and you are entitled to spend it however you like without changing the set of achievable $y_p$ . Lagrange's choice $u_1'y_1+u_2'y_2=0$ is not the only legal constraint — infinitely many would do — but it is the one engineered so that the second derivatives of the coefficients never appear, leaving a first-order algebraic system instead of a harder second-order one. A different constraint would still give a correct particular solution, just via messier algebra. The freedom is real; the constraint is a smart convention, not a law.

How is variation of parameters related to reduction of order?

Reduction of order is variation of parameters with a single known solution. If you know only $y_1$ , you can still write $y=u(x)y_1(x)$ — the ansatz with one drifting coefficient — and because there is only one unknown function, there is no redundant freedom and no auxiliary constraint to impose. Substituting into the ODE yields a first-order equation for $u'$ , which you solve to recover the missing second solution $y_2$ and, with the same stroke, a particular solution. So the two methods are one idea at different ranks: reduction of order is the rank-one case, variation of parameters the full-rank generalization once a complete fundamental system is available.

Why does a periodic forcing like sec x produce a non-periodic answer?

Because $\sec x$ contains spectral content that overlaps the oscillator's natural frequency. Driving a lossless oscillator at (or overlapping) its resonance pumps energy in coherently, and the amplitude grows without bound — the mathematical signature is a secular term, here $x\sin x$ , whose envelope rises linearly with $x$ . This is the same mechanism that makes $y''+\omega_0^2y=\cos\omega_0 t$ produce $\tfrac{t}{2\omega_0}\sin\omega_0 t$ . The bounded part of the answer, $\cos x\ln|\cos x|$ , is the ordinary forced response; the unbounded $x\sin x$ part is resonance made visible. Real systems have damping, which caps the growth, but the undamped model exposes the resonant tendency cleanly.

What happens to the arbitrary constants from the indefinite integrals?

Nothing harmful. When you integrate $u_1'$ and $u_2'$ you may add arbitrary constants $A$ and $B$ ; these contribute $A\,y_1+B\,y_2$ to $y_p$ , which is just another homogeneous solution. Since the general solution already carries a complementary function $c_1y_1+c_2y_2$ , those extra pieces merge into it and shift the meaning of $c_1,c_2$ but change nothing physical. This is why textbooks drop the constants when computing $u_1,u_2$ : any specific antiderivative gives a valid particular solution, and the family of all particular solutions differs only by homogeneous terms. Choosing definite limits from a base point $x_0$ is the systematic way to fix them, and it is what yields the Green's-function kernel with zero initial conditions.

Does the method extend to systems and higher-order equations?

Yes, cleanly. For an $n$ -th order equation you promote all $n$ constants to functions and impose $n-1$ constraints that zero the successive lower-derivative combinations $\sum_k u_k'y_k^{(j)}=0$ for $j=0,\dots,n-2$ ; the final equation $\sum_k u_k'y_k^{(n-1)}=f$ carries the forcing. The resulting $n\times n$ linear system has the full Wronskian as its determinant and is solved by Cramer's rule. For first-order systems $\mathbf{y}'=A(x)\mathbf{y}+\mathbf{f}$ the same idea reads $\mathbf{y}_p=\Phi(x)\!\int\Phi^{-1}(t)\mathbf{f}(t)\,\mathrm{d}t$ , where $\Phi$ is a fundamental matrix — the matrix Wronskian — and $\Phi^{-1}$ plays the role of dividing by $W$ . The structure is identical at every rank.