Particular Solutions by Variation of Parameters

In this section we learn how to find a particular solution of any linear differential equation, provided we already know the general solution of the corresponding homogeneous equation. Unlike undetermined coefficients, this method is not limited to constant coefficients and places no restriction on the forcing function. Its cost is that integrals must be evaluated rather than coefficients matched.

Quick Reference

Item Statement
Equation (order 2) y'' + P(x)y' + Q(x)y = R(x), in standard form
Ansatz y p = u 1 ( x ) y 1 ( x ) + u 2 ( x ) y 2 ( x )
Auxiliary condition u_1'y_1 + u_2'y_2 = 0
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}
Cramer's rule u_1' = \dfrac{-y_2R(x)}{W}, \qquad u_2' = \dfrac{y_1R(x)}{W}
Closed formula y p = y 1 ( x ) y 2 ( x ) R ( x ) W ( x ) d x + y 2 ( x ) y 1 ( x ) R ( x ) W ( x ) d x
Order n u_i' = (-1)^{n-i}\dfrac{W_i\,f(x)}{W},   1 i n
W i Wronskian of the set obtained by deleting y i from { y 1 , , y n }
Constants Omit the constants of integration; only one particular solution is needed

The Second-Order Theorem

If y 1 ( x ) and y 2 ( x ) are two linearly independent solutions of the homogeneous equation

L[y] = \frac{d^{2}y}{dx^{2}} + P(x)\frac{dy}{dx} + Q(x)y = 0 , \tag{i}

then we can always find a particular solution of the nonhomogeneous equation

\frac{d^{2}y}{dx^{2}} + P(x)\frac{dy}{dx} + Q(x)y = R(x) \tag{ii}

by the method of variation of parameters. The method always works provided that the general solution of the corresponding homogeneous equation is already known.

If P ( x ) , Q ( x ) , and R ( x ) are continuous on an open interval I , and if y 1 ( x ) and y 2 ( x ) are two linearly independent solutions of

y'' + P(x)y' + Q(x)y = 0 ,

then a particular solution of

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

is

y p = y 1 ( x ) y 2 ( x ) R ( x ) W ( x ) d x + y 2 ( x ) y 1 ( x ) R ( x ) W ( x ) d x

where W ( x ) = W [ y 1 , y 2 ] ( x ) is the Wronskian.

  • Since we need only one particular solution, we can omit the constants of integration in the above integrals.
  • Sometimes the expression for y p is written with definite integrals, y p = y 1 ( x ) x 0 x y 2 ( t ) R ( t ) W ( t ) d t + y 2 ( x ) x 0 x y 1 ( t ) R ( t ) W ( t ) d t , where x 0 is an arbitrarily chosen point in I . This form has the pleasant feature that y_p(x_0) = y_p'(x_0) = 0.

Derivation of the Method

Full derivation

If y 1 ( x ) and y 2 ( x ) are two linearly independent solutions of (i), the general solution of the homogeneous equation is

y h ( x ) = c 1 y 1 ( x ) + c 2 y 2 ( x ) ,

where c 1 and c 2 are two parameters. In the method of variation of parameters we replace those constants by two unknown functions u 1 ( x ) and u 2 ( x ) , and try to choose them so that

y_p = u_1(x)y_1(x) + u_2(x)y_2(x) \tag{iii}

becomes a particular solution. Letting the parameters vary is what gives the method its name.

There are two unknown functions, so we need two equations to determine them. One equation is the requirement that y p satisfy (ii). We are free to choose the second, and the choice below is what makes the computation tractable.

Differentiating (iii),

y_p' = u_1y_1' + u_2y_2' + \left(u_1'y_1 + u_2'y_2\right) .

If we assume that

u_1'y_1 + u_2'y_2 = 0, \qquad x \in I \tag{iv}

then the parenthesized group vanishes and

\begin{aligned} y_p' &= u_1y_1' + u_2y_2' \\ y_p'' &= u_1'y_1' + u_1y_1'' + u_2'y_2' + u_2y_2'' . \end{aligned}

The point of condition (iv) is now visible: without it, y_p'' would contain u_1'' and u_2'', and the system would be a pair of second-order equations rather than a pair of algebraic ones.

Substituting y p , y_p', and y_p'' into (ii) and grouping,

\begin{aligned} u_1&\left[y_1'' + P(x)y_1' + Q(x)y_1\right] + u_2\left[y_2'' + P(x)y_2' + Q(x)y_2\right] \\ &\quad + u_1'y_1' + u_2'y_2' = R(x) \end{aligned}

Since y 1 and y 2 are solutions of the homogeneous equation, both square brackets equal zero, and the equation reduces to

u_1'y_1' + u_2'y_2' = R(x) . \tag{v}

Equations (iv) and (v) form a system of two linear equations for the two unknown derivatives u_1' and u_2':

\left\{ \begin{aligned} &y_1u_1' + y_2u_2' = 0 \\[6pt] &y_1'u_1' + y_2'u_2' = R(x) \end{aligned} \right.

or in matrix form

\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} \tag{vi}

Notice that since y 1 and y 2 are linearly independent solutions, the determinant of the coefficient matrix, which is their Wronskian, is never zero, so this system always has a unique solution. Everything hinges on that fact, and it is exactly what Section 3.3 established.

By Cramer's rule,

u_1' = \frac{\begin{vmatrix} 0 & y_2 \\ R(x) & y_2' \end{vmatrix}}{\begin{vmatrix} y_1 & y_2 \\ y_1' & y_2' \end{vmatrix}} = \frac{-y_2R(x)}{W\left(y_1,y_2\right)}u_2' = \frac{\begin{vmatrix} y_1 & 0 \\ y_1' & R(x) \end{vmatrix}}{\begin{vmatrix} y_1 & y_2 \\ y_1' & y_2' \end{vmatrix}} = \frac{y_1R(x)}{W\left(y_1,y_2\right)}

We find u 1 ( x ) and u 2 ( x ) by integrating these results.

Standard form is mandatory. The derivation assumed the coefficient of y'' is 1 , so R ( x ) is the right-hand side after dividing by the leading coefficient. For x^{2}y'' - 2xy' + 2y = x^{3}\ln x, the correct R is x ln x , not x 3 ln x .

Worked Examples

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

Solution

Homogeneous solutions. The characteristic equation r 2 1 = 0 gives y 1 = e x and y 2 = e x , with

W = | e x e x e x e x | = 1 1 = 2 .

The derivatives. With R ( x ) = 2 1 + e x ,

u_1' = \frac{-y_2R}{W} = \frac{-e^{-x}\cdot\dfrac{2}{1+e^{x}}}{-2} = \frac{e^{-x}}{1+e^{x}} ,u_2' = \frac{y_1R}{W} = \frac{e^{x}\cdot\dfrac{2}{1+e^{x}}}{-2} = -\frac{e^{x}}{1+e^{x}} .

Integrate. The second is immediate:

u 2 = ln ( 1 + e x ) .

For the first, substitute t = e x , so d x = d t t :

u 1 = 1 / t 1 + t d t t = d t t 2 ( 1 + t ) = ( 1 t 2 1 t + 1 1 + t ) d t = 1 t ln t + ln ( 1 + t ) ,

that is,

u 1 = e x x + ln ( 1 + e x ) .

Assemble.

\begin{aligned} y_p &= u_1e^{x} + u_2e^{-x} \\ &= e^{x}\left[-e^{-x} - x + \ln\left(1+e^{x}\right)\right] - e^{-x}\ln\left(1+e^{x}\right) \\ &= -1 - xe^{x} + \left(e^{x} - e^{-x}\right)\ln\left(1+e^{x}\right) \end{aligned}

The general solution is

y = c 1 e x + c 2 e x 1 x e x + ( e x e x ) ln ( 1 + e x ) .

Solve x^{2}y'' - 2xy' + 2y = x^{3}\ln x for x > 0 . Note that the coefficients are not constant, so undetermined coefficients is unavailable.

Solution

Homogeneous solutions. The homogeneous equation is of Cauchy-Euler type (Section 4.5). Trying y = x r gives

r ( r 1 ) 2 r + 2 = r 2 3 r + 2 = ( r 1 ) ( r 2 ) = 0 ,

so y 1 = x and y 2 = x 2 , with

W = | x x 2 1 2 x | = 2 x 2 x 2 = x 2 .

Standard form. Divide the full equation by x 2 :

y'' - \frac{2}{x}y' + \frac{2}{x^{2}}y = x\ln x ,

