Converting to a system

Given a single ordinary differential equation, one method of finding numerical solutions entails transferring it into an equivalent system of differential equations of the first order. When a single differential equation has an isolated highest derivative, it is always possible to transfer the differential equation into an equivalent system of differential equations of the first order. This can be achieved by denoting sequential derivatives of the unknown variable with a new dependent variable with the exception of the last derivative, which is used to incorporate the given single differential equation. Other options are also possible for such a conversion, which will be clear from the examples.

The word "equivalent" indicates that both a single differential equation and the corresponding system of equations have the same set of solutions (so no solution is lost or added). The opposite is not always true, and some systems of n ordinary differential equations cannot be reduced to an equivalent single n-th order differential equation. The general solution of the single n-th order differential equation depends on n arbitrary constants and it belongs to the space of functions having n continuous derivatives, which is usually denoted by Cⁿ. A system of n differential equations has the same properties and its general solution also depends on n arbitrary constants. Therefore, a single n-th order ODE can be equivalent to the system of (n+1) first order ODEs only when the latter does not require extra differentiation of the unknown function. However, a single n-th order ODE cannot be equivalent to the (n−1) system of first order ODEs because the latter contains (n−1) arbitrary constants and requires utilization of C^n-1. A key illustration highlighting this is the pendulum equation, presented in Example 5, which does not use the derivative operator, and where a single second order equation is transformed into four first order equations representing an equivalent system. On the other hand, the derivative operator is used in Example 1, showing that the resulting system of equations is not equivalent to the original equation. Therefore, the key takeaway is to avoid using the derivative operator because it is an unbounded operation.

There are a couple of reasons why it is convenient to transfer a single differential equation of an order higher than one to a system of differential equations of first order. First of all, it is simpler to analyze theoretically the first order vector differential equation than deal with a higher order differential equation because the latter can be naturally included into a broader topic: systems of differential equations. Second, it is easier to use and implement a universal numerical algorithm for a first order derivative operator instead of utilizing special numerical procedures for higher order derivatives. It is more accurate to monitor errors of numerical calculations for a system of first order differential equations than for a single equation. Also, there are many situations in which we not only want to know the solution to an ODE (acronym for Ordinary Differential Equations), but also the derivative or acceleration of the solution. This information can be naturally extracted from the solution to the system of first order differential equations (with no extra work).

When treating the independent variable as a time, it is common to use Newton's notation \( \dot{y} \) for the derivative instead of Lagrange notation y' or Leibniz's notation \( \displaystyle {\text d}y/{\text d}t . \) Sometimes, we also utilize Euler's notation for the derivative operator: \( \displaystyle \texttt{D} = {\text d}/{\text d}t . \)

Example 1: Let us start with a simple case of constant coefficient linear differential equation when all manipulations become transparent. So we consider the following initial value problem for a linear second order differential equation

\begin{equation*} %\label{EqConv.1} x''(t) + 2x'(t) - 3x(t) = 0, \qquad x(0) =3, \quad x' (0) =-1. \tag{1.1} \end{equation*}

This is a homogeneous, linear, second-order ODE with constant coefficients. As a brief refresher from the previous course, recall that its solution can be computed by finding the roots of the corresponding characteristic equation:

\[ \lambda^2 + 2\,\lambda - 3 = 0. \tag{1.2} \]

There are two distinct, real roots to the quadratic (characteristic) equation: λ₁ = 1 and λ₂ = -3, and the general solution is a linear combination of two exponent functions:

\[ x(t) = c_{1}e^{\lambda_{1}t} + c_{2}e^{\lambda_{2}t} = c_{1}e^{t} + c_{2}e^{-3t} , \]

where c₁ and c₂ are arbitrary constants. To satisfy the initial conditions, we have to choose these constants so that x(0) = c₁ + c₂ = 3 and x'(0) = c₁ -3 c₂ = -1. This yields c₁ = 2 and c₂ = 1; and the solution becomes

