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

We recall the definition of a root multiplicity. A real or complex number α is called a root of multiplicity k of the polynomial p(x) if there is a polynomial s(x) such that $$s(\alpha ) \ne 0$$ and $$p(x) = \left( x- \alpha \right)^k s(x) .$$ If k=1, then α is called a simple root. If $$k \ge 2,$$ then α is called a multiple root.

# Reduction of order for constant coefficient equations

We start with constant coefficient linear homogeneous differential equations. Suppose that in a constant coefficient linear differential operator

$L \left[ \texttt{D} \right] y =0 ,$
the operator $$L \left[ \texttt{D} \right] = a_n \texttt{D}^n + a_{n-1} \texttt{D}^{n-1} + \cdots + a_1 \texttt{D} + a_0 ,$$ with $$\texttt{D} = {\text d}/{\text d}x \quad\mbox{and} \quad a_n \ne 0,$$ has repeated factors. That is, the characteristic polynomial $$L \left( \lambda \right) = a_n \lambda^n + a_{n-1} \lambda^{n-1} + \cdots + a_0$$ has repeated roots. For example, the second order constant coefficient linear differential operator $$L \left[ \texttt{D} \right] = a\,\texttt{D}^2 + b\,\texttt{D} +c$$ has a repeated factor if and only if the corresponding characteristic equation $$a\,\lambda^2 + b\,\lambda + c =0$$ has a double root:
$\lambda_1 = \lambda_2 = -b/(2a) .$
In other words, the quadratic polynomial can be factored $$a\,\lambda^2 + b\,\lambda + c = a \left( \lambda - \lambda_1 \right)^2$$ if and only if its discriminant $$b^2 - 4ac$$ is zero. In this case we have only one solution of exponential form:
$y_1 (x) = e^{-bx/(2a)}$
for the differential equation
$a\, y'' + b\, y' + c\,y =0, \qquad b^2 = 4ac.$

To find another linearly independent solution to the above homogeneous equation, we use the method of reduction of order, credited to Jacob Bernoulli (1655--1705). Setting

