Preface

This section presents some applications of power series method in numerical approximation of solutions to systems of ordinary differential equations.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part IV of the course APMA0340
Introduction to Linear Algebra

Glossary

Power Series Method

The most differential equations can’t be solved explicitly in terms of finite combinations of simple familiar functions. In this section, we develop an algorithm for solving a certan class of differential equations based on the power series method. In general, such a solution assumes a power series with unknown coefficients, then substitutes that solution into the differential equation to find a recurrence relation for the coefficients.

Power series solution method has been traditionally used to solve linear differential equations. Although the power series method is not, generally speaking, suitable for direct integration of many differential equations, it has many modifications and it is a part of theoretical analysis and practical numerical procedures. The main obstacle in its direct application lies in difficulty of determination of the radius of convergence. It is well known that a power series \( S(x) = \sum_{n\ge 0} a_n \left( x- x_0 \right)^n \) converges within a circle \( \left\vert x - x_0 \right\vert < R \) of radius R. However, its value is determined by behavior of coefficients \( a_n \) at infinity. Since \( R^{-1} = \lim_{n\to \infty} \sqrt[n]{a_n} , \) the radius of convergence depends on how fast the coefficients grow. However, if sum-function S(x) is unknown because in our applications it is usually the unknown solution of the differential equation, it is impossible to find the radius of convergence.

Nevertheless, the power series method can be used practically at least is some small neighborhood of the center. We are going to demonstrate this approach in a sequence of examples.

 

Example: We analyze the rotation of a free rigid body that is suspended in space, for example, a satellite or a planet. Recall that a rigid body is a collection of finite points that the distance between the points is fixed. So we consider a body that is rotated about some point, which could take as the center of mass. Then the position of a particle in the body can be written as