\begin{equation*} %\label{EqConv.2} x(t) = 2\, e^t + e^{-3t} \qquad \Longrightarrow \qquad x'(t) = 2\, e^t -3\, e^{-3t} . \tag{1.3} \end{equation*}

DSolve[{x''[t] + 2*x'[t] - 3*x[t] == 0, x[0] == 3, x'[0] == -1}, x[t], t]

Now we convert the given problem to a vector equation by letting x₁ = x and its first derivative we denote by x₂. Note that second derivative is known from the given equation \( x'' = 3\,x(t) - 2\, x' (t) = 3\, x_1 - 2\, x_2 . \) This leads to the system of first order differential equations

\[ \begin{split} x'_1 &= x_2 = x', \\ x'_2 &= x'' = -2\, x_2 + 3\,x_1 . \end{split} \]

Upon introducing a column vector x = [ x₁, x₂]^T (here "T" stands for transposition; see first chapter of this tutorial), we can rewrite the system of ODEs in a more appropriate vector form:

\begin{equation*} %\label{EqConv.3} \frac{\text d}{{\text d}t} \begin{bmatrix} x_1 (t) \\ x_2 (t) \end{bmatrix} = \begin{bmatrix} 0& \phantom{-}1 \\ 3 & -2 \end{bmatrix} \begin{bmatrix} x_1 (t) \\ x_2 (t) \end{bmatrix} , \qquad \begin{bmatrix} x_1 (0) \\ x_2 (0) \end{bmatrix} = \begin{bmatrix} \phantom{-}3 \\ -1 \end{bmatrix} . \tag{1.4} \end{equation*}

DSolve[{x1'[t] == x2[t], x2'[t] == -2*x2[t] + 3*x1[t], x1[0] == 3, x2[0] == -1}, {x1[t], x2[t]}, t]

{{x1[t] -> E^(-3 t) (1 + 2 E^(4 t)), x2[t] -> E^(-3 t) (-3 + 2 E^(4 t))}}

With matrix notation,

\[ {\bf x} (t) = \begin{bmatrix} x_1 (t) \\ x_2 (t) \end{bmatrix} , \qquad {\bf A} = \begin{bmatrix} 0& \phantom{-}1 \\ 3 & -2 \end{bmatrix} , \qquad {\bf x}_0 = \begin{bmatrix} \phantom{-}3 \\ -1 \end{bmatrix} , \]

we rewrite the system of ODEs and corresponding initial conditions in compact form:

\[ \frac{\text d}{{\text d}t}\,{\bf x}(t) = {\bf A}\,{\bf x}(t) , \qquad {\bf x}(0) = {\bf x}_0 . \]

In order to find the solution to this vector differential equation, we first need to find the eigenvalues of the corresponding matrix (which is called the companion matrix):

\[ \det \left( \lambda {\bf I} - {\bf A} \right) =0 \qquad \Longleftrightarrow \qquad \det \left( \lambda \begin{bmatrix} 1&0 \\ 0&1 \end{bmatrix} - \begin{bmatrix} 0&\phantom{-}1 \\ 3&-2 \end{bmatrix} \right) = \det \begin{bmatrix} \lambda & -1 \\ -3& \lambda +2 \end{bmatrix} = \lambda^2 + 2\,\lambda - 3 =0 , \]

where I is the identity matrix. From observation, the eigenvalues of matrix A are the same as the roots of the characteristic equation λ² + 2 λ −3 = 0 for the given second order differential equation.

Sometimes, we need to monitor the residual of numerical calculations by introducing an additional auxiliary variable y₃ = x'', the second derivative of the unknown solution x(t), and keep other variables: y₁ = x and y₂ = x'. Then its derivative, according to the given differential equation, will be:

\[ y'_3 = x''' = \frac{\text d}{{\text d}t} \left( x'' \right) = \frac{\text d}{{\text d}t} \left( -2\,y_2 + 3\, y_1 \right) = -2\,y'_2 + 3\, y'_1 = -2 \left( -2\, y_2 + 3\, y_1 \right) + 3\, y_2 = 7\, y_2 -6\, y_1 . \]

This allows us to compose the system of three equations. There are two options to achieve it. First, we recall that the derivative of y₂ is our new dependent variable y₃. So we get

\begin{equation*} %\begin{split} \begin{cases} y'_1 &= y_2 = x' , \\ y'_2 &= y_3 = x'' , \\ y'_3 &= -2\, y_3 +3 \, y_2 = x''' ; \end{cases} %\end{split} \qquad\quad %\begin{split} \begin{cases} y_1 (0) &= 3 \\ y_2 (0) &= -1 \\ y_3 (0) &= 11 . \end{cases} %\end{split} \tag{1.5} \end{equation*}

Therefore, we obtain a vector differential equation \( \dot{\bf y} = {\bf C}\,{\bf y}(t) \) that is equivalent (we hope?) to the the original IVP (1.1)

\[ \frac{\text d}{{\text d}t} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} = \begin{bmatrix} 0&1&\phantom{-}0 \\ 0&0&\phantom{-}1 \\ 0&3&-2 \end{bmatrix} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} , \qquad \begin{bmatrix} y_1 (0) \\ y_2 (0) \\ y_3 (0) \end{bmatrix} =\begin{bmatrix} \phantom{1}3 \\ -1 \\ 11 \end{bmatrix} , \qquad\mbox{with}\quad {\bf C} = \begin{bmatrix} 0&1&\phantom{-}0 \\ 0&0&\phantom{-}1 \\ 0&3&-2 \end{bmatrix} . \]

The initial condition for this function also follows from the given differential equation (1.1):

\[ y_3 = x'' = 3\, x(t) - 2\, x' (t) \qquad \Longrightarrow \qquad y_3 (0) = 3\,x(0) -2\, x' (0) = 3\cdot 3 -2 \cdot (-1) = 11 . \]

Note that the derivative of y₃ is expressed through the original function x(t) as \( y'_3 = 7\, y_2 -6\, y_1 . \) So another option is to rewrite the given system of ODEs as

\begin{equation*} %\label{EqConv.4} %\begin{split} \begin{cases} y'_1 &= y_2 = x' , \\ y'_2 &= y_3 = x'' , \\ y'_3 &= -6\, y_1 +7\, y_2 = x''' ; \end{cases} %\end{split} \qquad\quad %\begin{split} \begin{cases} y_1 (0) &= 3 \\ y_2 (0) &= -1 \\ y_3 (0) &= 11 . \end{cases} %\end{split} \tag{1.6} \end{equation*}

In vector form, it becomes:

\[ \frac{\text d}{{\text d}t} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} = \begin{bmatrix} \phantom{-}0&1&0 \\ \phantom{-}0&0&1 \\ -6&7&0 \end{bmatrix} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} , \qquad \begin{bmatrix} y_1 (0) \\ y_2 (0) \\ y_3 (0) \end{bmatrix} =\begin{bmatrix} \phantom{1}3 \\ -1 \\ 11 \end{bmatrix} , \qquad\mbox{with}\quad {\bf B} = \begin{bmatrix} \phantom{-}0&1&0 \\ \phantom{-}0&0&1 \\ -6&7&0 \end{bmatrix} . \]

