# Preface

This is a tutorial made solely for the purpose of education and it was designed for students taking Applied Math 0330. It is primarily for students who have very little experience or have never used Mathematica before and would like to learn more of the basics for this computer algebra system. As a friendly reminder, don't forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

# Variation of Parameters

We present the most general and powerful method for solving nonhomogeneous linear differential equations---variation of parameters method. It can be used for arbitrary driving functions in opposite, for instance, to the method of undetermined coefficients that requires a specific form of input functions and could be applied mostly for constant coefficient equations. Variation of parameters method can be used to solve arbitrary linear differential equation with integrable input functions, including piecewise continuous functions. However, practical application of the method may be limited because it needs a lot of calculations and integrations and requires explicit knowledge of a fundamental set of solutions for the associated homogeneous equation.

The method of variation of parameters was introduced by Leonhard Euler (1707--1783) and completed by his follower Joseph-Louis Lagrange (1736--1813). However, the variation of parameters method is actually an extension for higher order differential equation the Bernoulli method that is used to solve linear equations and the Bernoulli equations. In his 1749 study of the motions of the earth, Euler obtained differential equations for the orbital elements; and in 1753 he applied the method to his study of the motions of the moon. Lagrange first used the variation of parameters method in 1766 and applied it to solve some problems from celestial mechanics. It should be noted that Euler and Lagrange applied this method to nonlinear differential equations and that, instead of varying the coefficients of linear combinations of solutions to homogeneous equations, they varied the constants of the unperturbed motions of the celestial bodies. During 1808--1810, Lagrange gave the method of variation of parameters its final form in a series of papers.

Joseph-Louis Lagrange born Giuseppe Lodovico Lagrangia or Giuseppe Ludovico De la Grange Tournier (1736--1813) in Turin, Piedmont-Sardinia (now Italy). Lagrange was of Italian and French descent. Lagrange’s life divides very naturally into three periods. The first comprises the years spent in his native Turin (1736--1766). The second is that of his work at the Prussian Academy of Sciences in Berlin, between 1766 and 1787. He succeeded the director of mathematics position from Leonhard Euler (who returned to Russia and strongly recommended Joseph-Louis) and gained the full support of d'Alembert. The thirds finds him in Paris, from 1787 until his death in 1813.

The first two periods were the most fruitful in terms of scientific activity, which began as early as 1754 with the discovery of the calculus of variations and continued with the application of the latter to mechanics in 1756. He also worked in celestial mechanics in this first period, stimulated by the competitions held by the Paris Academy of Sciences in 1764 and 1766. The Berlin period was productive in mechanics as well as in differential and integral calculus. Yet during that time Lagrange distinguished himself primarily in the numerical and algebraic solution of equations, and even more in the theory of numbers. ■

We start the method with a second order nonhomogeneous linear differential equation

