Every solution of a nonhomogeneous linear equation $L [y] = f$ has the form $y = y_{h} + y_{p}$ : the general solution of the associated homogeneous equation plus any one solution of the full equation. Variation of parameters is the method that supplies that one particular solution, and unlike the method of undetermined coefficients, it works for every forcing function $f$ .

Quick Reference

Item	Statement
Structure theorem	The general solution of $L [y] = f$ is $y = y_{h} + y_{p}$
$y_{h}$	General solution of $L [y] = 0$ , containing $n$ arbitrary constants
$y_{p}$	Any single solution of $L [y] = f$ , with no arbitrary constants
Second-order ansatz	$y_{p} = u_{1} (x) y_{1} (x) + u_{2} (x) y_{2} (x)$
Defining system	\begin{bmatrix} y_1 & y_2 \\ y_1' & y_2' \end{bmatrix}\begin{bmatrix} u_1' \\ u_2' \end{bmatrix} = \begin{bmatrix} 0 \\ R(x) \end{bmatrix}
Second-order formulas	u_1' = -\dfrac{R(x)y_2}{W}, \qquad u_2' = \dfrac{R(x)y_1}{W}
$n$ th order formula	u_k' = \dfrac{W_k(x)}{W(x)}R(x), with $W_{k}$ the Wronskian determinant having column $k$ replaced by $(0, \dots, 0, 1)^{T}$
Requirement	The equation must be in standard form, so the coefficient of $y^{(n)}$ is $1$

The Structure of the General Solution

Let $y_{p}$ be any particular solution of $L [y] = f (x)$ on an interval $I$ , and let $y_{h} = c_{1} ϕ_{1} + \dots + c_{n} ϕ_{n}$ be the general solution of $L [y] = 0$ . Then every solution of $L [y] = f$ on $I$ has the form

y = y_{h} + y_{p} .

Proof

First, any such $y$ is a solution, since by linearity

L [y_{h} + y_{p}] = L [y_{h}] + L [y_{p}] = 0 + f = f .

Conversely, let $ψ$ be any solution of $L [y] = f$ . Then

L [ψ - y_{p}] = L [ψ] - L [y_{p}] = f - f = 0,

so $ψ - y_{p}$ solves the homogeneous equation. By the fundamental set theorem of Section 3.3, it equals $c_{1} ϕ_{1} + \dots + c_{n} ϕ_{n}$ for some constants, that is, $ψ - y_{p} = y_{h}$ . Therefore $ψ = y_{h} + y_{p}$ .

The consequence is a clean division of labor. Finding all solutions of the nonhomogeneous equation requires finding all solutions of the homogeneous one, plus just one solution of the full equation. The word "any" in the theorem is important: it does not matter which particular solution you find, since two different choices differ by a homogeneous solution and get absorbed into $y_{h}$ .

Variation of Parameters

Suppose we know a fundamental set ${y_{1}, y_{2}}$ for the second-order equation

y'' + P(x)y' + Q(x)y = R(x) .

The homogeneous general solution is $c_{1} y_{1} + c_{2} y_{2}$ with constant $c_{1}$ and $c_{2}$ . The idea of variation of parameters, met already in Section 2.6, is to let those constants vary:

y_{p} = u_{1} (x) y_{1} (x) + u_{2} (x) y_{2} (x) .

Substituting this into the equation gives one condition on two unknown functions, so we are free to impose a second condition of our own choosing. Choosing u_1'y_1 + u_2'y_2 = 0 keeps all second derivatives of $u_{1}$ and $u_{2}$ out of the calculation, and the two conditions together become the system

\begin{bmatrix} y_1 & y_2 \\ y_1' & y_2' \end{bmatrix} \begin{bmatrix} u_1' \\ u_2' \end{bmatrix} = \begin{bmatrix} 0 \\ R(x) \end{bmatrix}

The coefficient matrix is the Wronskian matrix, whose determinant $W (y_{1}, y_{2})$ is nonzero because the solutions are independent. Solving by Cramer's rule,

u_1' = -\frac{R(x)\,y_2}{W\left(y_1,y_2\right)}, \qquad u_2' = \frac{R(x)\,y_1}{W\left(y_1,y_2\right)}

Integrating each gives $u_{1}$ and $u_{2}$ , and assembling them gives $y_{p}$ . The derivation is carried out in full in Chapter 4; what matters here is the pattern and the fact that it never fails, provided the integrals can be evaluated.

Standard form is mandatory. The formulas assume the coefficient of y'' is $1$ . If the equation reads a_0(x)y'' + \cdots = b(x), divide through first, so that $R (x) = b (x) / a_{0} (x)$ . Forgetting this step is the most common source of wrong answers.

Worked Examples

Find the general solution of

y'' - 2y' + y = \frac{e^{x}}{1+x^{2}} .

Note that this equation cannot be solved by the method of undetermined coefficients.

Solution

The reason undetermined coefficients fails is that $\frac{e^{x}}{1 + x^{2}}$ is not of the special form (polynomial times exponential times sine or cosine) that the method requires. Variation of parameters has no such restriction.

Step 1: the homogeneous solutions. The characteristic equation is

r^{2} - 2 r + 1 = (r - 1)^{2} = 0 .

Because this has the double root $r_{1, 2} = 1$ , a fundamental set of solutions is ${e^{x}, x e^{x}}$ :

y_{1} = e^{x} and y_{2} = x e^{x} .

Step 2: set up the system. To find a particular solution, write
$y_{p} = u_{1} (x) e^{x} + u_{2} (x) x e^{x} .$ The derivatives of $u_{1}$ and $u_{2}$ are determined by $[\begin{matrix} y_{1} & y_{2} y_{1}^{'} & y_{2}^{'} \end{matrix}] [\begin{matrix} u_{1}^{'} u_{2}^{'} \end{matrix}] < / h 1 >< h 1 > [\begin{matrix} 0 R (x) \end{matrix}] ⟹ [\begin{matrix} e^{x} & x e^{x} e^{x} & e^{x} + x e^{x} \end{matrix}] [\begin{matrix} u_{1}^{'} u_{2}^{'} \end{matrix}] < / h 1 > [\begin{matrix} 0 \frac{e^{x}}{1 + x^{2}} \end{matrix}] $ whose solution is $ u_{1}^{'} = - \frac{R (x) y_{2}}{W (y_{1}, y_{2})}, u_{2}^{'} = \frac{R (x) y_{1}}{W (y_{1}, y_{2})} . $ Step 3: the Wronskian. $ W (y_{1}, y_{2}) = | \begin{matrix} e^{x} & x e^{x} e^{x} & e^{x} + x e^{x} \end{matrix} | = e^{2 x} + x e^{2 x} - x e^{2 x} = e^{2 x} . $ Step 4: compute u_1' and u_2'. $ u_{1}^{'} = - \frac{1}{e^{2 x}} \cdot \frac{e^{x}}{1 + x^{2}} \cdot x e^{x} = - \frac{x}{1 + x^{2}} $ and $ u_{2}^{'} = \frac{1}{e^{2 x}} \cdot \frac{e^{x}}{1 + x^{2}} \cdot e^{x} = \frac{1}{1 + x^{2}} . $ The exponentials cancel completely, leaving two elementary integrals. Step 5: integrate.
\begin{aligned} u_{1}(x) &= \int -\frac{x}{1+x^{2}},dx \ &= -\int\frac{dt}{2t} \qquad \left(t = 1+x^{2},\ dt = 2x,dx\right) \ &= -\frac{1}{2}\ln|t| = -\frac{1}{2}\ln\left(1+x^{2}\right) [6pt] u_{2}(x) &= \int\frac{dx}{1+x^{2}} = \tan^{-1}x . \end{aligned}
The constants of integration are omitted, since any particular solution will do. Step 6: assemble.
\begin{aligned} y_{p} &= u_{1}(x)y_{1}(x) + u_{2}(x)y_{2}(x) \ &= -\frac{1}{2}e^{x}\ln\left(1+x^{2}\right) + xe^{x}\tan^{-1}x \end{aligned}
Therefore the general solution of the given equation is
\begin{aligned} y &= y_{h} + y_{p} = c_{1}y_{1} + c_{2}y_{2} + y_{p} \ &= c_{1}e^{x} + c_{2}xe^{x} - \frac{1}{2}e^{x}\ln\left(1+x^{2}\right) + xe^{x}\tan^{-1}x . \end{aligned}