so R ( x ) = x ln x , not x 3 ln x .

The derivatives.

u_1' = \frac{-x^{2}\cdot x\ln x}{x^{2}} = -x\ln x, \qquad u_2' = \frac{x\cdot x\ln x}{x^{2}} = \ln x .

Integrate, both by parts:

u 1 = x ln x d x = ( x 2 2 ln x x 2 4 ) = x 2 2 ln x + x 2 4 u 2 = ln x d x = x ln x x

Assemble.

\begin{aligned} y_p &= u_1x + u_2x^{2} \\ &= -\frac{x^{3}}{2}\ln x + \frac{x^{3}}{4} + x^{3}\ln x - x^{3} \\ &= \frac{x^{3}}{2}\ln x - \frac{3x^{3}}{4} \end{aligned}

The general solution is

y = c 1 x + c 2 x 2 + x 3 2 ln x 3 4 x 3 .

Check. With y p = 1 2 x 3 ln x 3 4 x 3 ,

y_p' = \frac{3x^{2}}{2}\ln x + \frac{x^{2}}{2} - \frac{9x^{2}}{4} = \frac{3x^{2}}{2}\ln x - \frac{7x^{2}}{4} ,y_p'' = 3x\ln x + \frac{3x}{2} - \frac{7x}{2} = 3x\ln x - 2x .

Substituting,

x 2 ( 3 x ln x 2 x ) 2 x ( 3 x 2 2 ln x 7 x 2 4 ) + 2 ( x 3 2 ln x 3 x 3 4 ) = ( 3 3 + 1 ) x 3 ln x + ( 2 + 7 2 3 2 ) x 3 = x 3 ln x .

Variation of Parameters for Third Order and Higher

The application of variation of parameters is not limited to second-order equations. We applied the same idea to a first-order linear equation in Section 2.6, where there was only one unknown function and no auxiliary condition was required.

If the nonhomogeneous linear differential equation is of order n > 2 ,

y ( n ) + P 1 ( x ) y ( n 1 ) + + P n ( x ) y = f ( x )

and y 1 ( x ) , y 2 ( x ) , , y n ( x ) are n linearly independent solutions of the corresponding homogeneous equation, then a particular solution is obtained by setting

y p = u 1 ( x ) y 1 ( x ) + u 2 ( x ) y 2 ( x ) + + u n ( x ) y n ( x )

and finding the derivatives of the unknown functions by solving

\begin{bmatrix} y_1 & y_2 & \cdots & y_n \\ y_1' & y_2' & \cdots & y_n' \\ \vdots & \vdots & \ddots & \vdots \\ y_1^{(n-2)} & y_2^{(n-2)} & \cdots & y_n^{(n-2)} \\ y_1^{(n-1)} & y_2^{(n-1)} & \cdots & y_n^{(n-1)} \end{bmatrix} \begin{bmatrix} u_1' \\ u_2' \\ \vdots \\ u_{n-1}' \\ u_n' \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ f(x) \end{bmatrix}

The n 1 zeros are the higher order analogues of condition (iv): each one kills the derivatives of the u k at one stage of the differentiation. Notice that the determinant of the coefficient matrix is W ( y 1 , y 2 , , y n ) , which is nonzero by independence.

Solving by Cramer's rule yields

u_i' = (-1)^{n-i}\,\frac{W_i\,f(x)}{W}, \qquad 1 \leq i \leq n

where W is the Wronskian of { y 1 , , y n } and W i is the Wronskian of the set obtained by deleting y i from { y 1 , , y n } .

The third-order case written out

For