$y(x) = v(x)\, y_1 (x) = v(x)\, e^{-bx/(2a)} ,$
we have
\begin{align*} y' &= v' (x)\, y_1 (x) + v(x)\, y'_1 (x) = \left[ v' (x) - \frac{b}{2a}\, v(x) \right] y_1 (x) , \\ y'' &= \left[ v'' (x) - \frac{b}{a}\, v' (x) + \left( \frac{b}{2a} \right)^2 \right] y_1 (x) . \end{align*}
Upon substituting the Bernoulli form into the differential equation, we get
\begin{align*} L \left[ \texttt{D} \right] y &= a \left[ v'' (x) - \frac{b}{a}\, v' (x) + \left( \frac{b}{2a} \right)^2 \right] y_1 (x) + b \left[ v' (x) - \frac{b}{2a}\, v(x) \right] y_1 (x) + c\, v(x)\, y_1 (x) \\ &= a\, v'' (x)\, y_1 (x) + v(x) \,L \left[ \texttt{D} \right] y_1 (x) = a\, v'' (x)\, y_1 (x) =0. \end{align*}
The latter can be divided by nonzero term $$a\,y_1 (x) = a\,e^{-bx/(2a)} .$$ This yields $$v'' (x) =0,$$ and after integration, we obtain $$v (x) = C_1 + C_2 x,$$ with arbitrary constants C1 and C2. Finally, substituting for v(x), we define the general solution
$y(x) = v(x)\, y_1 (x) = \left( C_1 + C_2 x \right) e^{-bx/(2a)} .$
Functions $$y_1 (x) = e^{-bx/(2a)} \quad \mbox{and} \quad y_2 (x) = x\,e^{-bx/(2a)}$$ are linearly independent because their Wronskian
$W(x) = \begin{vmatrix} e^{\gamma x} & x\, e^{\gamma x} \\ \gamma \, e^{\gamma x} & \left( x\,\gamma + 1 \right) e^{\gamma x} \end{vmatrix} = e^{2\gamma x} \ne 0 .$
Therefore, these two functions form a fundamental set of solutions whenever γ is. ■

This method can be extended for arbitrary case when a linear constant coefficient differential operator has a multiple $$\left( \texttt{D} - \gamma \right)^m$$ of multiplicity m. Then to such multiple correspond m linearly independent solutions:

$y_1 (x) = e^{-x\gamma} , \ y_2 (x) = x\, e^{-x\gamma} , \ \ldots , \ y_m (x) = x^{m-1} e^{-x\gamma} .$
The strict derivation of the above statement is made in annihilation section.

Theorem: Let $$a_0 , a_1 , \ldots , a_n$$ be n+1 real (or complex) numbers with $$a_n \ne 0 ,$$ and y(x) be a n times continuously differentiable function on some interval |a,b|. Then y(x) is a solution of the n-th order linear differential equation with constant coefficients

$L \left[ \texttt{D} \right] y \equiv a_n \texttt{D}^n y + a_{n-1} \texttt{D}^{n-1} y + \cdots + a_1 \texttt{D} \, y+ a_0 \, y =0 , \qquad \texttt{D} = {\text d}/{\text d} x ,$
if and only if
$y (x) = \sum_{j=1}^m e^{-x\lambda_j x} \, P_j (x) ,$
where $$\lambda_1 , \lambda_2 , \ldots , \lambda_m$$ are distinct roots of the characteristic polynomial $$\sum_{k=1}^n a_k \lambda^k$$ with multiplicities $$m_1 , m_2 , \ldots , m_m ,$$ respectively, and $$P_k (x)$$ is a polynomial of degree $$m_k -1 .$$

We consider a partcular case of second order differential equation---the motion for undriven damped harmonic oscillator (coefficients μ and ω0 > 0 are assumed constants):

$\ddot{y} + 2\mu\,\dot{y} + \omega_0^2 y =0 ,$
subject to the initial conditions
$y(0) = y_0 , \qquad \dot{y}(0) = v_0 .$
We show that its solution can be obtained from the solution
\begin{align*} y(t) &= \frac{y_0 \lambda_2 - v_0}{\lambda_2 - \lambda_1}\, e^{\lambda_1 t} - \frac{y_0 \lambda_1 - v_0}{\lambda_2 - \lambda_1}\, e^{\lambda_2 t} \\ &= y_0 \, e^{\lambda_1 t} + \left( y_0 \lambda_1 - v_0 \right) \left[ \frac{e^{\lambda_1 t} - e^{\lambda_2 t}}{\lambda_2 - \lambda_1} \right] \end{align*}
by taking the limit when $$\lambda_1 \ \to \, \lambda_2 .$$ Here
$\lambda_1 , \ \lambda_2 = - \mu \pm \eta = -\mu \pm \left( \mu^2 - \omega_0^2 \right)^{1/2} , \qquad \eta = \left( \mu^2 - \omega_0^2 \right)^{1/2}$
are two distinct roots of the characteristic equation $$\lambda^2 + 2\mu\,\lambda + \omega_0^2 =0 .$$ Let assume, for simplicity, that $$\eta = \left( \mu^2 - \omega_0^2 \right)^{1/2} \ge 0 ,$$ that is $$\mu \ge \omega_0 .$$ Then rewriting the solution in the equivalent form:
$y(t) = y_0 e^{-\mu\, t} \,\cosh (\eta \,t) + \left( \mu y_0 + v_0 \right) e^{-\mu\,t} \,\frac{\sinh (\eta\, t)}{\eta} , \qquad \eta = \left( \mu^2 - \omega_0^2 \right)^{1/2} \ge 0,$
and taking the limit as $$\eta \, \to \, 0,$$ one gets
$y(t) = e^{-\mu\, t} \left[ y_0 + \left( \mu y_0 + v_0 \right) t \right] .$

Example: Let us consider the differential equation

$y'' -4\,y' +4\,y =0 .$
The characteristic equation $$\lambda^2 - 4\,\lambda +4 = \left( \lambda -2 \right)^2 =0$$ has a double root $$\lambda =2 .$$ Hence we get one exponential solution $$y_1 (x) = e^{2x}$$ and another one is a multiple of the latter: $$y_2 (x) = x\,e^{2x} .$$ So the general solution becomes
$y(x) = C_1 e^{2x} + C_2 x\,e^{2x} ,$
with some arbitrary constants C1 and C2.
DSolve[y''[x]-4 y'[x]+4 y[x]==0,y[x],x]
soln[x_] = Expand[y[x]/.%[[1]]]
y1[x_] = Coefficient[soln[x],C[1]]
y2[x_] = Coefficient[soln[x], C[2]]
basis = {y1[x], y2[x]}
WronskianDet = Det[{basis, D[basis, x]}]
Out[22]= E^(4 x)

Example: Consider the initial value problem

$y'' +6\,y' +9\,y =0 , \qquad y(0) = -1, \quad y' =(0) = 5.$
The characteristic polynomial $$\lambda^2 + 6\,\lambda +9 = = \left( \lambda +3 \right)^2$$ has one double root $$\lambda =-3 .$$ Therefore, the general solution of the given homogeneous differential equation becomes
$y(x) = C_1 e^{-3x} + C_2 x\, e^{-3x} .$
To satisfy the initial conditions, we have to choose arbitrary constants that are solutions of the system of equations:
$\begin{split} y(0) = C_1 =-1, \\ y'(0) = -3\, C_1 + C_2 =5. \end{split}$
So $$C_1 =-1, \quad C_2 =2$$ and the required solution becomes
$y(x) = - e^{-3x} + 2\, x\, e^{-3x} .$

soln = DSolve[{x''[t] + 6 x'[t] + 9 x[t] == 0, x[0] == -1,
x'[0] == 5}, x[t], t]
s2[t_] = x[t] /. soln[[1]]
Plot[s2[t], {t, 0.2, 3}]

Out[11]= {{x[t] -> E^(-3 t) (-1 + 2 t)}}
Out[12]= E^(-3 t) (-1 + 2 t)
Out[13]=

When does the extreme occur?
Solve[s2'[t] == 0, t]
out[48]= {{t -> 5/6}}

Where is the extreme?
s2[t /. %]
Out[49]= {2/(3 E^(5/2))}
N[%]
Out[50]= {0.0547233}

Example: Let us consider the fourth order linear differential operator:

$L \left[ \texttt{D} \right] y \equiv \texttt{D}^4 y - 8\, \texttt{D}^{3} y + 18\, \texttt{D}^2 \, y -27 \, y =0 , \qquad \texttt{D} = {\text d}/{\text d} x .$
The corresponding characteristic polynomial
$L \left( \lambda \right) = \lambda^4 - 8\, \lambda^{3} + 18\, \lambda^2 -27 = \left( \lambda -3 \right)^3 \left( \lambda +1 \right)$
has one simple root $$\lambda =-1 \0 and one triple root \( \lambda =3 .$$ Therefore, the general solution to the fourth order differential equation $$L \left[ \texttt{D} \right] y =0$$ is
$y(x) = C_1 e^{-x} + \left( C_2 + C_3 x + C_4 x^2 \right) e^{3x} ,$
where Ck, $$k=1,2,3,4 ,$$ are arbitrary constants.

# Reduction of order for variable coefficient equations

The reduction of order technique, which applies to arbitrary linear differential equations, allows us to go beyond equations with constant coefficients, provided that we already know one solution. For sake of clarity, we start with a second order linear differential equation with variable coefficients:

$y'' +p(x)\,y' +q(x)\,y =0 ,$
where p(x) and q(x) are some continuous functions on some interval |a,b|. Suppose that we know one its solution $$y = y_1 (x) \ne 0.$$ This means that
$L \left[ \texttt{D} \right] y_1 =0 , \qquad\mbox{where} \quad L \left[ \texttt{D} \right] = \texttt{D}^2 + p(x)\,\texttt{D} + q(x) .$
Then we seek another linearly independent solution in the form $$y(x) = v(x)\, y_1 (x) ,$$ where unknown function v(x) s determined upon substitution of this form into the given differential equation. First, we use the product rule to obtain
\begin{align*} y' &= v' (x)\, y_1 (x) + v(x)\, y'_1 (x) , \\ y'' &= v'' (x)\, y_1 (x) +2\, v' (x)\, y_1 (x) + v(x)\, y''_1 (x) . \end{align*}
This implies that we must have
\begin{align*} L \left[ \texttt{D} \right] v(x)\,y_1 &= v'' (x)\, y_1 (x) +2\, v' (x)\, y'_1 (x) + v(x)\, y''_1 (x) + p(x) \left[ v' (x)\, y_1 (x) + v(x)\, y'_1 (x) \right] + q(x)\,v(x)\,y_1 (x) \\ &= v'' (x)\, y_1 (x) +2\, v' (x)\, y'_1 (x) + p(x)\,v' (x)\, y_1 (x) + v(x)\, L \left[ \texttt{D} \right] y_1 =0. \end{align*}
Since $$L \left[ \texttt{D} \right] y_1 \equiv 0 ,$$ the function v(x) must satisfy the equation
$v'' (x)\, y_1 (x) + v' (x) \left[ 2\,y'_1 (x) + p(x) \, y_1 (x) \right] =0 .$
This is a second order differential equation where dependent function is missing. If we set $$u= v' ,$$ we reduce it to a first order differential equation:
$u' (x)\, y_1 (x) = - u (x) \left[ 2\,y'_1 (x) + p(x) \, y_1 (x) \right] \qquad\Longrightarrow \qquad \frac{{\text d} u}{u} = - \frac{1}{y_1 (x)} \left[ 2\,y'_1 (x) + p(x) \, y_1 (x) \right] {\text d} x .$
Integrating both sides, we get
$\ln u = - 2\int \frac{y'_1 (x)}{y_1 (x)} \,{\text d}x - \int p(x) \, {\text d} x = \ln y_1^{-2} - \int p(x) \, {\text d} x .$
This allows us to determine u(x) explicitly:
$v' = u = y_1^{-2} \, \exp \left\{ - \int p(x) \, {\text d} x \right\} ,$
where we dropped an arbitrary constant because we are after another linearly independent solution. Next integration yields
$v = \int y_1^{-2} \, \exp \left\{ - \int p(x) \, {\text d} x \right\} {\text d} x \qquad\Longrightarrow \qquad y_2 (x) = y_1 (x) \, \int y_1^{-2} \, \exp \left\{ - \int p(x) \, {\text d} x \right\} {\text d} x .$

Example: Find a second linearly independent solution to the differential equation

$t\,y'' - y' + 4t^3\,y =0 , \qquad\mbox{for} \quad t>0,$
provided that $$y_1 (t) = \sin \left( t^2 \right)$$ is a solution. First, let's rewrite our differential equation by diving by t to get:
$y'' - \frac{1}{t}\,y' + 4t^2\,y =0 .$
We seek another solution in the form $$y= v(t)\, y_1 (t) = v(t)\, \sin \left( t^2 \right) .$$ Then for u=v'(t) we have the separable equation
$\frac{u'}{u} = - 2\,\frac{y'_1}{y_1} + \frac{1}{t} \qquad\Longrightarrow \qquad \ln u = \ln \left( t\,y_1^{-2} \right) .$
Then
$v' = u = \frac{t}{y_1} \qquad\Longrightarrow \qquad u = \int \frac{t}{\sin^2 t^2} \, {\text d}t = \frac{1}{2}\, \cot \left( t^2 \right) .$
So another linearly independent solution becomes
$y_2 (t) = v(t)\,y_1 (t) = \frac{1}{2}\,\cos \left( t^2 \right) .$
Of course, the constant multiple (1/2) could be dropped.