$L \left[ x, \texttt{D} \right] y \equiv y'' + p(x)\, y' + q(x)\, y = f(x) ,$
where $$\texttt{D} = {\text d}/ {\text d}x$$ is the derivative operator, p(x), q(x), and f(x) are given functions provided that p and q are continuous in some interval; the input function f(x) is assumed to be integrable within this interval (so it could be piecewise continuous). The variation of parameters method is applicable only when the fundamental set of solutions for associated homogeneous equation $$y'' + p(x)\, y' + q(x)\, y = 0$$ is known. Thus, we assume that y1 and y2 are known two linearly independent solutions of the homogeneous equation:
$L \left[ x, \texttt{D} \right] y_1 \equiv y''_1 + p(x)\, y'_1 + q(x)\, y_1 = 0 \qquad \mbox{and} \qquad L \left[ x, \texttt{D} \right] y_2 \equiv y''_2 + p(x)\, y'_2 + q(x)\, y_2 = 0 .$
Then we seek a particular solution of the nonhomogeneous differential equation as the sum of products:
$y(x) = A(x)\, y_1 (x) + B(x)\, y_2 (x) ,$
where A(x) and B(x) are some two smooth functions to be determined. Upon differentiating the above function, we get
$y'(x) = A'(x)\, y_1 (x) + B'(x)\, y_2 (x) + A(x)\, y'_1 + B(x)\, y'_2 .$
Since A(x) and B(x) are arbitrary functions, we impose the condition on them:
$A'(x)\, y_1 (x) + B'(x)\, y_2 (x) =0 .$
Since $$A'(x)\, y_1 (x) + B'(x)\, y_2 (x) =0 ,$$ then the first derivative $$y' (x) = A(x)\, y'_1 + B(x)\, y'_2$$ will look like a derivative of the sum $$y (x) = A\, y_1 + B\, y_2 ,$$ where the derivatives of A and B do not participate in the product rule. This allows us to find the second derivative
$y''(x) = A'(x)\, y'_1 (x) + B'(x)\, y'_2 (x) + A(x)\, y''_1 + B(x)\, y''_2 .$
Substituting the first two derivatives into the differential equation, we get
$L \left[ x, \texttt{D} \right] y = A'\,y'_1 + B'\,y'_2 + A\,y''_1 + B\,y''_2 + p(x) \left( A\,y'_1 + B\,y'_2 \right) + q(x) \left( A\,y_1 + B\,y_2 \right) =f .$
Grouping all terms associated with A(x) and all terms associated with B(x), we obtain
$A(x) \left[ y''_1 + p(x)\, y'_1 + q(x)\, y_1 \right] + B(x) \left[ y''_2 + p(x)\, y'_2 + q(x)\, y_2 \right] ,$
which are the results of the linear operator $$L \left[ x, \texttt{D} \right]$$ acting on y1 and y2, respectively. Then we rearrange terms to obtain
$L \left[ x, \texttt{D} \right] y = L \left[ x, \texttt{D} \right] y_1 + L \left[ x, \texttt{D} \right] y_2 + A'(x)\, y'_1 (x) + B'(x)\, y'_2 (x) = f(x) .$
Since y1 and y2 are two linearly independent solutions of the homogeneous equation, $$L \left[ x, \texttt{D} \right] y_1 =0$$ and $$L \left[ x, \texttt{D} \right] y_2 =0 ,$$ we get
$A'(x)\, y'_1 (x) + B'(x)\, y'_2 (x) = f(x) .$
Therefore, we have two algebraic equations with two unknowns:
$\begin{split} A'(x)\, y_1 (x) + B'(x)\, y_2 (x) &=0 , \\ A'(x)\, y'_1 (x) + B'(x)\, y'_2 (x) &= f(x) . \end{split}$
The above system of algebraic equations with respect to derivatives of unknown functions A' and B' is usually referred to as the Lagrange system. Solving the system, we get
$A'(x) = - \frac{f(x)\, y_2 (x)}{y_1 y'_2 - y'_1 y_2} = - \frac{f(x)\, y_2 (x)}{W(x)} = \frac{1}{W(x)} \, \det \begin{bmatrix} 0& y_2 \\ f(x) & y'_2 \end{bmatrix} , \qquad B'(x) = \frac{f(x)\, y_1 (x)}{y_1 y'_2 - y'_1 y_2} = \frac{f(x)\, y_1 (x)}{W(x)} = \frac{1}{W(x)} \, \det \begin{bmatrix} y_1 & 0 \\ y'_1 & f(x) \end{bmatrix} ,$
where $$W(x) = y_1 y'_2 - y_2 y'_1 = \det \begin{bmatrix} y_1 & y_2 \\ y'_1 & y'_2 \end{bmatrix}$$ is the Wronskian of two linearly independent solutions of the associated homogeneous equation $$L \left[ x, \texttt{D} \right] y =0 .$$ Integrating, we obtain
\begin{align*} A(x) &= - \int \frac{f(x)\, y_2 (x)}{W(x)}\,{\text d} x + C_1 = - \int_{x_0}^x \frac{f(\xi )\, y_2 (\xi )}{W(\xi )}\,{\text d} \xi + C_1 , \\ B(x) &= \int \frac{f(x)\, y_1 (x)}{W(x)}\,{\text d} x + C_2 = \int_{x_0}^x \frac{f(\xi )\, y_1 (\xi )}{W(\xi )} \,{\text d} \xi + C_2 , \end{align*}
where C1 and C2 are constants of integration. Here x0 is an arbitrary point from the domain where the Wronskian is not zero, and ξ is a dummy variable of integration. Being able to determine the coefficients A(x) and B(x), we construct the solution of the given nonhomogeneous differential equation
$y(x) = A(x)\, y_1 (x) + B(x)\, y_2 (x) =-y_1 (x) \int \frac{f(x)\, y_2 (x)}{W(x)}\,{\text d} x + y_2 (x) \int \frac{f(x)\, y_1 (x)}{W(x)}\,{\text d} x + C_1 y_1 (x) + C_2 y_2 (x) .$
We can unite two integrals into one:
$y(x) = \int_{x_0}^x \frac{y_2 (x )\, y_1 (\xi ) - y_1 (x )\, y_2 (\xi )}{W(\xi )} \,f(\xi ) \,{\text d} \xi + C_1 y_1 (x) + C_2 y_2 (x) = \int_{x_0}^x G(x, \xi )\, f(\xi ) \,{\text d} \xi + C_1 y_1 (x) + C_2 y_2 (x) ,$
where
$G(x, \xi ) = \frac{y_2 (x )\, y_1 (\xi ) - y_1 (x )\, y_2 (\xi )}{y_1 (\xi ) y'_2 (\xi ) - y_2 (\xi ) y'_1 (\xi )}$
is called the Green function for the differential equation $$L \left[ x, \texttt{D} \right] y = f .$$ The function is named after the British mathematical physicist George Green (1793--1841) who wrote An Essay on the Application of Mathematical Analysis to the Theories of Electricity and Magnetism (Green, 1828). The essay introduced several important concepts, among them a theorem similar to the modern Green's theorem, the idea of potential functions as currently used in physics, and the concept of what are now called Green's functions. Green's life story is remarkable in that he was almost entirely self-taught. He received only about one year of formal schooling as a child, between the ages of 8 and 9.