Variation of Parameters for Higher Order Equations

The method extends to any order. For an $n$ th order equation in standard form with fundamental set ${y_{1}, \dots, y_{n}}$ , write $y_{p} = u_{1} y_{1} + \dots + u_{n} y_{n}$ and impose $n - 1$ auxiliary conditions of the same kind as before. Cramer's rule then gives

u_k'(x) = \frac{W_k(x)}{W(x)}\,R(x)

where $W$ is the full Wronskian and $W_{k}$ is the determinant obtained by replacing the $k$ th column of the Wronskian matrix with ${(0, 0, \dots, 1)}^{T}$ . For $n = 3$ this works out to

u_1' = \frac{W_1}{W}R, \qquad u_2' = -\frac{W_2}{W}R, \qquad u_3' = \frac{W_3}{W}R

where the alternating signs come from the cofactor expansion, and

W_1 = \begin{vmatrix} y_2 & y_3 \\ y_2' & y_3' \end{vmatrix}, \qquad W_2 = \begin{vmatrix} y_1 & y_3 \\ y_1' & y_3' \end{vmatrix}, \qquad W_3 = \begin{vmatrix} y_1 & y_2 \\ y_1' & y_2' \end{vmatrix} .

Find the general solution of

y''' - 4y'' + 3y' = \frac{e^{2x}}{1+e^{x}} .

Solution

Step 1: the characteristic equation.

\begin{aligned} r^{3} - 4r^{2} + 3r &= r\left(r^{2} - 4r + 3\right) = r(r-1)(r-3) = 0 \\ &\Longrightarrow \quad r = 0 \ \text{ or } \ r = 1 \ \text{ or } \ r = 3 \end{aligned}

Therefore three linearly independent solutions of the corresponding homogeneous equation are

y_{1} = 1, y_{2} = e^{x}, y_{3} = e^{3 x} .

Note that the root $r = 0$ contributes the constant solution $y = e^{0 \cdot x} = 1$ .

Step 2: the Wronskians.

W (y_{1}, y_{2}, y_{3}) = | \begin{matrix} 1 & e^{x} & e^{3 x} \\ 0 & e^{x} & 3 e^{3 x} \\ 0 & e^{x} & 9 e^{3 x} \end{matrix} | = | \begin{matrix} e^{x} & 3 e^{3 x} \\ e^{x} & 9 e^{3 x} \end{matrix} | = 9 e^{4 x} - 3 e^{4 x} = 6 e^{4 x}

where the first column, having a single nonzero entry, was expanded along.

W_{1} = | \begin{matrix} e^{x} & e^{3 x} \\ e^{x} & 3 e^{3 x} \end{matrix} | = 3 e^{4 x} - e^{4 x} = 2 e^{4 x}

W_{2} = | \begin{matrix} 1 & e^{3 x} \\ 0 & 3 e^{3 x} \end{matrix} | = 3 e^{3 x} W_{3} = | \begin{matrix} 1 & e^{x} \\ 0 & e^{x} \end{matrix} | = e^{x}

Step 3: the derivatives u_k'.

u_{1}' = \frac{W_{1}}{W}R(x) = \frac{2e^{4x}}{6e^{4x}}\cdot\frac{e^{2x}}{1+e^{x}} = \frac{e^{2x}}{3\left(1+e^{x}\right)}u_{2}' = -\frac{W_{2}}{W}R(x) = -\frac{3e^{3x}}{6e^{4x}}\cdot\frac{e^{2x}}{1+e^{x}} = -\frac{1}{2}\cdot\frac{e^{2x}}{e^{x}+e^{2x}} = -\frac{1}{2}\cdot\frac{e^{x}}{1+e^{x}}u_{3}' = \frac{W_{3}}{W}R(x) = \frac{e^{x}}{6e^{4x}}\cdot\frac{e^{2x}}{1+e^{x}} = \frac{e^{2x}}{6\left(e^{3x}+e^{4x}\right)} = \frac{1}{6\left(e^{x}+e^{2x}\right)}

Step 4: integrate. Every integral yields to the substitution $t = e^{x}$ , for which $d t = e^{x} d x$ , that is, $d x = \frac{d t}{t}$ .

\begin{aligned} u_{1} &= \frac{1}{3}\int\frac{e^{2x}}{1+e^{x}}\,dx = \frac{1}{3}\int\frac{t^{2}}{1+t}\cdot\frac{dt}{t} = \frac{1}{3}\int\frac{t}{1+t}\,dt \\ &= \frac{1}{3}\int\frac{t+1-1}{1+t}\,dt = \frac{1}{3}\int\left(1 - \frac{1}{1+t}\right)dt \\ &= \frac{1}{3}\left(t - \ln|1+t|\right) = \frac{1}{3}\left(e^{x} - \ln\left(1+e^{x}\right)\right) \end{aligned}\begin{aligned} u_{2} &= -\frac{1}{2}\int\frac{e^{x}}{1+e^{x}}\,dx = -\frac{1}{2}\int\frac{t}{1+t}\cdot\frac{dt}{t} = -\frac{1}{2}\int\frac{dt}{1+t} \\ &= -\frac{1}{2}\ln|1+t| = -\frac{1}{2}\ln\left(1+e^{x}\right) \end{aligned}\begin{aligned} u_{3} &= \frac{1}{6}\int\frac{dx}{e^{x}+e^{2x}} = \frac{1}{6}\int\frac{1}{t+t^{2}}\cdot\frac{dt}{t} = \frac{1}{6}\int\frac{dt}{t^{2}(1+t)} \\ &= \frac{1}{6}\int\left(\frac{1}{t^{2}} - \frac{1}{t} + \frac{1}{1+t}\right)dt \qquad \text{(partial fractions)} \\ &= \frac{1}{6}\left(-\frac{1}{t} - \ln|t| + \ln|1+t|\right) \\ &= \frac{1}{6}\left[-e^{-x} - x + \ln\left(1+e^{x}\right)\right] \end{aligned}

The partial fraction decomposition used above is

\frac{1}{t^{2} (1 + t)} = \frac{1}{t^{2}} - \frac{1}{t} + \frac{1}{1 + t},

which you can verify by combining the right side over the common denominator $t^{2} (1 + t)$ .

Step 5: assemble.

\begin{aligned} y_{p} &= u_{1}y_{1} + u_{2}y_{2} + u_{3}y_{3} \\ &= \frac{1}{3}\left[e^{x} - \ln\left(1+e^{x}\right)\right] - \frac{1}{2}e^{x}\ln\left(1+e^{x}\right) + \frac{1}{6}e^{3x}\left[-e^{-x} - x + \ln\left(1+e^{x}\right)\right] \end{aligned}

and the general solution is

y = c_{1} + c_{2} e^{x} + c_{3} e^{3 x} + y_{p} .

Note that the term $\frac{1}{3} e^{x}$ inside $y_{p}$ is itself a homogeneous solution, so it can be absorbed into $c_{2} e^{x}$ if you prefer a tidier answer. Extra homogeneous pieces appearing inside $y_{p}$ are harmless and are exactly what the "any particular solution" clause of the structure theorem permits.

Exercises

Use variation of parameters to solve y'' + y = \sec x on $- \frac{π}{2} < x < \frac{π}{2}$ .

Solution

The homogeneous solutions are $y_{1} = \cos x$ and $y_{2} = \sin x$ , with $W = \cos^{2} x + \sin^{2} x = 1$ .

u_1' = -\frac{\sec x\sin x}{1} = -\tan x \quad \Longrightarrow \quad u_1 = \ln|\cos x|u_2' = \frac{\sec x\cos x}{1} = 1 \quad \Longrightarrow \quad u_2 = x

Therefore

y_{p} = \cos x \ln | \cos x | + x \sin x

and the general solution is

y = c_{1} \cos x + c_{2} \sin x + \cos x \ln | \cos x | + x \sin x .

Undetermined coefficients could not have produced this, since $\sec x$ is not of the required form.

Solve y'' - 3y' + 2y = \dfrac{1}{1 + e^{-x}}.

Solution

The characteristic equation $r^{2} - 3 r + 2 = (r - 1) (r - 2)$ gives $y_{1} = e^{x}$ , $y_{2} = e^{2 x}$ , and

W = | \begin{matrix} e^{x} & e^{2 x} \\ e^{x} & 2 e^{2 x} \end{matrix} | = e^{3 x} .

Writing $R (x) = \frac{1}{1 + e^{- x}} = \frac{e^{x}}{1 + e^{x}}$ ,

u_1' = -\frac{R\,e^{2x}}{e^{3x}} = -\frac{R}{e^{x}} = -\frac{1}{1+e^{x}}, \qquad u_2' = \frac{R\,e^{x}}{e^{3x}} = \frac{1}{e^{x}\left(1+e^{x}\right)} .

With $t = e^{x}$ and $d x = d t / t$ ,

u_{1} = - \int \frac{d t}{t (1 + t)} = - \int (\frac{1}{t} - \frac{1}{1 + t}) d t = - \ln t + \ln (1 + t) = \ln \frac{1 + e^{x}}{e^{x}},

u_{2} = \int \frac{d t}{t^{2} (1 + t)} = - \frac{1}{t} - \ln t + \ln (1 + t) = - e^{- x} - x + \ln (1 + e^{x}) .

Therefore

y = c_{1} e^{x} + c_{2} e^{2 x} + e^{x} \ln (1 + e^{- x}) + e^{2 x} [- e^{- x} - x + \ln (1 + e^{x})] .

Suppose $y_{p_{1}}$ solves $L [y] = f_{1}$ and $y_{p_{2}}$ solves $L [y] = f_{2}$ . Write down a particular solution of $L [y] = 3 f_{1} - 5 f_{2}$ .

Solution

By linearity, $L [3 y_{p_{1}} - 5 y_{p_{2}}] = 3 L [y_{p_{1}}] - 5 L [y_{p_{2}}] = 3 f_{1} - 5 f_{2}$ , so

y_{p} = 3 y_{p_{1}} - 5 y_{p_{2}} .

This is the superposition principle for nonhomogeneous equations, and it lets you split a complicated forcing function into simple pieces and handle each separately.

Frequently Asked Questions

Why is the general solution

y_{h} + y_{p}

Because the difference of any two solutions of $L [y] = f$ solves $L [y] = 0$ . So once you have one solution $y_{p}$ , every other solution differs from it by a homogeneous solution, and the homogeneous solutions are exactly $y_{h}$ . Nothing more and nothing less.

Does it matter which particular solution I use?

No. Two different particular solutions differ by a homogeneous solution, which is absorbed into $y_{h}$ by relabeling the arbitrary constants. This is why you can safely drop the constants of integration when computing $u_{1}$ and $u_{2}$ .

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

Use undetermined coefficients when $f (x)$ is a polynomial, an exponential, a sine or cosine, or a product of these, because guessing is faster. Use variation of parameters for everything else: $\sec x$ , $\ln x$ , $\frac{e^{x}}{1 + x^{2}}$ , or any forcing function whose derivatives keep producing new forms. Variation of parameters always works; undetermined coefficients does not.

Why must the equation be in standard form?

Because the formulas take $R (x)$ to be the right-hand side of an equation whose leading coefficient is $1$ . If you apply them to 2y'' + \cdots = g(x) without dividing first, every $u_{k}$ comes out twice too large. Divide by the leading coefficient before you begin.

What if the integrals for

u_{1}

and

u_{2}

cannot be evaluated?

Leave them as integrals. Writing $u_{1} (x) = - \int \frac{R (t) y_{2} (t)}{W (t)} d t$ is a legitimate and complete answer, and it can be evaluated numerically or bounded analytically. The method itself never fails; only the final antidifferentiation might.

Why does the extra condition u_1'y_1 + u_2'y_2 = 0 not lose solutions?

Because you are looking for one particular solution, not all of them. Two unknown functions offer more freedom than one equation constrains, and the extra condition simply picks a convenient point in that freedom. It is chosen precisely to eliminate u_1'' and u_2'' from the computation, which is what keeps the resulting system algebraic rather than differential.

The General Solution and Variation of Parameters

Quick Reference

The Structure of the General Solution

Variation of Parameters

Worked Examples

Variation of Parameters for Higher Order Equations

Exercises

Frequently Asked Questions

Quick Links

Social Media