Solving linear systems can sometimes be a necessary component of solving nonlinear systems.
This post is part of the book Justin Math: Linear Algebra. Suggested citation: Skycak, J. (2019). Higher-Order Variation of Parameters. In Justin Math: Linear Algebra. https://justinmath.com/higher-order-variation-of-parameters/
Want to get notified about new posts? Join the mailing list and follow on X/Twitter.
Until this point, we have been working exclusively with linear systems. However, solving linear systems can sometimes be a necessary component of solving nonlinear systems.
Second-Order Variation of Parameters
For example, recall the variation of parameters method for solving a second-order differential equation of the form
$\begin{align*} y'' + a_1(x)y' + a_0(x)y = f(x). \end{align*}$
Variation of parameters proceeds by first guessing a solution of the form
$\begin{align*} y_f(x) = u_1(x)y_1(x) + u_2(x)y_2(x), \end{align*}$
where $y_1$ and $y_2$ are the two zero solutions of the differential equation
$\begin{align*} y'^\prime+a_1(x)y'+a_0(x)y=0, \end{align*}$
and $u_1(x)$ and $u_2(x)$ are some unknown multiplier functions that we solve for by setting up a system of equations.
To set up the first equation in our system, we force
$\begin{align*} y_f'(x) = u_1(x)y_1'(x)+u_2(x)y_2'(x) \end{align*}$
and equate it to the true derivative of $y_f$:
$\begin{align*} y_f' &= u_1 y_1' + u_2 y_2' \\ (u_1y_1 + u_2y_2)' &= u_1 y_1' + u_2 y_2' \\ u_1'y_1 + u_1y_1' + u_2'y_2 + u_2y_2' &= u_1 y_1' + u_2 y_2' \\ u_1'y_1 + u_2'y_2 &= 0 \end{align*}$
The second equation comes from substituting our guess for $y_f$ into the differential equation and simplifying, using the fact that $y_1$ and $y_2$ are the zero solutions.
$\begin{align*} f &= y_f'' + a_1 y_f' + a_0 y_f \\ f &= (u_1 y_1' + u_2 y_2')' + a_1 (u_1 y_1' + u_2 y_2') + a_0 (u_1 y_1 + u_2 y_2) \\ f &= (u_1' y_1' + u_1 y_1'' + u_2' y_2' + u_2 y_2'') + a_1 (u_1 y_1' + u_2 y_2') + a_0 (u_1 y_1 + u_2 y_2) \\ f &= (u_1)(y_1''+a_1y_1'+a_0y_1) + (u_2)(y_2''+a_1y_2'+a_0y_2) + u_1'y_1' + u_2'y_2' \\ f &= (u_1)(0) + (u_2)(0) + u_1'y_1' + u_2'y_2' \\ f &= u_1'y_1' + u_2'y_2' \end{align*}$
This results in a square system of 2 equations.
$\begin{align*} \begin{cases} u_1'y_1 + u_2'y_2 = 0 \\ u_1'y_1' + u_2'y_2' = f \end{cases} \end{align*}$
In 2 dimensions, we can easily solve for $u_1’$ and $u_2’$ using elimination, obtaining the result below.
$\begin{align*} u_1' &= -\frac{y_2f}{y_1y_2'-y_2y_1'} \\ u_2' &= \frac{y_1f}{y_1y_2'-y_2y_1'} \end{align*}$
Then, we can simply integrate and substitute these back into our particular solution.
$\begin{align*} y_f &= u_1y_1 + u_2y_y \\ &= - y_1 \int \frac{y_2f}{y_1y_2'-y_2y_1'} \, dx + y_2 \int \frac{y_1f}{y_1y_2'-y_2y_1'} \, dx \end{align*}$
Higher-Order Variation of Parameters
When we wish to use variation of parameters to find the particular solution of an Nth order differential equation
$\begin{align*} y^n + a_1(x)y^{n-1} + \ldots + a_n(x)y = f(x) \end{align*}$
we guess a solution of the form
$\begin{align*} y_f(x) = u_1(x)y_{1}(x) + u_2(x)y_{2}(x) + \ldots + u_n(x)y_{n}(x) \end{align*}$
and force
$\begin{align*} \begin{matrix} y_f'(x) &=& u_1(x)y_1'(x)& +& u_2(x)y_2'(x) &+& \ldots &+& u_n(x)y_n'(x) \\ y_f''(x) &=& u_1(x)y_1''(x) &+& u_2(x)y_2''(x) &+& \ldots &+& u_n(x)y_n''(x) \\ &\vdots& & & & & & & \\ y_f^{(n-1)}(x) &=& u_1(x)y_1^{(n-1)}(x) &+& u_2(x)y_2^{(n-1)}(x) &+& \ldots &+& u_n(x)y^{(n-1)}_n(x). \end{matrix} \end{align*}$
By equating each derivative with the true derivative of $y_f$ up to order N, we can set up a system of equations.
$\begin{align*} \begin{matrix} u_1'y_1 &+& u_2'y_2 &+& \ldots &+& u_n' y_n &=& 0 \\ u_1'y_1' &+& u_2'y_2' &+& \ldots &+& u_n'y_n' &=& 0 \\ & & & & & & &\vdots& \\ u_1'y^{(n-2)}_1 &+& u_2'y^{(n-2)}_2 &+& \ldots &+& u_n'y^{(n-2)}_n &=& 0 \\ u_1'y^{(n-1)}_1 &+& u_2'y^{(n-1)}_2 &+& \ldots &+& u_n'y^{(n-1)}_n &=& f \end{matrix} \end{align*}$
This system is difficult to solve by elimination. But now we can use Cramer’s rule! First, let’s write our system more compactly, using the notation
$\begin{align*} y_i^{(0:n-1)} = \left< y_i, y_i', \ldots, y_i^{(n-1)} \right>. \end{align*}$
The system becomes
$\begin{align*} u_1' y_1^{(0:n-1)} + u_2' y_2^{(0:n-1)} + \ldots + u_n' y_n^{(0:n-1)} = \left< 0, \ldots, 0, f \right>. \end{align*}$
According to Cramer’s rule, each $u_i’$ is given by
$\begin{align*} u_i' = \frac{\mbox{det} \left( y^{(0:n-1)}_1, \ldots, y^{(0:n-1)}_{i-1}, \left< 0, \ldots, 0, f \right>, y^{(0:n-1)}_{i+1}, \ldots, y^{(0:n-1)}_n \right) }{\mbox{det} \left( y^{(0:n-1)}_1, y^{(0:n-1)}_2, \ldots, y^{(0:n-1)}_n \right) } \, . \end{align*}$
The denominator of this fraction is also known as the Wronskian, denoted
$\begin{align*} W(y_1, y_2, \ldots, y_n) = \mbox{det} \left( y^{(0:n-1)}_1, y^{(0:n-1)}_2, \ldots, y^{(0:n-1)}_n \right). \end{align*}$
If we define $W_{i,f}$ as
$\begin{align*} W_{i,f}(y_1, y_2,\ldots, y_n) = \mbox{det} \left( y^{(0:n-1)}_1, \ldots, y^{(0:n-1)}_{i-1}, \left< 0, \ldots, 0, f \right>, y^{(0:n-1)}_{i+1}, \ldots, y^{(0:n-1)}_n \right) \end{align*}$
then we have
$\begin{align*} u_i' = \frac{W_{i,f} \left( y_1, y_2, \ldots, y_n \right)}{W \left( y_1, y_2, \ldots, y_n \right)}. \end{align*}$
Finally, we can write the particular solution to the differential equation by integrating and substituting into our initial guess.
$\begin{align*} y_f = y_1 \int \frac{W_{1,f} \left( y_1, y_2, \ldots, y_n \right)}{W \left( y_1, y_2, \ldots, y_n \right)} \, dx + \ldots + y_n \int \frac{W_{n,f} \left( y_1, y_2, \ldots, y_n \right)}{W \left( y_1, y_2, \ldots, y_n \right)} \, dx \end{align*}$
Demonstration
Let’s illustrate this method on a simple example. To make it easier to find the zero solutions, we’ll choose an example with constant coefficients, but remember that this method works even when the coefficients are functions themselves.
$\begin{align*} y'''-2y''-y'+2y=e^{3x} \end{align*}$
We start off by finding the zero solutions, i.e. those that satisfy the equation whose right-hand side is zero.
$\begin{align*} y'''-2y''-y'+2y=0 \end{align*}$
We do this by finding the roots of the characteristic polynomial $p(r)=r^3-2r^2-r+2$. We can find the roots via factoring by grouping.
$\begin{align*} 0 &= r^3 - 2r^2 - r + 2 \\ 0 &= r^2(r-2)-1(r-2) \\ 0 &= (r^2-1)(r-2) \\ 0 &= (r+1)(r-1)(r-2) \\ r &= -1, 1, 2 \end{align*}$
These roots correspond to the following zero solutions:
$\begin{align*} y_{-1} = C_{-1}e^{-x} \hspace{.5cm} y_{1} = C_{1}e^{x} \hspace{.5cm} y_{2} = C_{2}e^{2x} \end{align*}$
It remains to find the particular solution. To use variation of parameters, we need three independent zero solutions, so we’ll choose the simplest ones from above: $e^{-x}, \, e^x, \, e^{2x}$.
Substituting these into the variation of parameters formula, we have a particular solution of the form
$\begin{align*} y_f &= e^{-x} \int \frac{W_{1,f} \left( e^{-x}, e^x, e^{2x} \right)}{W \left( e^{-x}, e^x, e^{2x} \right)} \, dx + e^{x} \int \frac{W_{2,f} \left( e^{-x}, e^x, e^{2x} \right)}{W \left( e^{-x}, e^x, e^{2x} \right)} \, dx \\ &\hspace{.5cm} + e^{2x} \int \frac{W_{3,f} \left( e^{-x}, e^x, e^{2x} \right)}{W \left( e^{-x}, e^x, e^{2x} \right)} \, dx \end{align*}$
with $f(x)=e^{3x}$. Now, it remains to do the computations. First, we compute the standard Wronskian in the denominator.
$\begin{align*} W \left( e^{-x}, e^x, e^{2x} \right) &= \mbox{det} \begin{pmatrix} \begin{pmatrix} e^{-x} \\ \left( e^{-x} \right)' \\ \left( e^{-x} \right) '' \end{pmatrix}, \begin{pmatrix} e^{x} \\ \left( e^{x} \right)' \\ \left( e^{x} \right) '' \end{pmatrix}, \begin{pmatrix} e^{2x} \\ \left( e^{2x} \right)' \\ \left( e^{2x} \right)'' \end{pmatrix} \end{pmatrix} \\ &= \mbox{det} \begin{pmatrix} \begin{pmatrix} e^{-x} \\ - e^{-x} \\ e^{-x} \end{pmatrix}, \begin{pmatrix} e^{x} \\ e^{x} \\ e^{x} \end{pmatrix}, \begin{pmatrix} e^{2x} \\ 2 e^{2x} \\ 4e^{2x} \end{pmatrix} \end{pmatrix} \\ &= \left( e^{-x} \right) \left( e^{x} \right) \left( e^{2x} \right) \mbox{det} \begin{pmatrix} \begin{pmatrix} 1 \\ -1 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 2 \\ 4 \end{pmatrix} \end{pmatrix} \\ &= \left( e^{-x} \right) \left( e^{x} \right) \left( e^{2x} \right) \left( 6 \right) \\ &= 6e^{2x} \end{align*}$
Next, we compute the modified Wronskians in the numerators.
$\begin{align*} W_{1,f} \left( e^{-x}, e^x, e^{2x} \right) &= \mbox{det} \begin{pmatrix} \begin{pmatrix} 0 \\ 0 \\ e^{3x} \end{pmatrix}, \begin{pmatrix} e^{x} \\ \left( e^{x} \right)' \\ \left( e^{x} \right)'' \end{pmatrix}, \begin{pmatrix} e^{2x} \\ \left( e^{2x} \right)' \\ \left( e^{2x} \right)'' \end{pmatrix} \end{pmatrix} \\ &= \mbox{det} \begin{pmatrix} \begin{pmatrix} 0 \\ 0 \\ e^{3x} \end{pmatrix}, \begin{pmatrix} e^{x} \\ e^{x} \\ e^{x} \end{pmatrix}, \begin{pmatrix} e^{2x} \\ 2 e^{2x} \\ 4e^{2x} \end{pmatrix} \end{pmatrix} \\ &= \left( e^{3x} \right) \left( e^{x} \right) \left( e^{2x} \right) \mbox{det} \begin{pmatrix} \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 2 \\ 4 \end{pmatrix} \end{pmatrix} \\ &= \left( e^{3x} \right) \left( e^{x} \right) \left( e^{2x} \right) \left( 1 \right) \\ &= e^{6x} \end{align*}$
$\begin{align*} W_{2,f} \left( e^{-x}, e^x, e^{2x} \right) &= \mbox{det} \begin{pmatrix} \begin{pmatrix} e^{-x} \\ \left( e^{-x} \right)' \\ \left( e^{-x} \right)'' \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ e^{3x} \end{pmatrix}, \begin{pmatrix} e^{2x} \\ \left( e^{2x} \right)' \\ \left( e^{2x} \right)'' \end{pmatrix} \end{pmatrix} \\ &= \mbox{det} \begin{pmatrix} \begin{pmatrix} e^{-x} \\ - e^{-x} \\ e^{-x} \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ e^{3x} \end{pmatrix}, \begin{pmatrix} e^{2x} \\ 2 e^{2x} \\ 4e^{2x} \end{pmatrix} \end{pmatrix} \\ &= \left( e^{-x} \right) \left( e^{3x} \right) \left( e^{2x} \right) \mbox{det} \begin{pmatrix} \begin{pmatrix} 1 \\ -1 \\ 1 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 2 \\ 4 \end{pmatrix} \end{pmatrix} \\ &= \left( e^{-x} \right) \left( e^{3x} \right) \left( e^{2x} \right) \left( -3 \right) \\ &= -3e^{4x} \end{align*}$
$\begin{align*} W_{3,f} \left( e^{-x}, e^x, e^{2x} \right) &= \mbox{det} \begin{pmatrix} \begin{pmatrix} e^{-x} \\ \left( e^{-x} \right)' \\ \left( e^{-x} \right)'' \end{pmatrix}, \begin{pmatrix} e^{x} \\ \left( e^{x} \right)' \\ \left( e^{x} \right)'' \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ e^{3x} \end{pmatrix} \end{pmatrix} \\ &= \mbox{det} \begin{pmatrix} \begin{pmatrix} e^{-x} \\ - e^{-x} \\ e^{-x} \end{pmatrix}, \begin{pmatrix} e^{x} \\ e^{x} \\ e^{x} \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ e^{3x} \end{pmatrix} \end{pmatrix} \\ &= \left( e^{-x} \right) \left( e^{x} \right) \left( e^{3x} \right) \mbox{det} \begin{pmatrix} \begin{pmatrix} 1 \\ -1 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \end{pmatrix} \\ &= \left( e^{-x} \right) \left( e^{x} \right) \left( e^{3x} \right) \left( 2 \right) \\ &= 2e^{3x} \end{align*}$
Lastly, we substitute back into the formula for the particular solution, and simplify.
$\begin{align*} y_f &= e^{-x} \int \frac{W_{1,f} \left( e^{-x}, e^x, e^{2x} \right)}{W \left( e^{-x}, e^x, e^{2x} \right)} \, dx + e^{x} \int \frac{W_{2,f} \left( e^{-x}, e^x, e^{2x} \right)}{W \left( e^{-x}, e^x, e^{2x} \right)} \, dx \\ &\hspace{.5cm} + e^{2x} \int \frac{W_{3,f} \left( e^{-x}, e^x, e^{2x} \right)}{W \left( e^{-x}, e^x, e^{2x} \right)} \, dx \\ &= e^{-x} \int \frac{ e^{6x} }{ 6e^{2x} } \, dx + e^{x} \int \frac{ -3e^{4x} }{ 6e^{2x} } \, dx + e^{2x} \int \frac{ 2e^{3x} }{ 6e^{2x} } \, dx \\ &= e^{-x} \int \frac{1}{6}e^{4x} \, dx + e^{x} \int -\frac{1}{2}e^{2x} \, dx + e^{2x} \int \frac{1}{3}e^{x} \, dx \\ &= e^{-x} \left( \frac{1}{24}e^{4x} \right) + e^{x} \left( -\frac{1}{4}e^{2x} \right) + e^{2x} \left( \frac{1}{3}e^{x} \right) \\ &= \frac{1}{24}e^{3x} + -\frac{1}{4}e^{3x} + \frac{1}{3}e^{3x} \\ &= \frac{1}{8}e^{3x} \end{align*}$
The full solution to the differential equation, then, is
$\begin{align*} y = C_{-1}e^{-x} + C_1e^x + C_2e^{2x} + \frac{1}{8}e^{3x} . \end{align*}$
Exercises
Solve the following differential equations using variation of parameters. (You can view the solution by clicking on the problem.)
$\begin{align*} 1) \hspace{.5cm} y'''-y''-4y'+4y=e^{2x} \end{align*}$
Solution:
$\begin{align*} y=C_{-2}e^{-2x}+C_1e^x + \left( C_2 + \frac{1}{4}x \right) e^{2x} \end{align*}$
$\begin{align*} 2) \hspace{.5cm} y'''-y''-5y'-3y=e^x \end{align*}$
Solution:
$\begin{align*} y=\left( C_{-1,0} + C_{-1,1} x \right) e^{-x} +C_3 e^{3x}- \frac{1}{8} e^x \end{align*}$
$\begin{align*} 3) \hspace{.5cm} y'''+y''-y'-y=\cos x \end{align*}$
Solution:
$\begin{align*} y=C_1 e^x + \left( C_{-1,0} + C_{-1,1} x \right)e^{-x} - \frac{1}{4} \left( \sin x + \cos x \right) \end{align*}$
$\begin{align*} 4) \hspace{.5cm} y'''-2y''-9y'+18y = \sin x \end{align*}$
Solution:
$\begin{align*} y=C_{-3}e^{-3x} + C_2 e^{2x} + C_3 e^{3x} + \frac{1}{25} \sin x + \frac{1}{50} \cos x \end{align*}$
$\begin{align*} 5) \hspace{.5cm} y'''-2y''+y'-2=\cos x \end{align*}$
Solution:
$\begin{align*} y=C_0 + \left( C_{1,0} + C_{1,1}x \right) e^x + 2x + \frac{1}{2} \cos x \end{align*}$
$\begin{align*} 6) \hspace{.5cm} y'''-y''+y'-4=\sin x \end{align*}$
Solution:
$\begin{align*} y=&\left[ C_1 \sin \left( \frac{\sqrt{15}}{2} x \right) + C_2 \cos \left( \frac{\sqrt{15}}{2} x \right) \right] e^{\frac{1}{2}x} \\ &+ x + \frac{1}{10} \sin x - \frac{3}{10} \cos x + C_3 \end{align*}$
This post is part of the book Justin Math: Linear Algebra. Suggested citation: Skycak, J. (2019). Higher-Order Variation of Parameters. In Justin Math: Linear Algebra. https://justinmath.com/higher-order-variation-of-parameters/
Want to get notified about new posts? Join the mailing list and follow on X/Twitter.