## Third order Differential Equations

Now we turn our attention to higher order differential equations. Let us consider the third order one:
$y''' + p(x)\, y'' + q(x) \,y' + r(x)\, y = f(x) .$
It is assumed that the fundamental set of solutions for the corresponding homogeneous equation $$y''' + p(x)\, y'' + q(x) \,y' + r(x)\, y = 0$$ is known to be
$\left\{ y_1 (x), \ y_2 (x), \ y_3 (x) \right\} .$
We seek a particular solution as a linear combination of the above linearly independent functions:
$y_p (x) = A_1 (x)\, y_1 (x) + A_2 (x)\,y_2 (x) + A_3 (x)\, y_3 (x) ,$
where the derivatives of unknown functions are determined from the Lagrange system of equations:
$\begin{split} A'_1 (x)\, y_1 (x) + A'_2 (x)\, y_2 (x) + A'_3 (x)\, y_3 (x) &=0 , \\ A'_1 (x)\, y'_1 (x) + A'_2 (x)\, y'_2 (x) + A'_3 (x)\, y'_3 (x) &=0 , \\ A'_1 (x)\, y''_1 (x) + A'_2 (x)\, y''_2 (x) + A'_3 (x)\, y''_3 (x)&= f(x) . \end{split}$
Using Cramer's rule, we find its solution:
$A'_1 (x) = \frac{1}{W(x)}\, \begin{vmatrix} 0&y_1 & y_3 \\ 0 & y'_2 & y'_3 \\ f(x) & y''_2 & y''_3 \end{vmatrix} , \quad A'_2 (x) = \frac{1}{W(x)}\, \begin{vmatrix} y_1 &0 & y_3 \\ y'_1 & 0 & y'_3 \\ y''_1 & f(x) & y''_3 \end{vmatrix} , \quad A'_3 (x) = \frac{1}{W(x)}\, \begin{vmatrix} y_1&y_2 & 0 \\ y'_1 & y'_2 & 0 \\ y''_1 & y''_2 & f(x) \end{vmatrix} ,$
where W(x) is the Wronskian of the fundamental set of solutions:
$W(x) = \det \begin{bmatrix} y_1 &y_2 & y_3 \\ y'_1 & y'_2 & y'_3 \\ y''_1 & y''_2 & y''_3 \end{bmatrix} .$

Example: Our first example concerns constant coefficient equation

$y'' + y = \tan (x).$
The associated homogeneous equation $$y'' + y =0$$ has two linearly independent solutions $$y_1 (x) = \cos x \quad \mbox{and} \quad y_2 (x) = \sin x .$$ The corresponding Wronskian is $$W(x) = y_1 y'_2 - y_2 y'_1 =1 .$$ We seek a particular solution in the form
$y(x) = A(x)\, y_1 (x) + B(x)\, y_2 (x) = A(x)\,\cos x + B(x)\,\sin x ,$
where the derivatives of coefficients A(x) and B(x) satisfy the Lagrange system of algebraic equation:

$\begin{split} A'(x)\, \cos (x) + B'(x)\, \sin (x) &=0 , \\ -A'(x)\, \sin (x) + B'(x)\, \cos (x) &= \tan (x) . \end{split}$
Its solution is
$A' (x) = - \frac{\sin^2 x}{\cos x} , \qquad B'(x) = \sin x .$
Integration yields
$A (x) = - \ln \left\vert \tan x + \sec x \right\vert + \sin x , \qquad B'(x) = -\cos x .$
This allows us to construct a particular solution
$y_p (x) = A(x)\, y_1 (x) + B(x)\, y_2 (x) = -\cos x\, \ln \left\vert \tan x + \sec x \right\vert .$
The homogeneous (or complementary) part of the solution is
$y_h (x) = C_1\, y_1 (x) + C_2\, y_2 (x) = C_1 \cos x + C_2 \sin x ,$
so the general solution becomes
$y(x) = y_p (x) + y_h (x) = -\cos x\, \ln \left\vert \tan x + \sec x \right\vert + C_1 \cos x + C_2 \sin x .$

Example: Consider nonhomogeneous differential equation with variable coefficients

$x\left( x-2 \right) y'' - \left( x^2 -2 \right) y' + 2\left( x-1 \right) y = (x-2)^2 .$
The corresponding homogeneous equation $$x\left( x-2 \right) y'' - \left( x^2 -2 \right) y' + 2\left( x-1 \right) y = 0$$ has two linearly independent solutions $$y_1 = x^2$$ and $$y_2 = e^x .$$ Therefore, we seek a particular solution of the driven equation in the form
$y_p (x) = A(x)\, y_1 + B(x)\, y_2 = A(x)\, x^2 + B(x) \, e^x ,$
where functions A(x) and B(x) are to be determined from the Lagrange system of equations:
$\begin{split} A'(x)\, x^2 + B'(x)\, e^x &=0 , \\ A'(x)\, 2x + B'(x)\, e^x &= (x-2)/x . \end{split}$
Solving the above system of algebraic equations, we find
$A' (x) = - \frac{1}{x^2} , \quad B' (x) = e^{-x} \qquad \Longrightarrow \qquad A(x) = \frac{1}{x}, \quad B(x) = - e^{-x}.$
Using these expressions, we determine a particular solution
$y_p (x) = A(x)\, y_1 + B(x)\, y_2 = x-1 .$

Example: Consider the third order differential equation

$y''' + y' = \tan x .$
The characteristic equation $$\lambda^3 + \lambda =0 ,$$ corresponding to the homogeneous equation $$y''' + y' = 0 ,$$ has three distinct roots, two of them are complex conjugate:
$\lambda_1 =0, \quad \lambda_2 = {\bf j}, \quad \lambda_3 = -{\bf j} .$
Then the complementary function (which is a solution of the homogeneous equation) becomes
$y_1 =1, \quad y_2 = \cos x, \quad y_3 = \sin x .$
Now let's try to find a particular solution to the nonhomogeneous differential equation as
$y_p (x) = A_1 (x)\,y_1 + A_2 (x)\, y_2 + A_3 (x)\, y_3 = A_1 (x) + A_2 (x)\,\cos x + A_3 (x)\, \sin x .$
The Lagrange system of algebraic equations for unknown derivatives becomes
$A'_1 (x) = \frac{1}{W(x)}\, \begin{vmatrix} 0&\cos x & \sin x \\ 0 & -\sin x & \cos x \\ \tan (x) & -\cos x & -\sin x \end{vmatrix} = \tan x, \quad A'_2 (x) = \frac{1}{W(x)}\, \begin{vmatrix} 1 &0 & \sin x \\ 0 & 0 & \cos x \\ 0 & \tan (x) & -\sin x \end{vmatrix} = - \sin x, \quad A'_3 (x) = \frac{1}{W(x)}\, \begin{vmatrix} 1&\cos x & 0 \\ 0 & -\sin x & 0 \\ 0 & -\cos x & \tan (x) \end{vmatrix} = - \frac{\sin^2 x}{\cos x} ,$
where W(x) is the Wronskian of the fundamental set of solutions:
$W(x) = \det \begin{bmatrix} 1 &\cos x & \sin x \\ 0 & -\sin x & \cos x \\ 0& -\cos x & -\sin x \end{bmatrix} = 1.$
Integration yields
$A_1 (x) = \ln |\sec x| , \quad A_2 (x) = \cos x , \quad A_3 (x) = \sin x + \ln | \sec x + \tan x | .$
Then a particular solution becomes
$y_p (x) = \ln |\sec x| + \cos x \, \ln | \sec x + \tan x | .$
Adding to it a complementary function $$y_h (x) = C_1 + C_2 \sin x + C_3 \cos x ,$$ where C1, C2, and C3 are arbitrary constants, we obtain the general solution:
$y(x) = y_p (x) + y_h (x) = \ln |\sec x| + \cos x \, \ln | \sec x + \tan x | + C_1 + C_2 \sin x + C_3 \cos x .$