DSolve[{y1'[t] == y2[t], y2'[t] == -2*y2[t] + 3*y1[t], y3'[t] == 7*y2[t] - 6*y1[t], y1[0] == 3, y2[0] == -1, y3[0] == 11}, {y1[t], y2[t], y3[t]}, t]

{{y1[t] -> E^(-3 t) (1 + 2 E^(4 t)), y2[t] -> E^(-3 t) (-3 + 2 E^(4 t)), y3[t] -> E^(-3 t) (9 + 2 E^(4 t))}}

Despite that the first two components (y₁ and y₂) of the solution to the initial value problem (1.6) are the same as (1.3),

\[ \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \end{bmatrix} = \begin{bmatrix} e^{-3t} + 2\,e^{t} \\ -3\,e^{-3t} + 2\,e^{t} \\ 9\,e^{-3t} + 2\,e^{t} \end{bmatrix} , \]

the system of diferential equations (1.6) is not equivalent to (1.1) because the corresponding matrix B has three eigenvalues {−3, 2, 1}. Recall that matrix A has two eigenvalues {−3, 1}. It becomes clear when we rewrite the system of three equations with matrix B in the following form with Euler's notation \( \texttt{D} = {\text d}/{\text d}t \) for derivative operator:

\[ \begin{split} \texttt{D} y_1 &= y_2 , \\ \texttt{D} y_2 &= y_3 , \\ \texttt{D} y_3 &= -6\,y_1 + 7\,y_2 . \end{split} \]

The system above contains three unbounded derivative operators, but the original differential equation (1.1) is of order two. So its general solution is

\[ y_1 (t) = x(t) = C_1 e^t + C_2 e^{-3t} + C_3 e^{2t} , \]

with some arbitrary constants C₁, C₂, and C₃. Upon using the appropriate initial conditions we eliminated unwanted additional constant C₃ from its general solution. So we obtain the same solution from the system of three equations as the original IVP, at least from theoretical point of view. However, matrix B has a dominant eigenvalue λ = 2>1 and the corresponding numerical algorithm should avoid it.

Then the difference of y₃ and \( 3\, x(t) - 2\,x' (t) = 3\,y_1 - 2\, y_2 \) gives the residual (error) of calculations.

We use the standard Mathematica command NDSolve to solve numerically the system of three equations (1.6) and plot the residual of calculations for the equation. Blue curve gives the error with the true second derivative, and the red curve gives the residual based on numerical calculation.

sol3 = NDSolve[{y1'[t] == y2[t], y2'[t] == y3[t], y3'[t] == -6*y1[t] + 7*y2[t], y1[0] == 3, y2[0] == -1, y3[0] == 11}, {y1, y2, y3}, {t, 0, 3}];
sol2 = NDSolve[{x1'[t] == x2[t], x2'[t] == 3*x1[t] - 2*x2[t], x1[0] == 3, x2[0] == -1}, {x1, x2}, {t, 0, 3}];
a1 = Plot[ Evaluate[y3[t] /. sol3] + Evaluate[(2*x2[t] - 3*x1[t]) /. sol2], {t, 0, 1}, PlotTheme -> "Web"]
a2 = Plot[9 E^(-3 t) + 2 E^t - Evaluate[(y3[t]) /. sol3], {t, 0, 1}, PlotTheme -> "Business"]

Figure 1. Residuals of the second derivative: the error between the second derivative y3(t) and the true value in blue.		Figure 2. The residual of the second derivative y3(t) = d²x/dt² and the numerical solution 3x2(t) - 2x1(t).

Note that introducing the system of three equations with artificial differentiation does not provide an equivalent system for the given equation because the derivative operator is unbounded.

The transformation to a system of first order ODEs is straight forward, but it is not unique. Let us try another set of new dependent variables:

\[ z_1 = x, \quad z_2 = \dot{x} - x \qquad\Longrightarrow \qquad \dot{z}_2 = \ddot{x} - \dot{x} = -3 \left( \dot{x} - x \right) . \]

This leads to the system of first order differential equations:

\[ \begin{split} \dot{z}_1 &= z_1 + z_2 \qquad \Longleftrightarrow \qquad \dot{x} = \dot{x} , \\ \dot{z}_2 &= -3\, z_2 . \end{split} \]

We can also rewrite the above system in vector form

\[ \frac{\text d}{{\text d}t} \begin{bmatrix} z_1 (t) \\ z_2 (t) \end{bmatrix} = \begin{bmatrix} 1& \phantom{-}1 \\ 0 & -3 \end{bmatrix} \begin{bmatrix} z_1 (t) \\ z_2 (t) \end{bmatrix} , \qquad \begin{bmatrix} z_1 (0) \\ z_2 (0) \end{bmatrix} = \begin{bmatrix} \phantom{-}3 \\ -4 \end{bmatrix} , \]

because \( z_2 (0) = \dot{x}(0) - x(0) = -1 -3 =-4 . \) The characteristic polynomial for the corresponding matrix is

\[ \det \left( \lambda \begin{bmatrix} 1&0 \\ 0&1 \end{bmatrix} - \begin{bmatrix} 1&\phantom{-}1 \\ 0&-3 \end{bmatrix} \right) = \det \begin{bmatrix} \lambda -1& -1 \\ 0& \lambda +3 \end{bmatrix} = \lambda^2 + 2\,\lambda - 3 =0 . \]

■

Example 2: With our next example, we show that variable coefficients do not prevent conversion of a single equation of the second order into a system of two equations of the first order. The differential equation for modeling non-linear electrical oscillations has the form

\[ \texttt{D}^2 x(t) + q(t)\, x(t) = 0 , \tag{2.1} \]

where \( \ddot{x} = \frac{{\text d}^2 x}{{\text d} t^2} = \texttt{D}^2 x(t) \) stands for the second derivative with respect to time variable t, and q(t) is a single-valued periodic function. This linear differential equation has practical importance in the analysis of time-variable circuits. If function q(t) is a periodic function, then Eq.(2.1) is called Hill's equation. It was discovered by an American astronomer and mathematician from Rutgers University George William Hill (1838--1914) in 1886. Hill's differential equation is widely applicable to many problems involving sawtooth variations, rectangular ripple analysis, exponential variations, etc.

A particular case of Hill's differential equation when the function q(t) has the form

\[ q(t) = a_0 - a_1 \cos \omega t \tag{2.2} \]

is known as Mathieu's equation. The problem of vibrations of an elliptical membrane was first seriously studied by the French mathematician Émile Léonard Mathieu (1835--1890) in 1868, in his celebrated paper, “Mémoire sur le mouvement vibratoire d’une membrane de forme elliptique”, J. Math. Pures Appl. 13, 137 (1868). (http://eudml.org/doc/234720). Mathieu proved that the vertical displacement produced from oscillations on an elliptical drumhead is governed, in general, by the 2D Helmholtz equation. Furthermore, Mathieu showed that when the 2D Helmholtz equation is transformed into elliptical coordinates, the differential equations, respectively, for the radial and angular coordinates of the elliptical system become separable. The result is two almost identical linear second-order differential equations, each of which is aptly named the “Mathieu equation”. So, to this end, Mathieu’s equation has the canonical form

\[ \texttt{D}^2 x(t) + \left( \lambda - 2q\,\cos 2t \right) x(t) = 0 , \tag{2.3} \]

where the separation constant λ and the physical parameter q are both real. Since a Wronskian of the Mathieu equation is a constant, if y₁(t) is a solution, then another linearly independent solution will be

\[ y_2 (t) = y_1 (t)\,\int \frac{{\text d}t}{y_1^2 (t)} . \]

Consider the initial value problem for the linear homogeneous equation

\[ \ddot{x} + \left( A - \sin 2t \right) x =0 , \qquad x(0) =1, \quad \dot{x} (0) =0 , \]

where A is a parameter. Here double superdot \( \ddot{x} = \frac{{\text d}^2 x}{{\text d} t^2} \) stands for the second derivative with respect to time variable t. Upon converting the given differential equation to the system of first order equations

\[ \begin{cases} \dot{x} =y , \\ \dot{y} = - \left( A - \sin 2t \right) x , \end{cases} \qquad x(0) =1, \quad y(0) =0 , \]

 

So we see the behavior of solutions depends on the value of parameter A: when it is greater than 1, solutions have different properties. Therefore, such value of A = 1 is called the bifurcation point (or node). Mathematica has a dedicated command NDSolve to find an a approximate solution to the initial value problem. Correspondingly, we use it to plot solutions of the given initial value problem for different values of parameter A.

sol1 = NDSolve[{x''[t] + (1/2 - Sin[2*t])*x[t] == 0, x[0] == 1, x'[0] == 0}, x, {t, 0, 30}, WorkingPrecision -> 24, InterpolationOrder -> All]
a1 = Plot[x[t] /. sol1, {t, 0, 30}, PlotStyle -> {Thick, Orange}]
sol = NDSolve[{x''[t] + (1 - Sin[2*t])*x[t] == 0, x[0] == 1, x'[0] == 0}, x, {t, 0, 30}, WorkingPrecision -> 24, InterpolationOrder -> All]
a2 = Plot[x[t] /. sol, {t, 0, 30}, PlotStyle -> {Thick, Blue}]
sol3 = NDSolve[{x''[t] + (2 - Sin[2*t])*x[t] == 0, x[0] == 1, x'[0] == 0}, x, {t, 0, 30}, WorkingPrecision -> 24, InterpolationOrder -> All]
a3 = Plot[x[t] /. sol3, {t, 0, 30}, PlotStyle -> {Thick, Black}]
Show[a1, a2, a3]

Solution for A = 1/2

Solution for A = 1

Solution for A = 2

All three solutions plotted together

When you give a setting for WorkingPrecision (which is an option for various numerical operations that specifies how many digits of precision should be maintained in internal computations), this typically defines an upper limit on the precision of the results from a computation. But within this constraint you can tell the Wolfram Language how much precision and accuracy you want to get. You should realize that for many kinds of numerical operations, increasing precision and accuracy goals by only a few digits can greatly increase the computation time required. Given a particular setting for WorkingPrecision, each of the functions for numerical operations in the Wolfram Language uses certain default settings for PrecisionGoal (which is an option for various numerical operations that specifies how many effective digits of precision should be sought in the final result) and AccuracyGoal (which is an option for various numerical operations which specifies how many effective digits of accuracy should be sought in the final result). Typical is the case of NDSolve, in which these default settings are equal to half the settings given for WorkingPrecision.

One can try to control the accuracy of numerical calculations by introducing the residual variable: \( z (t) = \ddot{x} = - \left( A - \sin 2t \right) x = \left( \sin 2t - A \right) x . \) Upon its differentiation, we get another system of differential equations with three unknowns:

\[ \begin{cases} \dot{x} =y , \\ \dot{y} = z , \\ \dot{z} = 2x\,\cos 2t + \left( \sin 2t - A \right) y , \end{cases} \qquad x(0) =1, \quad y(0) (0) =0 , \quad z(0) = -A. \]

This system of 3 differential equations is not equivalent to the original Mathieu equation because its general solution depends on 3 arbitrary constants and its existence requires a solution to be from the space C³. To eliminate extra constant, we need to impose the third initial condition, so the initial value problem will have the same solution as the original IVP for the Mathieu equation.

⁎ ✱ ✲ ✳ ✺ ✻ ✼ ✽ ❋
◾

We can assign to equation \eqref{EqConv.5} the differential operator \( L\left[ t, \texttt{D} \right] = \texttt{D}^n + a_{n-1} (t)\, \texttt{D}^{n-1} + \cdots + a_1 (t)\,\texttt{D} + a_0 (x)\,\texttt{I} , \) where \( \texttt{D} = {\text d}/{\text d}t \) and \( \texttt{D}^0 = \texttt{I} \) are the derivative operator and the identical operator, respectively. Then we rewrite the differential equation in a more concise form:

The same approach can be used to convert a nonlinear differential equation to an equivalent system of first order equations. It should be noted that nonlinear systems of differential equations can be solved only numerically (except for some rare cases). When a numerical solver is applied, it is common to introduce an auxiliary variable---the residual---in order to control accuracy. This approach is demonstrated in the following series of examples.

Example 4: Consider the initial value problem for the ideal pendulum equation

\[ \ddot{\theta} + \omega^2 \sin \theta =0 , \qquad \theta (0) = a = \frac{\pi}{3} \approx 1.0472 , \quad \dot{\theta} (0) = v = -\frac{1}{10}, \tag{4.1} \]

where the initial displacement 𝑎 and velocity v are specified, and ω² = (g/ℓ), with g being local acceleration due to gravity and ℓ being length of the rod containing a point mass. This equation is called the ideal pendulum equation because friction is not taken into account, the weight of the wire is neglected, and oscillation is assumed to be two-dimensional. When resistance is presented and assumed to be proportional to the velocity of the bob, the equation becomes

\[ \ddot{\theta} + k\,\dot{\theta} + \omega^2 \sin \theta =0 , \qquad \theta (0) = a , \quad \dot{\theta} (0) = v , \tag{4.2} \]

where k is a positive coefficient.

These problems (for ideal pendulum and pendulum with resistance) can be converted to equivalent systems of differential equations of the first order:

\[ \begin{cases} \dot{x} = y , \\ \dot{y} = - \omega^2 \sin x , \end{cases} \qquad x(0) =a, \quad y(0) = v ; \tag{4.3} \]

and

\[ \begin{cases} \dot{x} = y , \\ \dot{y} = - \omega^2 \sin x - k\,y , \end{cases} \qquad x(0) =a, \quad y(0) = v , \tag{4.4} \]

respectively. We plot its phase portraits using Mathematica:

LineIntegralConvolutionPlot[{{y, -Sin[x]}, {"noise", 1000, 1000}}, {x, -6, 6}, {y, -2, 2}, ColorFunction -> Hue, LightingAngle -> 0, LineIntegralConvolutionScale -> 3, Frame -> False]

Phase portrait for undriven ideal pendulum

Phase portrait for undriven pendulum with resistance

 

The ideal pendulum equation can be converted to a system of three equations upon introducing the residual variable:

\[ \begin{cases} \dot{x} = y , \\ \dot{y} = z , \\ \dot{z} = - \omega^2 y\,\cos x , \end{cases} \qquad x(0) =a, \quad y(0) = v , \quad z(0) = - \omega^2 \sin a . \]

So z is actually the second derivative of the original unknown function x(t). Therefore, the last equation is obtained by differentiation of the original pendulum equation with respect to time variable t. WE control approximation of the pendulum differential equation with the following Mathematica code:

		Since the function z(t) is supposed to be the second derivative, we can control the accuracy of Mathematica approximation of the pendulum equation by evaluating \( \ddot{\theta} + \sin \theta . \) NDSolve[{x'[t] == y[t], y'[t] == z[t], z'[t] == -y[t]*Cos[x[t]], x[0] == Pi/3, y[0] == -1/10, z[0] == -Sin[Pi/3]}, {x, y, z}, {t, 0, 30}]; Plot[Evaluate[(z[t] + Sin[x[t]]) /. s], {t, 0, 10}, PlotRange -> All]
Plot of \( \ddot{\theta} + \sin \theta . \)		Mathematica code

   ◾

Example 5: In the previous example, we converted the ideal pendulum equation \( \ddot{\theta} + \sin \theta =0 \) into the system of first order equations upon substituting y₁ = θ and let y₂ to be its derivative. Now we may define a new variable y₃ = sinθ. According to the chain rule, its derivative is

\[ \frac{{\text d}y_3}{{\text d}t} = \left(\frac{{\text d}y_1}{{\text d}t} \right) \cos y_1 = y_2 \cos \theta . \]

Then we introduce a fourth dependent variable y₄ = cosθ and obtain for it

\[ \frac{{\text d}y_4}{{\text d}t} = - \left(\frac{{\text d}y_1}{{\text d}t} \right) \sin y_1 = - y_2 \,y_3 . \]

Finally, we combine all equations together:

\begin{align*} \dot{y}_1 &= y_2 , \\ \dot{y}_2 &= -\omega^2 y_3 , \\ \dot{y}_3 &= y_2\, y_4 , \tag{5.1} \\ \dot{y}_4 &= -y_2\, y_3 . \end{align*}

The initial conditions for y₃ and y₄ will be related to those imposed on y₁ and y₂:

\[ y_1 (0) = a = \frac{\pi}{3}, \quad y_2 (0) = v = -\frac{1}{10}, \quad y_3 (0) = \sin a = \frac{\sqrt{3}}{2} \approx 0.866025, \quad y_4 (0) = \cos a = \frac{1}{2} . \tag{5.2} \]

Picard's iteration (see section in this tutorial) for this system of equations is

\[ {\bf \phi}_{m+1} (t) = \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \\ y_4 (t) \end{bmatrix}_{(m+1)} = \begin{bmatrix} a \\ v \\ \sin a \\ \cos a \end{bmatrix} + \int_0^t {\text d}s \begin{bmatrix} y_2 (s) \\ - \omega^2 y_3 (s) \\ y_2 (s)\,y_4 (s) \\ -y_2 (s)\,y_3 (s) \end{bmatrix}_{(m)} , \qquad m= 0,1,2,\ldots . \tag{5.3} \]

Since the right-hand side in Eq.(5.3) is a polynomial, we can apply Picard's iteration. Performing a couple of iterations, we obtain

\[ \theta_{(3)} (t) = a + tv - \frac{t^2}{2}\, \sin a - \frac{t^3}{6}\,v\,\cos a . \]

a = Pi/3; v = -1/10;
s = NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == Pi/3, x'[0] == -1/10}, x, {t, 0, 20}];
phi0 = {a, v, Sqrt[3]/2, 1/2};
phi1[t_] = phi0 + t*{v, -w^2 *Sqrt[3]/2, v/2, -v*Sqrt[3]/2} phi2[t_] = phi0 + Integrate[{phi1[t][[2]], -w^2*phi1[t][[3]], phi1[t][[2]]*phi1[t][[4]], -phi1[t][[2]]*phi1[t][[3]]}, {t, 0, t}];
phi3[t_] = phi0 + Integrate[{phi2[t][[2]], -w^2*phi2[t][[3]], phi2[t][[2]]*phi2[t][[4]], -phi2[t][[2]]*phi2[t][[3]]}, {t, 0, t}];
phi4[t_] = phi0 + Integrate[{phi3[t][[2]], -w^2*phi3[t][[3]], phi3[t][[2]]*phi3[t][[4]], -phi3[t][[2]]*phi3[t][[3]]}, {t, 0, t}];
phi[t_] = Pi/3 - t/10 - (Sqrt[3] t^2)/4 + t^3/120 + ( 17 Sqrt[3] t^4)/1600 + t^5/800;

a + t v - 1/6 t^3 v Cos[a] - 1/2 t^2 Sin[a]

Now we plot the fourth approximation along with the true solution:

Plot[{Evaluate[x[t] /. s], phi[t]}, {t, 0, 5}, PlotTheme -> "Web", Epilog -> {Text[Style["approximation", 25, Blue], {3.1, 2}]}]

   ◾

Example 6: Consider the initial value problem for Duffing's equation

\[ \ddot{x} + k\,\dot{x} + x - \varepsilon\,x^3 = f(t) , \qquad x(0) =a, \quad \dot{x} = v , \tag{6.1} \]

where f(t) is a periodic forcing function, which can be chosen as \( f(t) = A\,\cos \omega t \) or \( f(t) = A\,\sin \omega t \) or their combination \( f(t) = A\,\sin \left( \omega t + \phi \right) . \) Here A is the amplitude, k is the dumping coefficient, and ε is a real number.

We convert it to the system of two differential equations

\[ \begin{cases} \dot{x} = y , \\ \dot{y} = - k\,y -x + \varepsilon\,x^3 + f(t) , \end{cases} \qquad x(0) =a, \quad y(0) = v . \tag{6.2} \]

Its direction field for undriven case is plotted with the following commands:

StreamPlot[{y, -y - x - x^3}, {x, -6, 6}, {y, -2, 2}, StreamScale -> Medium]

The residual system of equations becomes

\[ \begin{cases} \dot{x} = y , \\ \dot{y} = z , \\ \dot{z} = - k\,z - y - 3\varepsilon \, x^2 y - \dot{f} , \end{cases} \qquad x(0) =a, \quad y(0) = v , \quad z(0) = - k\,v - a - \varepsilon \, a^3 + f(0). \]

   ◾