MATLAB TUTORIAL for the First Course. Part 4: Reduction Higher Order ODEs

Prof. Vladimir A. Dobrushkin

This tutorial contains many matlab scripts.
You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to

Email Vladimir Dobrushkin

Sometimes higher order differential equations can be reduced to lower and first order equations. We consider two classes of equations when this is possible.

 

Dependent Variable Missing

For a differential equation of the form \( y^{(n)} = f(x, y^{(n-1)} ) , \) the substitution \( p = y^{(n-1)} , \ p' = y^{(n)} \) leads to a first order equation of the form \( p' = f(x,p) . \) If this equation can be solved for p, then y can be obtained by integrating \( {\text d}^{n-1} y/{\text d}x^{n-1} = p . \) Note that one arbitrary constant is obtained in solving the first order equation for p, and n-1 constants are introduced in the integration for y.

Example: Consider the differential equation with dependent variable missing:

\[ x^2 \frac{{\text d}^ 2 y}{{\text d}x^2} + x\,\frac{{\text d} y}{{\text d}x} = 3x^2 +1 , \qquad x> 0. \]

Upon setting p=y', we reduce the given differential equation to a linear differential equation of first order
\[ x^2 p' + x\,p = 3x^2 +1 . \]
We seek its solution using Bernoulli method in the form: p = uv, where u is a solution of the homogeneous equation
\[ x^2 u' + x\,u = 0 \qquad \Longrightarrow \qquad \frac{{\text d}u}{u} = -2\,\frac{{\text d}x}{x} . \]
Integration yields \( u = x^{-2} . \) Then for v we get
\[ x^2 u\, v' = 3x^2 +1 \qquad \Longrightarrow \qquad v = x^3 + x + C_1 , \]
where C1 is a constant of integration. Then we solve for p and obtain
\[ y' = p = x^{-2} v = x + x^{-1} + C)1 x^{-2} \qquad \Longrightarrow \qquad y = \frac{x^2}{2} + \ln x + C_1 x^{-1} + C_2 , \]

 

Independent Variable Missing

Consider second order differential equation of the form \( y'' = f(y,y') , \) in which the independent variable t does not appear explicitly. If we let \( v = y' = \dot{y} , \) then we obtain \( {\text d} v/{\text d}t = f(y,v) . \) Since the right side of this equation depends on y and v, rather than on t and v, this equation contains too many variables. However, if we think of y as the independent variable, then by chain rule, \( {\text d} v/{\text d}t = \left( {\text d} v/{\text d}y \right) \left( {\text d} y/{\text d}t \right) = v \left( {\text d} v/{\text d}y \right) . \) Hence the original differential equation can be written as \( v \left( {\text d} v/{\text d}y \right) = f(y,v). \) Provided that this first order equation can be solved, we obtain v as a function of y. A relation between y and t results from solving \( {\text d} y/{\text d}t = v(y) , \) which is a separable equation. Again, there are two arbitrary constants in the final answer.

Example: Consider the differential equation \( y'' + y \left( y' \right)^3 =0 . \) Upon setting \( y' =v , \) we reduce the given differential equation of the second order to another one: \( v \, \frac{{\text d}v}{{\text d}y} + y \left( v \right)^3 =0 . \) Assuming that \( v \ne 0 , \) we reduce it to the first order equation

\[ \frac{{\text d}v}{{\text d}y} + y \, v^2 \qquad \Longrightarrow \qquad - \frac{{\text d}v}{v^2} = y\,{\text d}y . \]
Integration yields
\[ \frac{1}{v} = \frac{y^2}{2} + C_1 \qquad \Longrightarrow \qquad v = \frac{2}{y^2 + C_1} = \frac{{\text d}y}{{\text d}t} . \]
Next separation of variable leads to
\[ 2\,{\text d}t = \left( y^2 + C_1 \right) {\text d}y \qquad \Longrightarrow \qquad 2t + C_2 = \frac{y^3}{3} + C_1 y , \]
which is the general solution of the given differential equation in implicit form.

Example: Dr. D.H. Parsons derived in 1963

  • D.H. Parsons, Exact Integration of the Space-Charge Equation for a Plane Diode: A Simplified Theory, Nature, Volume 200, October, 126--127, 1963.
the following differential equation for the potential V in a planar vacuum diod
\[ \frac{{\text d}^2 V}{{\text d} x^2} = \frac{a}{\sqrt{V + \varepsilon^2}} \quad (0 < x<\ell ), \qquad V(0) =0, \quad V(\ell ) = V_{\ell} , \]
where \( \displaystyle a = \frac{4\pi I}{\left( 2e/m \right)^{1/2}} ; \) here m is the mass and e is the magnitude of an electron; I is the current per unit plate area.

 The above second order differential equation admits two integrating factors, one less obvious than the other:

\[ \mu (x) = - 2\, \frac{{\text d}V}{{\text d}x} \qquad \mbox{and} \qquad \mu (x) = 3\left[ \left( \frac{{\text d}V}{{\text d}x} \right)^2 - 2\,\sqrt{V= \varepsilon^2} \right] . \]
Multiplication by each integrating factor will reduce the given second order differential equation to an exact equation. Then simple integration yields
\begin{eqnarray*} 4a\,\sqrt{V + \varepsilon^2} - \left( \frac{{\text d}V}{{\text d}x} \right)^2 &=& A >0, \\ \left( \frac{{\text d}V}{{\text d}x} \right)^3 - 6a \,\frac{{\text d}V}{{\text d}x}\,\sqrt{V + \varepsilon^2} + 3a^2 a &=& B , \end{eqnarray*}
where A and B are arbitrary constants.

 In many practical problems, a second order differential equation can be reduced to an exact equation. For example, consider the equation of motion for the damped harmonic oscillator:

\[ \ddot{y} + 2\mu\,\dot{y} + \omega_0^2 y =0 , \]
where μ and ω0 are positive constants. Let
\[ \lambda_1 = \mu + \sqrt{\mu^2 - \omega_0^2} > \lambda_2 = \mu - \sqrt{\mu^2 - \omega_0^2} \]
be the roots of the characteristic equation \( \lambda^2 + 2\mu\, \lambda + \omega_0^2 =0 , \) that is, \( \lambda_1 + \lambda_2 = 2\mu \) and \( \lambda_1 \lambda_2 = \omega_0^2 . \) Multiplying both sides of the harmonic oscillator equation by \( \mu (y,t) = e^{2\mu t} \left( \dot{y} + \mu y \right) , \) we obtain an exact equation:
\[ \frac{1}{2} \, \frac{{\text d} A}{{\text d} t} =0 , \qquad\mbox{where} \quad A = e^{2\mu t} \left[ \left( \dot{y} + \mu\, y \right)^2 + \left( \omega_0^2 - \mu^2 \right) y^2 \right] . \]
Therefore, A is a constant of motion (which is also called the first integral of the equation). Note that A is similar to the energy integral for a conservatve system. The connection is easily seen if μ = 0, corresponding to the undamped oscilaltor. Following Bohlin
  • K. Bohlin, Integrationsmethode f̈ur lineare Differential-gleichungen, Astron. Iakttag. Undersokn. Stockholms Observ., Volume 9, 3--6, 1908.
it can be shown that the second order differential equation for damped oscillator admits another first integral
\[ \frac{\left( \dot{y} + \lambda_1 y \right)^{\lambda_1}}{\left( \dot{y} + \lambda_2 y \right)^{\lambda_2}} =B >0 \]
or
\[ b = \ln B = \lambda_1 \ln\left( \dot{y} + \lambda_1 y \right) - \lambda_2 \ln \left( \dot{y} + \lambda_2 y \right) . \]
This result is easily established by considering
\[ \frac{{\text d} b}{{\text d} t} = \frac{\left( \lambda_1 - \lambda_2 \right) \dot{y}}{\left( \dot{y} + \lambda_1 y \right) \left( \dot{y} + \lambda_2 y \right)} \left[ \ddot{y} + \left( \lambda_1 + \lambda_2 \right) \dot{y} + \lambda_1 \lambda_2 y \right] , \]
which vanishes by vertue of the equation of motion.

 The equation of motion for the damped harmonic oscillator may be derived from the Lagrangian

\[ {\cal L} (q, \dot{q},t) = e^{2\mu\,t} \, \frac{1}{2} \left[ \dot{q}^2 - \omega_0^2 q^2 \right] , \]
as one can easily verify. The canonical momentum associated with q is
\[ p = \frac{\partial {\cal L}}{\partial \dot{q}} = \dot{q} \, e^{2\mu\,t} . \]
The Hamiltonian is then
\[ H(q,p,t) = p\dot{q} - {\cal L} = \frac{1}{2}\, e^{-2\mu\,t} \,p^2 + \frac{\omega_0^2}{2} \,q^2 e^{2\mu\,t} . \]
Since this hamiltonian depends explicitly on time, it does not represent a conserved quantity. This is just what we should expect because we are dealing with a dissipative system.

 Let us simplify the Hamiltonian with the transformation

\[ Q = q \,e^{2\mu\,t} , \qquad P= p\,e^{-2\mu\,t} . \]
The transformed Hamiltonian is
\[ K(Q,P,t) = H(q,p,t) + \frac{\partial}{\partial t} \,q \,P\,e^{2\mu\,t} . \]
Upon some simplifications, we get
\[ K(Q,P,t) = \frac{P^2}{2} + \frac{\omega_0^2}{2} \,Q^2 + \mu\,P\,Q . \]
This formula is rather remarkable because the new hamiltonian does not depend explicitly on time, and so is conservatived. If we express the new Hamiltonian in old variables, we realize that the expression
\[ e^{-2\mu\,t} \, \frac{p^2}{2} + \frac{\omega_0^2}{2} \, q^2 e^{2\mu\,t} + \mu\,pq \]
is an integral of the motion.