Variation of Parameters in General

Why the quick methods run out

By this point in the rung you can solve a great many nth-order linear equations with forcing on the right. Spot the form of the forcing, write a trial solution with unknown coefficients — or, more systematically, hit the whole equation with an annihilator that kills the forcing — and out drops a particular solution. These methods are fast and they feel almost mechanical, which is exactly their charm. But notice the small print on every example: the forcing was always a polynomial, an exponential e^(rx), a sine or cosine, or a product of those.

That is no accident — it is the whole limitation. The annihilator method only works on forcings that *are themselves* solutions of some constant-coefficient homogeneous equation, because only those have an operator that annihilates them. So sec(x), tan(x), ln(x), 1/x, or x^x have no annihilator at all: no finite polynomial in D sends them to zero. And the instant a coefficient varies — a Cauchy-Euler equation, or anything wilder — there is no characteristic polynomial to begin with, so undetermined coefficients has nothing to stand on. You need a method that does not care what the forcing looks like.

The one idea, scaled up to order n

You already met variation of parameters at second order, and the central trick has not changed one bit. Start from a fundamental set y1, y2, ..., yn of the homogeneous equation, whose general solution is c1 y1 + ... + cn yn. The constants c1, ..., cn are dead weight — every choice of them still solves the homogeneous problem, never the forced one. So let the constants vary: replace each constant ck by an unknown function uk(x), and look for a particular solution of the shape y_p = u1 y1 + u2 y2 + ... + un yn. It is the very same instinct as reduction of order, where a single constant was promoted to a function; here all n of them are.

But now there are n unknown functions u1, ..., un, and a single equation L[y_p] = g(x) cannot pin down n functions — it is one constraint for n choices. The freedom is real, and the genius of the method is to *spend* that freedom wisely. As you differentiate y_p over and over, each derivative throws off a cluster of terms involving the uk'. Variation of parameters chooses, at every stage but the last, to make those clusters vanish: it imposes n-1 simplifying conditions that force the uk'-terms to cancel at each derivative level up to the (n-2)th. This keeps the expressions from exploding, and it is exactly enough extra equations to make the system solvable.

The Wronskian system that pops out

When the dust settles, those n-1 cancellation conditions plus the one genuine equation L[y_p] = g combine into a clean linear system for the n unknowns u1', ..., un'. And the coefficient matrix is not some new monster — it is precisely the Wronskian matrix of your fundamental set: the rows are y1...yn, then y1'...yn', on down to the (n-1)th derivatives in the bottom row. The right-hand side is almost all zeros, with the forcing sitting alone at the bottom. That structure is worth staring at: the same Wronskian that tested whether your solutions were independent now does the work of building the particular solution.

Solve for u1', ..., un' :   [ y1      y2     ...  yn     ] [u1']   [ 0      ]
                            [ y1'     y2'    ...  yn'    ] [u2'] = [ 0      ]
                            [  :       :           :     ] [ : ]   [  :     ]
                            [ y1^(n-1) y2^(n-1) ... yn^(n-1)] [un']   [ g(x)   ]

Cramer's rule:   uk' = Wk(x) / W(x)        ( W = Wronskian determinant )
Integrate:       uk  = INT [ Wk(x) / W(x) ] dx
Assemble:        y_p = u1*y1 + u2*y2 + ... + un*yn

( Wk = W with column k replaced by the column [0, 0, ..., g(x)] )

Solving by Cramer's rule makes the shape beautifully explicit: uk' = Wk(x) / W(x), where W is the Wronskian determinant and Wk is that same determinant with the k-th column swapped out for the column of zeros-then-g. Integrate each uk' to recover uk, build y_p, and you are done. Notice that W lives in every denominator — and it never vanishes precisely because y1, ..., yn form a fundamental set (independent solutions have a nonzero Wronskian). The very condition that made them a good basis is the condition that keeps this method from dividing by zero.

Walking through it once

Here is the recipe distilled to its steps. The first three are exactly the second-order procedure you already know, written for n functions instead of two; the only real cost of going higher is that the determinants grow and the bookkeeping gets heavier.

1. Put the equation in standard form (leading coefficient 1) and solve the homogeneous part first: find a fundamental set y1, ..., yn. Variation of parameters is a particular-solution machine — it assumes you already have the homogeneous solutions.
2. Set up the Wronskian-matrix system: rows of y's and their derivatives up to order n-1, right-hand side all zeros except g(x) in the last slot.
3. Solve for u1', ..., un' by Cramer's rule: uk' = Wk(x) / W(x), where W is the Wronskian and Wk replaces column k with the right-hand side.
4. Integrate each uk' to get uk. Drop the integration constants — they only rebuild the homogeneous solution, which you will add back separately.
5. Assemble y_p = u1 y1 + ... + un yn, then write the full answer y = (c1 y1 + ... + cn yn) + y_p — homogeneous plus particular, the same structure as ever.

A flavour of why it earns its keep: suppose you need a particular solution of y'' + y = sec(x). No annihilator touches sec(x), so undetermined coefficients is dead on arrival. But variation of parameters does not blink. The homogeneous solutions are y1 = cos(x), y2 = sin(x), with Wronskian W = 1, so u1' = -sin(x) sec(x) = -tan(x) and u2' = cos(x) sec(x) = 1. Integrating, u1 = ln|cos(x)| and u2 = x, giving y_p = cos(x) ln|cos(x)| + x sin(x). That logarithm and that x sin(x) are exactly the kind of answer the quick methods can never produce — and they fell out of one Wronskian system.

What it costs, and where it sends you next

Be clear-eyed about the price. First, variation of parameters guarantees a particular solution *as an integral* — it does not promise that integral is elementary. The formula uk' = Wk/W is always correct, but the integral of Wk/W may simply have no closed form, in which case the honest answer is the solution left as an integral. That is not a failure; an integral is a perfectly legitimate way to define a function. Second, the bookkeeping at order n is genuinely heavy: an n-by-n Wronskian, n determinants Wk, and n integrals. The method always works in principle; whether it works *pleasantly* depends on the equation.

It is worth naming what variation of parameters is *not*, too. Like everything in this rung, it lives entirely inside linear theory: the whole construction leans on superposition — on y_p being a combination of homogeneous solutions — which collapses the moment a nonlinear term like (y')^2 appears. And it is not a homogeneous-solver: it consumes a fundamental set, it does not produce one. For variable-coefficient equations that means you may still need reduction of order, a known special function, or series methods just to obtain the yk before this method can even begin.

With this guide the higher-order rung closes a complete circle. You can write down the n-dimensional solution space, read its basis off the degree-n characteristic polynomial, wield differential operators and the annihilator on the friendly forcings, tame the variable-coefficient Cauchy-Euler case, and now — with variation of parameters — handle *any* integrable forcing whatsoever. The next rungs trade these hand methods for new machinery: the Laplace transform, which turns differentiation into algebra and swallows step and impulse forcings whole, and series methods for the equations where no elementary fundamental set exists at all. You are leaving the land of closed forms with the strongest general tool the elementary theory has to offer.