y''' + P_1(x)y'' + P_2(x)y' + P_3(x)y = f(x)

with y 1 , y 2 , y 3 three linearly independent solutions of the homogeneous equation, a particular solution is

y p = u 1 ( x ) y 1 ( x ) + u 2 ( x ) y 2 ( x ) + u 3 ( x ) y 3 ( x )

where

u_1' = \frac{W_1f(x)}{W}, \qquad u_2' = -\frac{W_2f(x)}{W}, \qquad u_3' = \frac{W_3f(x)}{W}

with

W = \begin{vmatrix} y_1 & y_2 & y_3 \\ y_1' & y_2' & y_3' \\ y_1'' & y_2'' & y_3'' \end{vmatrix}

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} .

Note the alternating signs, which come from the cofactor expansion and are already built into the factor ( 1 ) n i of the general formula. A fully worked third-order example appears in Section 3.4.

Exercises

Solve y'' + y = \csc x on 0 < x < π .

Solution

Here y 1 = cos x , y 2 = sin x , and W = 1 .

u_1' = -\sin x\csc x = -1 \quad \Longrightarrow \quad u_1 = -xu_2' = \cos x\csc x = \cot x \quad \Longrightarrow \quad u_2 = \ln|\sin x|

Therefore

y = c 1 cos x + c 2 sin x x cos x + sin x ln ( sin x ) .

Solve y'' - 2y' + y = \dfrac{e^{x}}{x} for x > 0 .

Solution

The double root r = 1 gives y 1 = e x , y 2 = x e x , and W = e 2 x as computed in Section 3.4.

u_1' = -\frac{xe^{x}\cdot e^{x}/x}{e^{2x}} = -1 \quad \Longrightarrow \quad u_1 = -xu_2' = \frac{e^{x}\cdot e^{x}/x}{e^{2x}} = \frac{1}{x} \quad \Longrightarrow \quad u_2 = \ln x

Therefore

y = c 1 e x + c 2 x e x x e x + x e x ln x .

The term x e x is itself a homogeneous solution and may be absorbed into c 2 x e x .

Solve x^{2}y'' + xy' - y = x^{2} for x > 0 , given the homogeneous solutions y 1 = x and y 2 = 1 x from Section 4.1.

Solution W = | x 1 / x 1 1 / x 2 | = 1 x 1 x = 2 x .

In standard form, R ( x ) = x 2 x 2 = 1 . So

u_1' = \frac{-\left(1/x\right)\cdot 1}{-2/x} = \frac{1}{2} \quad \Longrightarrow \quad u_1 = \frac{x}{2}u_2' = \frac{x\cdot 1}{-2/x} = -\frac{x^{2}}{2} \quad \Longrightarrow \quad u_2 = -\frac{x^{3}}{6}

Therefore

y p = x 2 x x 3 6 1 x = x 2 2 x 2 6 = x 2 3

and

y = c 1 x + c 2 x + x 2 3 .

Check: x 2 2 3 + x 2 x 3 x 2 3 = 2 x 2 3 + 2 x 2 3 x 2 3 = x 2 . ✓

Frequently Asked Questions

When should I use variation of parameters?

Whenever undetermined coefficients does not apply: variable coefficients, or a forcing function like sec x , ln x , 1 x , or e x 1 + x 2 . It also works in every case where undetermined coefficients works, so it is the safe fallback, at the cost of computing integrals.


Why impose the condition u_1'y_1 + u_2'y_2 = 0?

Two unknown functions give more freedom than one equation constrains, so one extra condition can be imposed at will. This particular choice makes u_1'' and u_2'' disappear from y_p'', turning what would be a system of second-order equations into a simple algebraic system. Any other choice would be legal but far less convenient.


Do I need constants of integration for u 1 and u 2 ?

No. Adding a constant k 1 to u 1 adds k 1 y 1 to y p , which is a homogeneous solution and gets absorbed into y h . Since any particular solution will do, take the simplest antiderivative each time.


What is R ( x ) exactly?

The right-hand side of the equation after it has been divided by the leading coefficient. For 3y'' + \cdots = 6x, take R ( x ) = 2 x . For x^{2}y'' + \cdots = x^{3}\ln x, take R ( x ) = x ln x . Skipping this division is the most frequent error.


What if I cannot evaluate the integrals?

Leave them. Writing y p in terms of definite integrals from x 0 to x is a complete and rigorous answer, suitable for numerical evaluation or for theoretical estimates. The method itself never fails; only the antidifferentiation may.


Does the method require constant coefficients?

No, and this is its main advantage. All it requires is a known fundamental set for the homogeneous equation. The Cauchy-Euler example above has variable coefficients and yields to the method without difficulty.


Why does the general formula have the factor ( 1 ) n i ?

It comes from expanding the Cramer's rule determinant along the last column, which is ( 0 , , 0 , f ) T . The cofactor of the entry in row n , column i , carries the sign ( 1 ) n + i , which equals ( 1 ) n i since the two exponents differ by the even number 2 i . For n = 2 this reproduces the minus sign in u_1' and the plus sign in u_2'.