\[ {\bf r}_i (t) = {\bf R}(t) + \Delta{\bf r}_i (t) , \]
where Δri is the position measured from the center of mass. For any point r in the body
\begin{align*} \frac{{\text d} {\bf r}}{{\text d}t} &=r_a \,\frac{{\text d}R_{ab}}{{\text d}t}\, {\bf e}_b \qquad \mbox{in the space frame} \\ &= r_a \,\frac{{\text d}{\bf e}_1 (t)}{{\text d}t} \qquad \mbox{in the body frame} . \end{align*}
The components of the orthogonal matrix \( \displaystyle \left[ R_{ab} \right] \) satisfy \( \displaystyle {\bf e}_a (t) = R_{ab} (t)\,{\bf e}_b . \) The angular velocity vector ω(t) has components ωaea. Then examining the kinetic energy, we get
\begin{align*} {\mbox K} &= \frac{1}{2}\,\sum_i m_i \dot{r}_i^2 \\ &= \sum_i m_i \left[ \frac{1}{2} \,\dot{\bf R}^2 + \dot{\bf R} \cdot \left( \omega \times \Delta {\bf r}_i \right) + \frac{1}{2} \left( \omega \times \Delta {\bf r}_i \right)^2 \right] \\ &= \frac{1}{2}\, M\,\dot{\bf R}^2 + \frac{1}{2} \,\omega_a I_{ab} \omega_b , \end{align*}
because \( \displaystyle \sum_i m_i \Delta {\bf r}_i = 0. \) So we find that the dynamics separates into the motion of the center of mass R, together with rotation about the center of mass. For rotational process, the angular momentum must be conserved:
\[ \frac{{\text d} {\bf L}}{{\text d}t} = 0 , \]
where \( \displaystyle {\bf L} = L_1 {\bf e}_1 + L_2 {\bf e}_2 + L_3 {\bf e}_3 = I_1 \omega_1 {\bf e}_1 + I_2 \omega_2 {\bf e}_2 + I_3 \omega_3 {\bf e}_3 , \) Ik are the principal moments of inertia and ωk are the components of the angular velocity about the principal axes. Using \( \displaystyle L_a = I_{ab} \omega_b , \) we apply to every point and obtain the equations of motion (Euler's equations)
\begin{align*} I_1 \dot{\omega}_1 &= \omega_3 L_2 - \omega_2 L_3 = \omega_2 \omega_3 \left( I_2 - I_3 \right) , \\ I_2 \dot{\omega}_2 &= \omega_1 L_3 - \omega_3 L_1 = \omega_1 \omega_3 \left( I_3 - I_1 \right) , \\ I_3 \dot{\omega}_3 &= \omega_2 L_1 - \omega_1 L_2 = \omega_1 \omega_2 \left( I_1 - I_2 \right) . \end{align*}
Dividing by moments of inertia, we obtain (Euler, 1749)
\begin{align*} \dot{\omega}_1 &= \frac{I_2 - I_3}{I_1} \, \omega_2 \omega_3 , \\ \dot{\omega}_2 &= \frac{I_3 - I_1}{I_2} \, \omega_1 \omega_3 , \\ \dot{\omega}_3 &= \frac{I_1 - I_2}{I_3} \, \omega_1 \omega_2 . \end{align*}
The spin of the Earth causes it to bulge at the equator so it is no longer a sphere but can be treated as a symmetric top. It is an oblate ellipsoid, with I3 > I1 and is spherical to roughly 1 part in 300, meaning
\[ \frac{I_1 - I_3}{I_1} \approx \frac{1}{300} . \]
Of course, we know the magnitude of the spin ω3 is just 1 (day)-1. This information is enough to calculate the frequency of the Earth’s wobble, which takes 433 days and is known as the Chandler wobble.

Suppose that the body is freely rotating symmetrically about the z-axis. Then Ixx = Iyy = I⊥. Likewise, we can write Izz = I∥ ≠ I⊥, in general. In this case, Euler's equations become

\begin{align*} \dot{\omega}_1 &= \frac{I_{\perp} - I_{\parallel}}{I_{\perp}} \, \omega_2 \omega_3 , \\ \dot{\omega}_2 &= \frac{I_{\parallel} - I_{\perp}}{I_{\perp}} \, \omega_1 \omega_3 , \\ \dot{\omega}_3 &= 0 . \end{align*}
Clearly, ω3 = ωz is a constant, so the above equations yield
\begin{align*} \dot{\omega}_x &= - \Omega \, \omega_y , \\ \dot{\omega}_y &= \Omega \, \omega_x , \end{align*}
where \( \displaystyle \omega_z \left( I_{\parallel}/I_{\perp} -1 \right) . \) Since this system of equations is linear, its solution is not hard to determine:
\begin{align*} \omega_x &= - \omega \, \cos \left( \Omega t \right) , \\ \omega_y &= \omega \, \sin \left(\Omega t \right) , \end{align*}
where ω is a constant. Thus, the projection of the angular velocity vector onto the x-y plane has the fixed length ω, and rotates steadily about the z-axis with angular velocity Ω.    ■

 

Example: The famous Lorenz system of equations is given as

\begin{align*} \dot{x} &= \sigma \left( y - x \right) , \\ \dot{y} &= r\,x -y -x\,z , \\ \dot{z} &= x\,y -b\,z , \end{align*}
where x, y, and z are respectively proportional to the convective velocity, the temperature difference between descending and ascending flows, and the mean convective heat flow; r is the bifurcation parameter, and the real coefficients σ and b are taken the numerical values r = 28, σ = 10, and b = 8/3 where the system exhibits chaotic behavior.

The explicit solution to the Lorenz system we assume to be represented by power series:

\begin{align*} x(t) &= \sum_{k\ge 0} \frac{a_k}{k!} \left( t - t^{\ast} \right)^k , \\ y(t) &= \sum_{k\ge 0} \frac{b_k}{k!} \left( t - t^{\ast} \right)^k , \\ z(t) &= \sum_{k\ge 0} \frac{c_k}{k!} \left( t - t^{\ast} \right)^k , \end{align*}
where the coefficients are given by the recurrence relations
\begin{align*} a_0 &= x(t^{\ast} ) , \qquad b_0 = y(t^{\ast} ) , \qquad c_0 = z(t^{\ast} ) , \\ a_k &= - \sigma \,a_{k-1} + \sigma\,b_{k-1} , \qquad k \ge 1, \\ b_k &= r\,a_{k-1} + b_{k-1} - (k-1)! \sum_{i=0}^{k-1} \frac{a_i c_{k-i-1}}{i!\,(k-i-1)!} , \qquad k \ge 1, \\ c_k &= -b\,c_{k-1} + (k-1)! \sum_{i=0}^{k-1} \frac{a_i b_{k-i-1}}{i! \,(k-i-1)!} , \qquad k\ge 1 . \end{align*}
These recurrence relations are used for algorithmic computations over time mesh.    ■

 

  1. Abbasbandy, S. and ervillier, C., Analytic continuations of Taylor series and the two-point boundary value problem of some nonlinear ordinary differential equations, [arXiv:1104.5073v1] (2011).
  2. Asaithambi, A., A series solution of the Falkner–Skan equation using the crocco–wang transformation, International Journal of Modern Physics, 2017, Vol. 28, No. 11, 1750139; https://doi.org/10.1142/S012918311750139X
  4. Noor, N.F.M., Hashim, I., Noorani, M.S.M., A Note On The Accuracy Of The Adomian Decomposition Method Applied To The Chaotic Lorenz System, Proceedings of the 2nd IMT-GT Regional Conference on mathematics and Applications, University Saint Malaysia, Penang, June 13--15, 2006.

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions