Undetermined coefficients breaks when the forcing gets ugly. This is the universal fallback — built from scratch, from pure logic. The constants become functions. Everything else follows.
We want the full solution to a nonhomogeneous second-order linear ODE:
We already know how to solve the associated homogeneous version (when the right side is zero). Suppose those two linearly independent solutions are \(y_1\) and \(y_2\). Then the homogeneous solution is:
The question is: how do we find a particular solution \(y_p\) to the full forced equation?
Undetermined coefficients works by guessing the shape of \(y_p\). That's great when differentiating the forcing function keeps producing the same "family" of functions — like \(e^x\), \(\sin x\), polynomials. But what happens when you differentiate \(\tan x\)?
You get \(\sec^2 x\). Differentiate that and you get something else entirely. The family never closes. There's no finite set of shapes to guess from. The photocopier is broken.
Forcing functions that generate infinitely many new shapes when differentiated:
\(\tan x,\quad \sec x,\quad \ln x,\quad \dfrac{1}{x},\quad e^x \ln x\)
Any continuous forcing function \(g(x)\). No restrictions on form. Universal method.
Here is the core insight. Think about the homogeneous solution:
The constants \(c_1, c_2\) are fixed numbers. They don't change. That's exactly what makes it the homogeneous solution — the fixed weights just combine two natural modes of the system.
But what if the system is being driven by an external force \(g(x)\)? A fixed weight can't respond to a moving driver. So the idea is:
For the forced equation, let those constants become functions of \(x\):
\[ c_1 \;\to\; u_1(x), \qquad c_2 \;\to\; u_2(x) \]So we try \(y_p = u_1(x)\,y_1(x) + u_2(x)\,y_2(x)\) — the same structure as the homogeneous solution, but with adaptive weights. The structure absorbs the forcing. The weights track it.
This is not a guess pulled from thin air. It is motivated by the structure of the problem. The homogeneous basis \(\{y_1, y_2\}\) already spans all solutions to the unforced system — so any particular solution to the forced system should live in the same neighborhood, just with coefficients that adapt to the forcing.
We want to find \(u_1(x)\) and \(u_2(x)\) such that
satisfies \(y'' + p(x)y' + q(x)y = g(x)\). To do that, we differentiate \(y_p\) twice and substitute.
Apply the product rule on every term:
Notice we have four terms. If we differentiate again as-is, the next derivative will have eight terms, and the algebra becomes a complete swamp. So we impose a deliberate simplifying condition:
We freely impose:
\[ u_1' y_1 + u_2' y_2 = 0 \]This is not a constraint forced on us by the physics — it is a choice we make to keep the algebra clean. We have two unknown functions \(u_1, u_2\), which means we need two equations to determine them. The ODE will give us one. This condition is the second. We're using our freedom wisely.
With that condition, the first derivative simplifies beautifully:
Now differentiate the simplified \(y_p'\):
So now we have a clean inventory of all three quantities we need:
Plug everything into \(y'' + py' + qy = g\):
The magic just happened. Because \(y_1\) and \(y_2\) are solutions to the homogeneous equation, the giant grouped expressions vanish identically. All the structure collapses, and we're left with a simple second equation.
We know \(y_1'' + py_1' + qy_1 = 0\) because \(y_1\) is a homogeneous solution — it was built to satisfy exactly that. Same for \(y_2\). So when we substitute \(y_p\) into the full ODE, those terms are identically zero. The forcing function \(g(x)\) is the only thing left.
We now have exactly two equations for two unknowns \(u_1'(x)\) and \(u_2'(x)\):
This is a straightforward linear system. Write it in matrix form:
The determinant of the coefficient matrix is exactly the Wronskian:
Think of \(y_1\) and \(y_2\) as two vectors in function space. The Wronskian measures how "different" they are — their linear independence. If \(W \neq 0\), the two functions span a 2D space and you can uniquely write any solution as a combination of them. If \(W = 0\), they're essentially the same solution in disguise, and the matrix can't be inverted.
This is precisely why \(W \neq 0\) is required: we need to invert the system. No independence → no inversion → no solution. Game over.
To solve a 2×2 system \(A\mathbf{x} = \mathbf{b}\), Cramer's rule says: replace each column of \(A\) with the right-hand side vector \(\mathbf{b}\) and divide by \(\det(A)\).
For \(u_1'\) — replace the first column:
For \(u_2'\) — replace the second column:
\(u_1'\) has a minus sign because the right-hand side \((0, g)^T\) hits column 1's position differently in the determinant expansion. \(u_2'\) does not. This asymmetry comes from the Cramer's rule expansion — it is baked into the structure, not a typo.
Now integrate \(u_1'\) and \(u_2'\) to get the weight functions:
Then the particular solution is:
And the full general solution is:
Valid when: \(y'' + p(x)y' + q(x)y = g(x)\), \(W(y_1,y_2) \neq 0\) on the interval
You might be skeptical: we just imposed the condition \(u_1' y_1 + u_2' y_2 = 0\). Isn't that cheating?
No. Here's the key insight:
We introduced two unknown functions \(u_1(x)\) and \(u_2(x)\). To uniquely determine two functions, we need exactly two equations. The ODE supplies one (after substitution). The gauge condition is the second — it's a free choice we make to tame the algebra.
This is the same spirit as choosing a coordinate system in physics. The physics doesn't care which coordinates you use — you choose the ones that make the problem cleanest. The gauge condition is our clever coordinate choice.
Without it, variation of parameters would still work in principle, but differentiating \(y_p'\) the second time would produce many extra terms involving \(u_1''\) and \(u_2''\), and the system would become far more complicated. The condition eliminates this mess without changing the solution space.
In physics, this kind of freedom-fixing is called a gauge choice — the same idea underlies how you fix the gauge in electromagnetism. The name "variation of parameters" could almost be "variation with gauge freedom."
When you see \(y'' + p(x)y' + q(x)y = g(x)\), do this:
Make sure the coefficient of \(y''\) is exactly 1. If it's \(a(x)\), divide through first. The \(g(x)\) in the formula is \(r(x)/a(x)\), not \(r(x)\).
Find two linearly independent solutions \(y_1\) and \(y_2\) from the characteristic equation or known methods.
\(W = y_1 y_2' - y_2 y_1'\). Check it's nonzero on your interval.
\(u_1' = -y_2 g / W\) and \(u_2' = y_1 g / W\). Don't drop the minus sign on \(u_1'\).
\(u_1 = \int u_1'\,dx\) and \(u_2 = \int u_2'\,dx\). Drop constants of integration — they'll merge into \(c_1, c_2\) in \(y_h\).
Simplify. If any pieces of \(y_p\) are already in \(y_h\), absorb them — they're redundant.
Include all arbitrary constants \(c_1, c_2\) from \(y_h\). Apply initial conditions if given.
This example is the canonical reason variation of parameters exists. The forcing function \(\tan x\) has no closed finite derivative family.
Step A — Identify Standard Form
Compare with \(y'' + p(x)y' + q(x)y = g(x)\):
The coefficient of \(y''\) is already 1. We're in standard form.
Step B — Homogeneous Solution
Solve \(y'' + y = 0\). Characteristic equation:
So \(y_1 = \cos x\) and \(y_2 = \sin x\).
Step C — Wronskian
The Wronskian equals 1 everywhere. This is as clean as it gets.
Step D — Compute \(u_1'\)
Now integrate:
Step E — Compute \(u_2'\)
Step F — Build \(y_p\)
Final Answer
The log term comes from the integral of \(\sec x\) — exactly the kind of "exotic" term that undetermined coefficients would never produce. This is variation of parameters being necessary, not just convenient.
This example shows something slick: the forcing function \(e^x\) is already part of the homogeneous solution. Undetermined coefficients would require a resonance fix (multiply the guess by \(x\)). Variation of parameters handles this automatically — no special casing needed.
Step A — Homogeneous Solution
So \(y_1 = e^x\) and \(y_2 = e^{-x}\). Notice \(g(x) = e^x = y_1\). This is resonance territory.
Step B — Wronskian
Step C — Compute \(u_1'\)
Step D — Compute \(u_2'\)
Step E — Build \(y_p\)
The term \(-\tfrac{1}{4}e^x\) is already a homogeneous solution (it's a multiple of \(y_1 = e^x\)). It belongs in \(y_h\), not \(y_p\). Since \(c_1 e^x\) already captures it, we can absorb this term there — leaving the cleanest possible particular solution.
Look at the form: \(x e^x\). This is exactly the "resonance guess" you'd have needed in undetermined coefficients. Variation of parameters produced it automatically from the integral \(\int \tfrac{1}{2}\,dx = \tfrac{x}{2}\).
Final Answer
In undetermined coefficients, resonance requires a special rule: if your guess overlaps with \(y_h\), multiply by \(x\). In variation of parameters, there is no special rule — the algebra produces the \(x\) factor naturally through the integration step. The method doesn't know or care about resonance; it just works.
This explorer solves \(y'' + \omega^2 y = A\sin(\beta x)\) using variation of parameters numerically. Adjust the parameters to see how the particular solution tracks the forcing.
━ \(y_h\) (homogeneous) ━ \(y_p\) (particular) ━ \(y\) (full solution)
The formulas \(u_1' = -y_2 g/W\) and \(u_2' = y_1 g/W\) require the equation to be in standard form with leading coefficient 1. If your equation is \(a(x)y'' + b(x)y' + c(x)y = r(x)\), you must divide through by \(a(x)\) first. The \(g\) in the formula is \(r(x)/a(x)\), not \(r(x)\).
The formula is \(u_1' = -y_2 g / W\). That minus sign is structural — it comes from the Cramer's rule expansion of the first column. \(u_2'\) has no minus. Every time you write these down, double-check the sign on \(u_1'\).
If \(y_1\) and \(y_2\) are linearly dependent (one is just a multiple of the other), the Wronskian is zero and the system cannot be solved. You need a genuine fundamental set — two genuinely different solutions spanning the solution space.
If any term in your \(y_p\) is of the form \(C \cdot y_1\) or \(C \cdot y_2\), absorb it into the \(c_1, c_2\) constants in \(y_h\). Particular solutions are not unique — you want the cleanest one. Keeping extra homogeneous pieces in \(y_p\) isn't wrong, but it's messy and confusing.
| Undetermined Coefficients | Variation of Parameters | |
|---|---|---|
| Speed | Fast when it works | More computation |
| Scope | Limited forcing functions | Any continuous \(g(x)\) |
| Resonance | Special rule required | Automatic |
| Variable coefficients | No | Yes (if \(y_1, y_2\) are known) |
Variation of parameters works because of this chain of logic:
① Standard form (leading coeff = 1)
② Identify \(y_1, y_2\) from homogeneous
③ Verify \(W \neq 0\)
④ Minus sign on \(u_1'\) — always
⑤ Drop integration constants (merge into \(c_1, c_2\))
⑥ Simplify \(y_p\) fully before adding \(y_h\)
Forcing is \(\tan x\), \(\sec x\), \(\ln x\), \(1/x\), or anything UC can't guess
Forcing has a finite derivative family: \(e^{ax}\), \(\sin bx\), \(\cos bx\), polynomials, products of these