Historically, Picard's iteration scheme was the first method to solve analytically nonlinear differential equations, and it was discussed in the first part of the course (see introductory section xv Picard). In this section, we widen this procedure for systems of first order differential equations written in normal form \( \dot{\bf x} = {\bf f}({\bf x}, t) . \) Although this method is rarely used for actual evaluations due to slow convergence and obstacles with performing explicit integration, it is very important for understanding material. Working with Picard's iterations and its refinements helps everyone to develop computational skills. This section also gives some motivated examples about transferring differential equations with complicated slope functions to systems of ODEs with polynomial driven terms that are suitable for Picard's iteration.

This section expands Picard's iteration process to systems of ordinary differential equations in normal form (when the derivative is isolated). Picard's iterations for a single differential equation \( {\text d}x/{\text d}t = f(t,x) \) was considered in detail in the first tutorial (see section for reference). Therefore, our main interest would be applying Picard's iteration to systems of first order ordinary differential equations in normal form (which means that the derivative is isolated)

\[ \begin{cases} \dot{x}_1 &= f_1 (x_1 , x_2 , \ldots , x_n , t) , \\ \dot{x}_2 &= f_2 (x_1 , x_2 , \ldots , x_n , t) , \\ \quad \vdots & \quad \vdots \\ \dot{x}_n &= f_n (x_1 , x_2 , \ldots , x_n , t) . \end{cases} \]

This explicit notation cannot be called compact and informative. Hence, we introduce n-dimensional column-vectors of unknown variables and the slope function in order to rewrite the system above in compact form. So the system of differential equations can be written in vector form as

\[ \frac{{\text d} {\bf x}}{{\text d}t} = {\bf f}({\bf x} , t ) , \]

where

\[ {\bf x} (t) = \begin{bmatrix} x_1 (t) \\ x_2 (t) \\ \vdots \\ x_n (t) \end{bmatrix} , \qquad {\bf f} (x_1 , x_2 , \ldots , x_n , t) = \begin{bmatrix} f_1 (x_1 , x_2 , \ldots , x_n ,t) \\ f_2 (x_1 , x_2 , \ldots , x_2 , t) \\ \vdots \\ f_n (x_1 , x_2 , \ldots , x_n , t) \end{bmatrix} \]

are n-column vectors. Note that the preceding system of equations contains n dependent variables \( x_1 (t), x_2 (t) , \ldots , x_n (t) \) while the independent variable is denoted by t, which may be associated with time. In engineering and physics, it is a custom to follow Isaac Newton and denote a derivative with respect to time variable t by dot: \( \dot{\bf x} = {\text d}{\bf x} / {\text d} t. \) When input function f satisfies some general conditions (usually required to be Lipschitz continuity in the dependent variable \( {\bf x} \) ), then Picard's iteration converges to a solution of the Volterra integral equation in some small neighborhood of the initial point. Thus, Picard's iteration is an essential part in proving existence of solutions for the initial value problems.

Note that Picard's iteration is not suitable for numerical calculations. The reason is not only slow convergence, but mostly it is impossible, in general, to perform explicit integration to obtain the next iteration. Although in case of a polynomial input function, integration can be performed explicitly (especially, with the aid of a computer algebra system), and the number of terms quickly grows as a snow ball. While we know that the resulting series converges eventually to the true solution, its range of convergence is too small to keep many iterations. In another section, we show how to bypass this obstacle.

If an initial position of the vector x(t) is known, we get an initial value problem:

\begin{equation} \label{EqPicard.1} \frac{{\text d} {\bf x}}{{\text d}t} = {\bf f}({\bf x},t ) , \qquad {\bf x} (t_0 ) = {\bf x}_0 , \end{equation}

where x₀ is an initial column vector. Many brilliant mathematicians participated in proving the existence of a solution to the given initial value problem more than 100 years ago. Their proof was based on what is now called the Picard's iteration, named after the French mathematician Charles Émile Picard (1856--1941) whose theories did much to advance research in analysis, algebraic geometry, and mechanics.

The Picard procedure is actually a practical extension of the Banach fixed point theorem, which is applicable to continuous contractible functions. Since any differential equation involves an unbounded derivative operator, the fixed point theorem is not suitable for it. To bypass this obstacle, Picard suggested to apply the (bounded) inverse operator L^-1 to the first derivative \( \texttt{D} . \) Recall that the inverse \( \texttt{D}^{-1} , \) called in mathematical literature as «antiderivative», is not an operator because it assigns to every input-function infinite many outputs. To restrict its output to a single one, we consider the differential operator on the set of functions (which becomes a vector space only when the differential equation and the initial condition are all homogeneous) with a specified initial condition f(x₀) = y₀. So the derivative operator on this set of functions we denote by L, and its inverse is a bounded integral operator. The first step in deriving Picard's iterations is to rewrite the given initial value problem in equivalent (this is true when the slope function f is continuous in Lipschitz sense) form as Volterra integral equation of second kind:

\begin{equation} \label{EqPicard.2} {\bf x} ( t ) ={\bf x}_0 + \int_{t_0}^t {\bf f} ( {\bf x} (s), s ) \,{\text d}s . \end{equation}

This integral equation is obtained upon integration of both sides of the differential equation \( {\text d}{\bf x} / {\text d} t = {\bf f}(t, {\bf x}), \) subject to the initial condition.

The equivalence follows from the Fundamental Theorem of Calculus. It suffices to find a continuous function \( {\bf x}(t) \) that satisfies the integral equation within the interval \( t_0 -h < t < t_0 +h , \) for some small value h since the right hand-side (integral) will be continuously differentiable in t.

Assuming that the vector-function \( {\bf f} ({\bf x}, t) \) satisfies the Lipschitz condition in x:

\[ \| {\bf f}({\bf x} , t ) - {\bf f}({\bf y} ,t ) \| \le L\,\| {\bf x} - {\bf y} \| , \]

where constant L is called the Lipschitz constant, then the initial value problem \eqref{EqPicard.1} has a unique solution in some neighborhood of the initial position.

Now we apply a technique known as Picard iteration to construct the required solution:

\begin{equation} \label{EqPicard.3} {\bf x}_{m+1} ( t ) ={\bf x}_0 + \int_{t_0}^t {\bf f} ({\bf x}_m (s) , s) \,{\text d}s , \qquad m=0,1,2,\ldots . \end{equation}

The initial approximation is chosen to be the initial value (constant): \( {\bf x}_0 (t) \equiv {\bf x}_0 .\) (The sign ≡ indicates the values are identically equivalent, so this function is the constant). When input function f(t, x) satisfies some (sufficient) conditions, it can be shown that this iterative sequence converges to the true solution uniformly on the interval \( [t_0 -h, t_0 +h ] \) when h is small enough. Actually, the necessary and sufficient conditions for existence of a solution to a vector initial value problem are still waiting for their discovery.

Of course, we can differentiate the recursive integral relation to obtain a sequence of the initial value problems

\begin{equation} \label{EqPicard.4} \frac{{\text d}{\bf x}_{m+1} ( t )}{{\text d}t} = {\bf f} (t, {\bf x}_m (t) , \qquad {\bf x}_m (t_0 ) = {\bf x}_0 , \qquad m=0,1,2,\ldots . \end{equation}

that is equivalent to the original recurrence \eqref{EqPicard.3}. However,this recurrence relation \eqref{EqPicard.4} involves the derivative operator \( \texttt{D} = {\text d}/{\text d}t . \) Since \( \texttt{D} \) is an unbounded operator, it is not possible to prove convergence of the iteration procedure \eqref{EqPicard.4}. On the other hand, the convergence of \eqref{EqPicard.3} involving a bounded integral opeprator can be accomplished directly (see Tutorial I).

Example 1: . Using Newton's laws of motion and gravitation you can in theory work out how a number of objects (say stars, planets and moons) should move due to the gravitational pull they exert on each other. If there are only two bodies involved the answer is straight-forward: they move along curves that are easy to describe (ellipses in the case of a planet orbiting a star). The problem of three or more bodies motivated the science community and Henri Poincaré solved this three body problem on November 30th 1889, unfortunately, withe a mistake. He revised his competition entry, still got the prize and was paid, and later wrote more extensively about what the mistake had revealed to him---the existence of chaos.

Let x_ij and v_ij represent the i-th position coordinate and the i-th velocity component, respectively, of particles j = 1, 2, …, N of respective masses m_j. Furthermore, let \( r_{jk} = \left\vert \sum_{i=1}^3 \left( x_{ij} - x_{ik} \right)^2 \right\vert^{1/2} \) denote the separation distances between the centers of particles j and k, where k = 1, 2, …, N. In classical form, the N-body problem consists of 6N coupled autonomous first-order ODEs, namely \begin{equation*} \label{EqNbody.1} \begin{split} \frac{{\text d} x_{ij}}{{\text d}t} = v_{ij} , \\ \frac{{\text d} v_{ij}}{{\text d}t} = \sum_{k\ne j} m_k \left( x_{ik} - x_{ij} \right) / r_{jk}^3 . \end{split} \tag{1.1} \end{equation*} By defining the inverse separation distances to be \( s_{jk} = 1 / r_{jk} \ (j\ne k) , \) we derive the following polynomial formulation of the N-body problem: \begin{equation*} \label{EqNbody.2} \begin{cases} \frac{{\text d} x_{ij}}{{\text d}t} = v_{ij} , \\ \frac{{\text d} v_{ij}}{{\text d}t} = \sum_{k\ne j} m_k \left( x_{ik} - x_{ij} \right) * s_{jk}^3 , \\ \frac{{\text d} s_{ij}}{{\text d}t} = - s_{jk}^3 \sum_{i=1}^3 \left( x_{ik} - x_{ij} \right) \left( v_{ik} - v_{ij} \right) . \end{cases} \tag{1.2} \end{equation*} Even when symmetry is exploited ( \( s_{jk} = s_{kj} \) ), we have paid a price for polynomial form; the number of equations n has been increased to n = 6N +N(N - 1)/2. The gain is that the system now readily admits a Maclaurin-series solution of arbitrary order.

Let x_ijl denote the l-th order coefficient of the Maclaurin series for x_ij. That is, \( x_{ij}(t) = \sum_{l\ge 0} x_{ijl} t^l . \) Let similar notation hold for v_ij and s_ij. In light of the theoretical considerations mentioned previously, the coefficients for the m-th term are obtained simply by integrating (1.2) with respect to t, to yield \begin{equation*} \label{EqNbody.3} x_{ijm} = v_{i,j,m-1}/m . \tag{1.3} \end{equation*} Similarly, we have \begin{equation*} \label{EqNbody.4} v_{ijm} = \sum_{k=1}^n m_k \sum_{l=0}^{m-1} \left( x_{pkl} - x_{ijl} \right) \left( s^3 \right)_{j,k,m-l-1} /m , \tag{1.4} \end{equation*} where \begin{equation*} \label{EqNbody.5} \left( s^2 \right)_{jkm} = \sum_{l=0}^m s_{jkl} s_{j,k,l-m} , \qquad \left( s^3 \right)_{jkm} = \sum_{l=0}^m \left( s^2 \right)_{ijk} s_{j,k,l-m} \tag{1.5} \end{equation*} are convolutions of sequences s with itself. Finally, by defining \begin{equation*} \label{EqNbody.6} a_{jkm} = \sum_{l=1}^m \sum_{i=1}^3 \left( x_{ijl} - x_{ikl} \right) \left( v_{i,j,m-l} - v_{i,k,m-l} \right) , \tag{1.6} \end{equation*} we obtain from \begin{equation*} \label{EqNbody.7} s_{jkm} = - \sum_{l=1}^{m-1} \left( s^3 \right)_{jkl} a_{j,k,m-l} /m . \tag{1.7} \end{equation*}

The three-body problem with close encounters is notoriously ill-conditioned because it admits chaotic solutions that manifest extreme sensitivity to initial conditions. The following example demonstrates the Picard's method for an ill-conditioned three-body problem comprised of the earth, its moon,and a ‘‘spacecraft’’.

The numerical values of three body problem are presented in the following table. Data are only approximately related to realistic values. For all bodies, the initial coordinate x₃ and velocity v₃ are assumed zero so that the orbits remain co-planar. Masses are scaled relative to the mass of the earth; distances are in earth radii. Velocity is normalized by the orbital velocity at the earth's surface, for which the orbital period is 2π. The integration time (T = 3200) is in excess of 500 orbital periods.

Body	Mass, M_j	x₁	x₂	v₁	v₂
Earth	1.0	0.0	0.0	0.0	0.0
Moon	0.0125	42.4	−42.4	0.095	0.095
Craft	0.00001	−1.0	0.0	0.0	1.4034

■

Picard's iteration for higher order differential equations

Now we turn our attention to application of Picard's iteration procedure for solving initial value problems for higher order (linear or nonlinear) single differential equations:

\begin{equation} \label{EqPicard.5} x^{(n)} (t) = f\left( t, x, x' , \ldots , x^{(n-1)} (t) \right) , \qquad x(0) = x_0 , \ x' (0) = x_1 , \ \ldots , \ x^{(n-1)} (0) = x_{n-1} , \end{equation}

where \( x^{(n)} (t) = \frac{{\text d}^n x}{{\text d}t^n} = \texttt{D}^n x(t) \) stands for the n-th derivative. There are usually two options to achieve this. First, we can transfer this single differential equation to an equivalent system of first order ODEs upon introducing new dependent variables:

\begin{equation} \label{EqPicard.6} y_1 (t) = x(t) , \quad y_2 (t) = x' (t) = y'_1 (t) , \quad \ldots , \quad y_{n} (t) = y'_{n-1} (t) = x^{(n-1)} (t) . \end{equation}

Rewriting this system of equations in a standard vector form \( \dot{\bf y} = {\bf f}(t, {\bf y}) , \) we can apply Picard's iterations for column-vector y.

Another approach is based on extension of Picard's idea for a higher order equation. It is based on the observation that the initial value problem for n-th order derivative

\begin{equation} \label{EqPicard.7} \texttt{D}^n x(t) = f(t, x, x' , \ldots , x^{(n-1)}) , \quad x(0) = x_0 , \ x' (0) = x_1 , \ \ldots , \ x^{(n-1)} (0) = x_{n-1} , \end{equation}

where \( \texttt{D} = \frac{{\text d}}{{\text d}t} \) is the derivative operator, has an explicit solution (which is unique)

\begin{equation} \label{EqPicard.8} x(t) = \frac{1}{(n-1)!} \int_0^t \left( t-s \right)^{n-1} f(s, x(s) , \ldots )\,{\text d}s + \sum_{k=0}^{n-1} \frac{x_k}{k!}\, t^k . \end{equation}

Upon denoting by ϕ_m(t) the m-th Picard iterate, we write the Picard iteration scheme:

\begin{equation} \label{EqPicard.9} \phi_{m+1} (t) = \frac{1}{(n-1)!} \int_0^t \left( t-s \right)^{n-1} f\left( s, \phi_m (s) , \ldots \right)\,{\text d}s + \sum_{k=0}^{n-1} \frac{a_{k}}{k!}\, t^k , \qquad m=n, n+1, \ldots , \end{equation}

where coefficients \( a_k \) and first n initial terms ϕ₀ , … , ϕ_n-1 are determined from the initial conditions.

Since the convergence of the successive approximations to the true solution is often slow, it could be beneficial to choose the initial approximation wisely. All approaches are demonstrated in a series of examples.

Example 2: Consider the initial value problem for the first order differential equation

\[ \frac{{\text d}y}{{\text d}x} = \sqrt{\frac{1}{4} + 2\,y} , \qquad y(0) = 1. \tag{2.1} \]

We apply Picard's iteration directly by solving the fixed point problem

\[ \phi_{n+1} (x) = 1 + \int_0^x \sqrt{\frac{1}{4} + 2\,\phi_n (t)}\,{\text d}t , \qquad n=0,1,2,\ldots . \tag{2.2} \]

If we start with the initial condition ϕ₀ = 1, we almost immediately come to a problem of explicit integration. Indeed, we calculate first the few terms:

\begin{align*} \phi_0 &= 1 , \\ \phi_1 (x) &= 1 +\int_0^x \sqrt{\frac{1}{4} + 2}\,{\text d}t = 1 + \frac{3}{2}\, x , \\ \phi_2 (x) &= 1 +\int_0^x \sqrt{\frac{1}{4} + 2\,\phi_1 (t)}\,{\text d}t = 1+ \int_0^x \sqrt{\frac{1}{4} + 2 +3\,t}\,{\text d}t = \frac{x}{\sqrt{3}}\, \sqrt{3 + 4\,x} + \frac{1}{4} \left( 1 + \sqrt{9 + 12\,x} \right) , \end{align*}

because SymPy helps:

Although, the next iteration Mathematica can handle by invoking special functions,

Simplify[1 + Assuming[x > 0, Integrate[ Sqrt[t Sqrt[1 + (4 t)/3] + 1/4 (1 + Sqrt[9 + 12 t])], {t, 0, x}]]

it will not overcome the further step.

This example illustrates that Picard's scheme is not practical because the corresponding iteration cannot be performed indefinitely---the functions are not explicitly integrable. In order to apply Picard's method, the slope function must be a polynomial---otherwise the results of integration cannot be expressed through elementary functions. However, in many practical cases, this obstacle can be bypassed with introduction of auxiliary dependent variables and converting the given equation to a system of differential equations.

To apply Picard's method, we convert the given initial value problem to the system of differential equations upon introducing auxiliary dependent variables:

\begin{align*} y_1 (x) &= y(x), \\ y_2 (x) &= \sqrt{\frac{1}{4} + 2\,y(x)} , \\ y'_1 &= y_2 (x) , \\ y'_2 (x) &= \frac{1}{2\,\sqrt{\frac{1}{4} + 2\,y(x)}} \, 2\, y' (x) = 1 . \end{align*}

We get an equivalent to the given autonomous differential equation of the first order the vector equation with polynomial right-hand side input:

\[ \frac{\text d}{{\text d}x} \begin{bmatrix} y_1 (x) \\ y_2 (x) \end{bmatrix} = \begin{bmatrix} y_2 (x) \\ 1 \end{bmatrix} , \qquad \begin{bmatrix} y_1 (0) \\ y_2 (0) \end{bmatrix} = \begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} . \tag{2.3} \]

Integrating the vector equation above allows us to rewrite the given initial value problem in equivalent vector form:

\[ \begin{bmatrix} y_1 (x) \\ y_2 (x) \end{bmatrix} = \begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} + \int_0^x \begin{bmatrix} y_2 (s) \\ 1 \end{bmatrix} {\text d} s . \]

The Picard's iteration scheme will be

\[ \begin{bmatrix} \phi_1 (x) \\ \phi_2 (x) \end{bmatrix}_{(m+1)} = \begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} + \int_0^x \begin{bmatrix} \phi_2 (s) \\ 1 \end{bmatrix}_{(m)} {\text d} s , \qquad m=0,1,2, \ldots . \]

Since Picard's iteration always starts with the initial condition, we get

\[ \begin{bmatrix} \phi_1 (x) \\ \phi_2 (x) \end{bmatrix}_{(1)} = \begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} + \int_0^x \begin{bmatrix} \frac{3}{2} \\ 1 \end{bmatrix} {\text d} s = \begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} + \begin{bmatrix} \frac{3}{2}\, x \\ x \end{bmatrix} = \begin{bmatrix} 1+\frac{3}{2}\, x \\ \frac{3}{2} + x \end{bmatrix} . \]

Next term will be

\[ \begin{bmatrix} \phi_1 (x) \\ \phi_2 (x) \end{bmatrix}_{(2)} = \begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} + \int_0^x \begin{bmatrix} \frac{3}{2} +s \\ 1 \end{bmatrix} {\text d} s = \begin{bmatrix} 1+ \frac{3}{2}\, x + \frac{x^2}{2} \\ \frac{3}{2} +x \end{bmatrix} . \]

The second component, ϕ₂, remains the same with every iteration. This component is used to generate the first component, ϕ₁, as well. As such, every following iteration will be the same, and the true solution becomes

\[ y(x) = 1 +\frac{3}{2}\, x + \frac{x^2}{2} . \]

Its derivative is

\[ y' (x) = \frac{3}{2} + x \qquad \Longrightarrow \qquad \frac{1}{4} + 2\, y(x) = \left( \frac{3}{2} + x \right)^2 . \]

■

Example 3: Consider the initial value problem for the simple pendulum equation

\[ \frac{{\text d}^2 \theta}{{\text d} t^2} + \sin \theta = 0 , \qquad \theta (0) = \frac{\pi}{4} , \quad \dot{\theta} (0) = \frac{1}{10} , \tag{3.1} \]

where dot stands for the derivative with respect to time variable t. We integrate the pendulum equation (3.1) twice, and use initial condition provided. We may try to apply Picard's iteration scheme

\[ \phi_{m+1} (t) = \frac{\pi}{4} + \frac{1}{10}\, t - \int_0^t \left( t-s \right) \sin \phi_m (s) \,{\text d}s , \qquad m= 0,1,2, \ldots . \tag{3.2} \]

The initial function is chosen to match the initial conditions: \( \phi_{0} (t) = \frac{\pi}{4} + \frac{t}{10} . \) Although we can find the next term

\[ \phi_{1} (t) = \frac{\pi}{4} + \frac{1}{10}\, t - \int_0^t \left( t-s \right) \sin \left( \frac{\pi}{4} + \frac{1}{10}\,s \right) {\text d}s = \frac{\pi}{4} + \frac{1}{10}\, t + 5\sqrt{2} \left( 10 +t \right) - 100\,\sin \frac{5\pi + 2t}{20} \]

with the aid of Python

the second term ϕ₂(t) is impossible to determine, even for such a powerful computer algebra system. We can bypass this obstacle by introducing auxiliary variables:

\begin{align*} y_1 (t) &= \theta (t), \\ y_2 (t) &= \dot{y}_1 = \dot{\theta} , \\ y_3 (t) &= \sin \theta (t) , \\ y_4 (t) &= \cos \theta (t) . \end{align*}

Then the pendulum equation and corresponding initial conditions will be equivalent to the following IVP:

\[ \frac{\text d}{{\text d}t} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \\ y_4 (t) \end{bmatrix} = \begin{bmatrix} y_2 (t) \\ -y_3 (t) \\ y_2 (t)\,y_4 (t) \\ -y_2 (t)\,y_3 (t) \end{bmatrix} , \qquad \begin{bmatrix} y_1 (0) \\ y_2 (0) \\ y_3 (0) \\ y_4 (0) \end{bmatrix} = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} . \tag{3.3} \]

Since this system of four differential equations has a polynomial slope part, we expect to perform all steps of Picard's iteration scheme:

\[ \begin{bmatrix} \phi_1 (t) \\ \phi_2 (t) \\ \phi_3 (t) \\ \phi_4 (t) \end{bmatrix}_{(m+1)} = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} + \int_0^t \begin{bmatrix} \phi_2 (t) \\ -\phi_3 (t) \\ \phi_2 (t)\,\phi_4 (t) \\ -\phi_2 (t)\,\phi_3 (t) \end{bmatrix}_{(m)} {\text d} s , \qquad m=0,1,2,\ldots . \tag{3.4} \]

The first term is

\[ \begin{bmatrix} \phi_1 (t) \\ \phi_2 (t) \\ \phi_3 (t) \\ \phi_4 (t) \end{bmatrix}_{(1)} = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} + \int_0^t \begin{bmatrix} \frac{1}{10} \\ -\frac{1}{\sqrt{2}} \\ \frac{1}{10} \,\frac{1}{\sqrt{2}} \\ -\frac{1}{10} \,\frac{1}{\sqrt{2}} \end{bmatrix} {\text d} s = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} + t \begin{bmatrix} \frac{1}{10} \\ -\frac{1}{\sqrt{2}} \\ \frac{1}{10} \,\frac{1}{\sqrt{2}} \\ -\frac{1}{10} \,\frac{1}{\sqrt{2}} \end{bmatrix} . \]

Next,

\[ \begin{bmatrix} \phi_1 (t) \\ \phi_2 (t) \\ \phi_3 (t) \\ \phi_4 (t) \end{bmatrix}_{(2)} = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} + \int_0^t \begin{bmatrix} \frac{1}{10} - \frac{s}{\sqrt{2}} \\ -\frac{1}{\sqrt{2}} + \frac{s}{10\,\sqrt{2}} \\ \left( \frac{-s}{\sqrt{2}} + \frac{1}{10} \right) \left( \frac{1}{\sqrt{2}} + \frac{s}{10\,\sqrt{2}} \right) \\ -\left( \frac{-s}{\sqrt{2}} + \frac{1}{10} \right) \left( \frac{1}{\sqrt{2}} - \frac{s}{10\,\sqrt{2}} \right) \end{bmatrix} {\text d} s = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} + \begin{bmatrix} \frac{t}{10} - \frac{t^2}{2\,\sqrt{2}} \\ -\frac{t}{\sqrt{2}} + \frac{t^2}{20\,\sqrt{2}} \\ \frac{t}{10\,\sqrt{2}} - \frac{t^2}{4} + \frac{t^2}{200\,\sqrt{2}} - \frac{t^3}{60} \\ -\frac{t}{10\,\sqrt{2}} + \frac{t^2}{4} - \frac{t^2}{200\,\sqrt{2}} - \frac{t^3}{60} \end{bmatrix} . \]

We can expand this procedure indefinitely and obtain

\[ \theta (t) = \frac{\pi}{4} + \frac{t}{10} -\frac{\sqrt{2}}{4}\, t^2 - \frac{\sqrt{2}}{120}\, t^3 + \left( \frac{1}{48} + \frac{\sqrt{2}}{4800} \right) t^4 + \left( \frac{\sqrt{2}}{240 000} - \frac{1}{1200} \right) t^5 + \left( \frac{1111\,\sqrt{2}}{1600 000} - \frac{11}{144 000} \right) t^6 + \cdots . \]

We check our answer with Mathematica:

AsymptoticDSolveValue[{theta''[t] + Sin[theta[t]] == 0, theta[0] == Pi/4, theta'[0] == 1/10}, theta[t], {t, 0, 7}]

\[Pi]/4 + t/10 - t^2/(2 Sqrt[2]) - t^3/( 60 Sqrt[2]) + ((100 + Sqrt[2]) t^4)/4800 + ((-200 + Sqrt[ 2]) t^5)/240000 + ( 11 (-100 + 909 Sqrt[2]) t^6)/14400000 + ((400 + 159999 Sqrt[2]) t^7)/1008000000

		Now we plot Picard's polynomial approximation along with a true solution. theta6[t_] = \[Pi]/4 + t/10 - t^2/(2 Sqrt[2]) - t^3/( 60 Sqrt[2]) + ((100 + Sqrt[2]) t^4)/4800 + ((-200 + Sqrt[2]) t^5)/ 240000 + (11 (-100 + 909 Sqrt[2]) t^6)/14400000; s = NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == Pi/4, x'[0] == 1/10}, x, {t, 0, 8}]; Plot[{Callout[theta6[t], "true", LabelStyle -> Medium], Callout[x[t] /. s, "Picard", LabelStyle -> Medium]}, {t, 0, 5.0}, PlotTheme -> "Business"]
Figure 1: Picard's approximation along with the true solution.		Mathematica code

■

Example 4: Let us apply Picard's iteration procedure to a constant coefficient homogeneous vector differential equation

\[ \dot{\bf x} (t) = {\bf A}\, {\bf x}(t) , \qquad {\bf x}(0) = {\bf x}_0 , \tag{4.1} \]

where A is n × n constant matrix and x(t) is unknown n-column vector. First, we transfer the given initial value problem to the integral equation

\[ {\bf x}(t) = {\bf x}_0 + \int_0^t {\bf A}\, {\bf x}(s)\,{\text d}s = {\bf x}_0 + {\bf A} \int_0^t {\bf x}(s)\,{\text d}s . \tag{4.2} \]

Note that matrix A is constant, so we can exchange the order of integration and multiplication by matrix A. Now we employ Picard's iteration scheme:

\[ {\bf \phi}_{m+1} (t) = {\bf x}_0 + {\bf A} \int_0^t {\bf \phi}_m (s)\,{\text d}ds . \tag{4.3} \]

To start iteration, we have to choose the initial value \( {\bf \phi}_0 = {\bf x}_0 . \) Otherwise, the integral equation will be not equivalent to the given initial value problem. Then next iteration will be

\[ {\bf \phi}_1 = {\bf x}_0 + \int_0^t {\bf A}\, {\bf x}_0 \,{\text d}s = {\bf x}_0 + {\bf A}\,\int_0^t {\bf x}_0 \,{\text d}s = {\bf x}_0 + {\bf A}\,t \, {\bf x}_0 = {\bf x}_0 + t\,{\bf A} \, {\bf x}_0 . \]

The second term will be

\[ {\bf \phi}_2 = {\bf x}_0 + \int_0^t {\bf A}\, {\bf x}_1 \,{\text d}s = {\bf x}_0 + {\bf A}\,\int_0^t \left( {\bf x}_0 + {\bf A}\,s \right) {\text d}s = {\bf x}_0 + t\, {\bf A}\,{\bf x}_0 + \frac{t^2}{2} \, {\bf A}^2 {\bf x}_0 . \]

The third term becomes

\[ {\bf \phi}_3 = {\bf x}_0 + \int_0^t {\bf A}\, {\bf x}_2 \,{\text d}s = {\bf x}_0 + {\bf A}\,\int_0^t \left( {\bf x}_0 + {\bf A}\,s + {\bf A}^2 \frac{s^2}{2} \right) {\text d}s = {\bf x}_0 + t\, {\bf A}\,{\bf x}_0 + \frac{t^2}{2} \, {\bf A}^2 {\bf x}_0 + \frac{t^3}{2\cdot 3} \, {\bf A}^3 {\bf x}_0 . \]

The pattern is obvious and we can write the n-th iterate as

\[ {\bf \phi}_n = {\bf x}_0 + t\, {\bf A}\,{\bf x}_0 + \frac{t^2}{2} \, {\bf A}^2 {\bf x}_0 + \cdots + \frac{t^n}{n!} \, {\bf A}^n {\bf x}_0 . \tag{4.4} \]

E. Picard proved that the sequence of n-column vector functions ϕ_m(t) converges to the solution x>(t) of the given IVP as m → ∞. From the explicit formula for ϕ_m(t), it follows that it is the m-th approximation of the exponential matrix function:

\[ e^{{\bf A}\,t} = {\bf I} + t\,{\bf A} + \frac{t^2 {\bf A}^2}{2!} + \cdots + \frac{t^m {\bf A}^m}{m!} + \cdots . \tag{4.5} \]

Since matrix A is constant, the series for the exponential matrix converges for every real t. ■

Example 5: . Consider the system of differential equations subject to the initial conditions

\[ \begin{split} \dot{x} &= -6\,x + 4\,y , \\ \dot{y} &= \phantom{-}5\,x -5\,y ; \end{split} \qquad x(0) = 2, \quad y(0) =0 . \tag{5.1} \]

We convert the given IVP into the integral equation:

\[ \begin{bmatrix} x(t) \\ y(t) \end{bmatrix} = \begin{bmatrix} 2 \\ 0 \end{bmatrix} + \int_0^t \begin{bmatrix} -6\,x(s) + 4\,y(s) \\ 5\,x(s) -5\,y(s) \end{bmatrix} {\text d} s . \tag{5.2} \]

It is useful, but not necessary to rewrite the system of ODEs (5.1) as a single second order differential equation

\[ \ddot{z} + 11\,\dot{z} + 10\, z =0, \qquad z(0) =2, \quad \dot{z} (0) =-22 , \tag{5.3} \]

where z = x - y. We can express variables x and y through z:

\[ \dot{z} = \dot{x} - \dot{y} = -6 \,x + 4\,y - 5\,x + 5\,y = -11\,x + 9\,y \]

upon solving the system of algebraic equations

\[ \begin{split} z &= x-y , \\ \dot{z} &= -11\,x + 9\,y . \end{split} \]

Python helps:

This yields

\[ x = - \frac{1}{2} \left( 9\, z + \dot{z} \right) \qquad\mbox{and}\qquad y = - \frac{1}{2} \left( 11\, z + \dot{z} \right) . \]

Picard's iteration for the system becomes

\[ \begin{bmatrix} x(t) \\ y(t) \end{bmatrix}_{(m+1)} = \begin{bmatrix} 2 \\ 0 \end{bmatrix} + \int_0^t \begin{bmatrix} -6\,x(s) + 4\,y(s) \\ 5\,x(s) - 5\,y(s) \end{bmatrix}_{(m)} \,{\text d}s , \qquad m=0,1,2,\ldots . \tag{5.4} \]

Calculating the first few iterations, we obtain

\[ \begin{split} \begin{bmatrix} x (t) \\ y(t) \end{bmatrix}_{(1)} &= \begin{bmatrix} 2 \\ 0 \end{bmatrix} + \int_0^t \begin{bmatrix} -6\,x(0) + 4\,y(0) \\ 5\,x(0) - 5\,y(0) \end{bmatrix} = \begin{bmatrix} 2 - 12\,t \\ 10\,t \end{bmatrix} , \\ \begin{bmatrix} x (t) \\ y(t) \end{bmatrix}_{(2)} &= \begin{bmatrix} 2 \\ 0 \end{bmatrix} + \int_0^t \begin{bmatrix} -6\,x_1 (s) + 4\,y_1 (s) \\ 5\,x_1 (s) - 5\,y_1 (s) \end{bmatrix} = \begin{bmatrix} 2 - 12\,t + 56\, t^2 \\ 10\,t - 55\, t^2 \end{bmatrix} , \\ \begin{bmatrix} x (t) \\ y(t) \end{bmatrix}_{(3)} &= \begin{bmatrix} 2 - 12\,t + 56\, t^2 - \frac{556}{3}\, t^3 \\ 10\,t - 55\, t^2 + 185\, t^3 \end{bmatrix} . \end{split} \]

Using the single differential equation, we have

\[ \phi_{m+1} (t) = 2 -22t - \int_0^t \left( t-s \right) \left[ 11\,\phi'_m (s) + 10\,\phi_m (s) \right] {\text d} s, \qquad \phi_0 = 2, \quad \dot{\phi}_0 = -22 , \qquad m=0,1,2,\ldots . \]

Again, the first few iterations are

\[ \begin{split} \phi_1 (t) &= 2 - 22\,t + 111\,t^2 \qquad \Longrightarrow \qquad \dot{\phi}_1 = -22 + 222\, t , \\ \phi_2 (t) &= 2 - 22\,t + 111\,t^2 - \frac{1111}{3}\, t^3 - \frac{185}{2}\, t^4 \qquad \Longrightarrow \qquad \dot{\phi}_1 = -22 + 222\, t - 1111\, t^2 - 370\, t^3 , \\ \phi_3 (t) &= 2 - 22\,t + 111\,t^2 - \frac{1111}{3}\, t^3 - \frac{185}{2}\, t^4 \qquad \Longrightarrow \qquad \dot{\phi}_1 = -22 + 222\, t - 1111\, t^2 - 370\, t^3 . \end{split} \]

As we see, single equation iteration is faster than the previous. We can compare our single equation iteration with Picard's scheme for the original system of first order equations:

\[ x_1 = - \frac{1}{2} \left( 9\, \phi_1 + \dot{\phi}_1 \right) = -2 + 12\,t + \frac{999}{2}\, t^2 \]

Based on the single differential equation \( \ddot{z} + 11\,\dot{z} + 10\, z =0, \) we transfer it to the system of first order equations using standard substitutions:

\[ \dot{z} = u, \qquad \dot{u} = -11\, u -10\, z . \]

This leads to the following Picard's scheme:

\[ \begin{bmatrix} z(t) \\ u(t) \end{bmatrix}_{(m+1)} = \begin{bmatrix} 2 \\ -22 \end{bmatrix} + \int_0^t \begin{bmatrix} u(s) \\ -11\,u(s) - 10\,z(s) \end{bmatrix}_{(m)} {\text d}s , \qquad m=0,1,2,\ldots . \]

A first few terms are

\[ \begin{split} \begin{bmatrix} z(t) \\ u(t) \end{bmatrix}_{(1)} &= \begin{bmatrix} 2 - 22\,t \\ -22 + 222\,t \end{bmatrix} , \\ \begin{bmatrix} z(t) \\ u(t) \end{bmatrix}_{(2)} &= \begin{bmatrix} 2 - 22\,t + 111\,t^2 \\ -22 + 222\,t - 1111\,t^2 \end{bmatrix} , \\ \begin{bmatrix} z(t) \\ u(t) \end{bmatrix}_{(3)} &= \begin{bmatrix} 2 - 22\,t + 111\,t^2 - \frac{1111}{3}\,t^3 \\ -22 + 222\,t - 1111\,t^2 - \frac{11111}{3}\, t^3 \end{bmatrix} . \end{split} \]

■

Example 6: Consider the initial value problem:

\[ \frac{{\text d}^2 y}{{\text d}x^2} + 2x\, \frac{{\text d}y}{{\text d}x} - 2\,y(x) =0 , \qquad y(0) =1, \quad y' (0) = -1. \tag{6.1} \]

Its solution according to Mathematica is

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

\[ y(x) = e^{-x^2 /2} -x + x\, \sqrt{\pi} \,\mbox{erf}(x) = e^{-t^2 /2} -x + 2\,x\,\int_0^x e^{-t^2 /2} \,{\text d}t , \tag{6.2} \]

where erf(x) = \( \frac{2}{\sqrt{\pi}} \, \int_0^x e^{-t^2 /2} \,{\text d}t \) is the error function. Before applying Picard's iteration scheme, we find its Maclaurin approximation:

\[ y(x) = 1 - x + x^2 - \frac{x^4}{6} + \frac{x^6}{30} - \frac{x^8}{168} + \frac{x^{10}}{1080} - \frac{x^{12}}{7920} + \cdots . \tag{6.3} \]

We need this expansion to check the quality of Picard's iterations.

The given differential equation has a Maclaurin power series solution that converges everywhere in ℂ because its coefficients are holomorphic functions:

\[ y(x) = 1 - x + \sum_{n\ge 2} c_n x^n . \]

The linear term 1−x is chosen to satisfy the initial conditions, and coefficients c_n are to be determined from the equation. Its derivatives are

\[ y'(x) = -1 + \sum_{n\ge 2} c_n n\, x^{n-1} , \qquad y'' (x) = \sum_{n\ge 2} c_n n \left( n-1 \right) x^{n-2} = \sum_{k\ge 0} c_{k+2} \left( k+2 \right) \left( k+1 \right) x^k . \]

Substitution into the given differential equation yields

\begin{align*} y'' + 2x\,y' -2\,y &= \sum_{k\ge 0} c_{k+2} \left( k+2 \right) \left( k+1 \right) x^k -2x + 2 \sum_{n\ge 2} c_n n\, x^{n} -2 + 2x - 2\sum_{n\ge 2} c_n x^n \\ &= 2 \left( c_2 -1 \right) + 6x \,c_3 + \sum_{n\ge 2} \left[ c_{n+2} \left( n+2 \right) \left( n+1 \right) + 2c_n n - 2\,c_n\right] x^n = 0. \end{align*}

Equating coefficients of like powers of x, we obtain the recurrence:

\[ c_{n+2} = \frac{2 - 2n}{\left( n+2 \right) \left( n+1 \right)} c_n = - \frac{2\left( n-1 \right)}{\left( n+2 \right) \left( n+1 \right)} \, c_n , \qquad n=2,3,\ldots ; \]

where the first two terms are

\[ c_2 = 1 , \qquad c_3 = 0. \]

Thus, we obtain the convergent power series solution

\[ y(x) = 1 - x + x^2 - \frac{x^4}{6} + \frac{x^5}{30} - \frac{x^8}{168} + \cdots . \]

▣

First, we transfer the given problem to the system of first order differential equations by introducing new variables: y₁ = y(x) and y₂ to be its derivative. Then we get

\[ \frac{\text d}{{\text d}x}\,{\bf y} \equiv \frac{\text d}{{\text d}x} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} y_2 \\ 2\,y_1 - 2x\, y_2 \end{bmatrix} \qquad \Longleftrightarrow \qquad \begin{split} y'_1 &= y_2 , \\ y'_2 &= 2\, y_1 -2x\,y_2 . \end{split} \tag{6.4} \]

Starting with the initial condition as the first approximation, we obtain the Picard's iteration formula:

\[ \phi_{m+1} (x) \equiv \begin{bmatrix} y_1 (x) \\ y_2 (x) \end{bmatrix}_{(m+1)} = \begin{bmatrix} \phantom{-}1 \\ -1 \end{bmatrix} + \int_0^x \begin{bmatrix} y_2 (t) \\ 2\, y_1 (t) -2t\,y_2 (t) \end{bmatrix}_{(m)} {\text d}t , \qquad \phi_0 = \begin{bmatrix} \phantom{-}1 \\ -1 \end{bmatrix} . \tag{6.5} \]

The third approximation ϕ₃ has noise terms of order four. Of course, they will be eliminated with next iterations, but new ones of higher degree appear with every iterative step. This is typical for Picard's iterations because the given differential equation has variable coefficients. Noise terms are not observed only for autonomous differential equations.

Now we calculate some first iterations directly from the second order differential equation without converting it to a system of first order equations. Integrating twice both sides of the given differential equation, we obtain

\[ y(x) = 1 - x + \int_0^t {\text d}t \int_0^t \left[ - 2y' (s) + 2y(s) \right] {\text d}s . \]

Upon exchanging the order of integration, we get

\[ \int_0^t {\text d}t \int_0^t \left[ - 2y' (s) + 2y(s) \right] {\text d}s = \int_0^t \left[ - 2y' (s) + 2y(s) \right] {\text d}s \int_s^t {\text d}t = \int_0^t \left( t-s \right)\left[ - 2y' (s) + 2y(s) \right] {\text d}s . \]

The corresponding Picard's scheme becomes

\[ \phi_{(m+1)} (x) = 1-x + \int_0^x \left( x- s \right) \left[ 2\,\phi_{(m)} (s) - 2s\,\phi'_{(m)} (s) \right] {\text d}s , \qquad m=0,1,2,\ldots . \]

We can simplify this scheme further by eliminating the derivative of the unknown function under the sign of integral. Integrating by parts the term containing the derivative of ϕ_m, we reduce Picard's scheme to the following

\[ \phi_{(m+1)} (x) = 1-x + 2 \int_0^x \left( 2x- 3s \right) \phi_{(m)} (s)\, {\text d}s , \qquad m=0,1,2,\ldots . \tag{6.6} \]

For[{n = 1, y[0][x_] = {1, -1}}, n < 4, n++, y[n][x_] = {1, -1} + Integrate[{y[n - 1][t][[2]], 2*y[n - 1][t][[1]] - 2*t*y[n - 1][t][[2]]}, {t, 0, x}]];
Do[Print["n= ", n, " y(n)= ", y[n][x]], {n, 0, 3}]

n= 0 y(n)= {1,-1}
n= 1 y(n)= {1-x,-1+2 x+x^2}
n= 2 y(n)= {1-x+x^2+x^3/3,-1+2 x-(4 x^3)/3-x^4/2}
n= 3 y(n)= {1-x+x^2-x^4/3-x^5/10,-1+2 x-(2 x^3)/3+x^4/6+(8 x^5)/15+x^6/6}

So we get the third approximation to be

\[ \phi_3 (x) = 1 - x + x^2 - \frac{x^4}{3} - \frac{x^5}{10} , \]

which we compare with the true solution, obtained with the aid of Maclaurin expansion

\[ y(x) = 1 - x + x^2 - \frac{x^4}{6} + \frac{x^6}{30} + \cdots . \]

If you do not trust our iteration code, you can chack manually step by step. Upon choosing the initial term in accordance to the initial conditions

\[ \phi_{(0)} (x) = 1-x , \]

we find the first iterate:

\[ \phi_{(1)} (x) = 1-x + 2 \int_0^x \left( 2x- 3s \right) \left( 1 - s \right) \, {\text d}s = 1-x+x^2 \]

The next iteration gives

\[ \phi_{(2)} (x) = 1-x + 2 \int_0^x \left( 2x- 3s \right) \left( 1 - s - s^2 \right) \, {\text d}s = 1-x+x^2 - \frac{x^4}{6} , \]

and so on.

■

Example 7: . Consider the initial value problem for the homogeneous Duffing equation:

\[ \ddot{y} + 2\eta\,\dot{y} + y + \varepsilon \, y^3 = 0 , \qquad y(0) = d, \quad \dot{y} (0) = v , \tag{7.1} \]

where d = 1 is the initial displacement and v =-1 is the initial velocity. For simplicity, we choose the following numerical values of parameters: η = 3 and ε = 1. First, we convert the given problem to a system of first order differential equations

\[ \frac{\text d}{{\text d}t} \begin{bmatrix} y_1 (t) \\ y_2 (t) \end{bmatrix} = \begin{bmatrix} y_2 \\ -2\eta\, y_2 - y_1 - \varepsilon\, y_1^3 \end{bmatrix} = \begin{bmatrix} y_2 \\ -6\, y_2 - y_1 - y_1^3 \end{bmatrix} , \qquad \begin{bmatrix} y_1 (0) \\ y_2 (0) \end{bmatrix} = \begin{bmatrix} \phantom{-}1 \\ -1 \end{bmatrix} . \tag{7.2} \]

Upon integration, we transfer this initial value problem into the Volterra integral equation:

\[ \begin{bmatrix} y_1 (t) \\ y_2 (t) \end{bmatrix} = \begin{bmatrix} \phantom{-}1 \\ -1 \end{bmatrix} + \int_0^t \begin{bmatrix} y_2 (s) \\ -6\, y_2 (s) - y_1 (s) - y_1^3 (s) \end{bmatrix} {\text d}s . \tag{7.3} \]

This allows us to set up the Picard iteration scheme explicitly:

\[ {\bf \phi}_0 = \begin{bmatrix} d \\ v \end{bmatrix} = \begin{bmatrix} \phantom{-}1 \\ -1 \end{bmatrix} , \qquad {\bf \phi}_{n+1} (t) = \begin{bmatrix} x_{n+1} (t) \\ y_{n+1} (t) \end{bmatrix} = \begin{bmatrix} d\\ v \end{bmatrix} + \int_0^t \, \begin{bmatrix} y_n (s) \\ - x_n (s) - \varepsilon\, x^3_n (s) - 2\eta\, y_n (s) \end{bmatrix}_n \, {\text d} s , \qquad n=0,1,2,\ldots . \tag{7.4} \]

We show each iteration separately.

\begin{align*} {\bf \phi}_1 &= \begin{bmatrix} \phantom{-}1 \\ -1 \end{bmatrix} + \int_0^t \, \begin{bmatrix} -1 \\ 6-1-1 \end{bmatrix} \, {\text d} s = \begin{bmatrix} 1-t \\ 4t-1 \end{bmatrix} , \\ {\bf \phi}_2 &= \begin{bmatrix} \phantom{-}1 \\ -1 \end{bmatrix} + \int_0^t \, \begin{bmatrix} 4s-1 \\ 6(1-4s) -(1-s) - (1-s)^3 \end{bmatrix} \, {\text d} s \\ &= \begin{bmatrix} \phantom{-}1 \\ -1 \end{bmatrix} + \int_0^t \begin{bmatrix} 4s-1 \\ \left[ 4 -20s -3 s^2 + s^3 \right] \end{bmatrix} \, {\text d} s = \begin{bmatrix} 1-t + 2t^2 \\ -1+ 4t - 10 t^2 - t^3 + \frac{t^4}{4} \end{bmatrix} , \\ {\bf \phi}_3 &= \begin{bmatrix} \phantom{-}1 \\ -1 \end{bmatrix} + \int_0^t \, \begin{bmatrix} -1+ 4s - 10 s^2 - s^3 + \frac{s^4}{4} \\ - \left( 1-s + 2s^2 \right) - \left( 1-s + 2s^2 \right)^3 - 6 \left( -1+ 4s - 10 s^2 - s^3 + \frac{s^4}{4} \right) \end{bmatrix} \, {\text d} s \\ &= \begin{bmatrix} 1-t + 2t^2 - \frac{10}{3}\, t^3 - \frac{t^4}{4} + \frac{t^5}{20} \\ -1 + 4t - 10 t^2 + \frac{49}{3}\, t^3 + \frac{19}{4}\, t^4 - \frac{39}{10}\, t^5 + 2\,t^6 - \frac{8}{7} \, t^7 \end{bmatrix} , \\ {\bf \phi}_4 &= \begin{bmatrix} 1 \\ -1 \end{bmatrix} + \int_0^t \, \begin{bmatrix} -1 + 4s - 10 s^2 + \frac{49}{3}\, s^3 + \frac{19}{4}\, s^4 - \frac{39}{10}\, s^5 + 2\,s^6 - \frac{8}{7} \, s^7 \\ 4- 20t + \cdots \end{bmatrix} {\text d} s \\ &= \begin{bmatrix} 1-t + 2t^2 - \frac{10}{3}\, t^3 + \frac{49}{12}\, t^4 + \frac{19}{20}\, t^5 - \frac{13}{20}\, t^6 + \frac{2}{7}\, t^7 - \frac{1}{7}\, t^8 \\ -1 + 4t - 10 t^2 + \frac{49}{3}\, t^3 - \frac{215}{12}\, t^4 - \frac{131}{10}\, t^5 + \frac{279}{20}\, t^6 - \frac{5357}{420} \, t^7 + \cdots \end{bmatrix} . \end{align*}

Of course, all these manipulations we delegate to Python. To check the quality of Picard's iteration scheme, we compare the expansions obtained with correct Maclaurin expansion:

\[ y(t) = 1-t + 3\,t^2 - \frac{19}{3}\,t^3 + \frac{41}{4}\, t^4 - \frac{833}{60}\,t^5 + \frac{511}{30}\,t^6 - \frac{12937}{630}\,t^7 + \cdots . \]


True solution of Duffing's equation and its 10-term approximation (in blue).		Python code

■

Example 8: Jacobi Elliptic functions.

The German mathematician Carl Gustav Jakob Jacobi (1804--1851) came from a Jewish family but he was given the French style name Jacques Simon at birth. His father, Simon Jacobi, was a banker and his family was prosperous. Carl was the second son of the family, the eldest being Moritz Jacobi who eventually became a famous physicist. In 1825 Carl obtained the degree of Doctor of Philosophy with a dissertation on the partial fraction decomposition of rational fractions defended before a commission led by Enno Dirksen. He followed immediately with his Habilitation and at the same time converted to Christianity. Now qualifying for teaching University classes, the 21-year-old Jacobi lectured in 1825/26 on the theory of curves and surfaces at the University of Berlin (Prussia).

In 1827 he became a professor and in 1829, a tenured professor of mathematics at Königsberg University, and held the chair until 1842. Jacobi was the first Jewish mathematician to be appointed professor at a German university. Jacobi suffered a breakdown from overwork in 1843. He then visited Italy for a few months to regain his health. On his return he moved to Berlin, where he lived as a royal pensioner until his death. One of Jacobi's greatest accomplishments was his theory of elliptic functions and their relation to the elliptic theta function.

Consider the initial value problem

\[ \dot{x} ( t ) = y\, z, \quad \dot{y} (t) = -x\,z, \quad \dot{z} (t) = -k^2 x\,y , \qquad x(0) =0, \quad y(0) =1, \quad z(0) =1, \tag{8.1} \]

where dots stand for derivatives with respect to t, and k denotes a positive constant \( \le 1 . \) The parameter k is known as the modulus; its complementary modulus is \( \kappa = \sqrt{1-k^2} . \) The solutions to the initial value problem above are usually called the basic Jacobi elliptic functions and denoted as

\[ x(t) = \mbox{sn} (t) , \qquad y(t) = \mbox{cn}(t) , \qquad z(t) = \mbox{dn} (t) . \]

These functions were introduced by Jacobi in 1829. They are found in the description of the motion of a pendulum (see also pendulum (mathematics)), as well as in the design of the electronic elliptic filters.

Using Picard's iteration we derive the sequence

\begin{eqnarray*} \tag{8.2} x_{(k+1)} (t) &=& \int_0^t y_{(k)} (s)\,z_{(k)} (s) \,{\text d}s , \\ y_{(k+1)} (t) &=& 1 - \int_0^t x_{(k)} (s)\,z_{(k)} (s) \,{\text d}s , \\ z_{(k+1)} (t) &=& 1 - k^2 \int_0^t x_{(k)} (s)\,y_{(k)} (s) \,{\text d}s , \end{eqnarray*}

where initially x₀ = 0, y₀ = 1, z₀ = 1. The first iteration becomes

\begin{eqnarray*} x_{(1)} (t) &=& \int_0^t y_{(0)} (s)\,z_{(0)} (s) \,{\text d}s = \int_0^t \,{\text d}s = t , \\ y_{(1)} (t) &=& 1 - \int_0^t x_{(0)} (s)\,z_{(0)} (s) \,{\text d}s =1, \\ z_{(1)} (t) &=& 1 - k^2 \int_0^t x_{(0)} (s)\,y_{(0)} (s) \,{\text d}s = 1 \end{eqnarray*}

because x₀ = 0 and y₀ = z₀ = 1.

Second iteration is

\begin{eqnarray*} x_{(2)} (t) &=& \int_0^t y_{(1)} (s)\,z_{(1)} (s) \,{\text d}s = \int_0^t \,{\text d}s = t , \\ y_{(2)} (t) &=& 1 - \int_0^t x_{(1)} (s)\,z_{(1)} (s) \,{\text d}s =1- \int_0^t s\, {\text d}s = 1 - \frac{t^2}{2} , \\ z_{(2)} (t) &=& 1 - k^2 \int_0^t x_{(1)} (s)\,y_{(1)} (s) \,{\text d}s = 1 - k^2 \int_0^t s\, {\text d}s = 1 - k^2 \frac{t^2}{2} . \end{eqnarray*}

If we proceed successively, we find that the coefficients of the various powers of t ultimately become stable, i.e., they remain unchanged in successive iterations. This leads to the representation of solutions as convergent power series, which are usually denoted by the symbols snt (or sn(t.m) for m=k²), cnt, and dnt, and are called sine amplitude, cosine amplitude, and delta amplitude, respectively. These functions were introduced by Carl Jacobi in 1829.

\begin{eqnarray*} \mbox{sn} (t,k) &=& t - \left( 1 + k^2 \right) \frac{t^3}{3!} + \left( 1 + 14 k^2 + k^4 \right) \frac{t^5}{5!} - \left( 1 + 135 k^2 + 135 t^4 + k^6 \right) \frac{t^7}{7!} + O\left( t^9 \right) , \\ \mbox{cn} (t,k) &=& 1 - \frac{t^2}{2!} + \left( 1+ 4k^2 \right) \frac{t^4}{4!} - \left( 1 + 44k^2 + 16 k^4 \right) \frac{t^6}{6!} + O\left( t^8 \right) , \\ \mbox{dn} (t,k) &=& 1 - k^2 \frac{t^2}{2!} + k^2 \left( 4 + k^2 \right) \frac{t^4}{4!} - k^2 \left( 16 + 44k^2 + k^4 \right) \frac{t^6}{6!} + O\left( t^8 \right) . \end{eqnarray*}

Each of Jacobi elliptic functions satisfy a single-variable equation:

sn(t) solves the differential equations
\[ \frac{{\text d}^2 y}{{\text d}t^2} + \left( 1 + k^2 \right) y - 2\,k^2 y^3 =0 \qquad\mbox{and} \qquad \left( \frac{{\text d} y}{{\text d}t} \right)^2 = \left( 1- y^2 \right) \left( 1 - k^2 y^2 \right) . \]
cn(t) solves the differential equations
\[ \frac{{\text d}^2 y}{{\text d}t^2} + \left( 1 -2\, k^2 \right) y + 2\,k^2 y^3 =0 \qquad\mbox{and} \qquad \left( \frac{{\text d} y}{{\text d}t} \right)^2 = \left( 1- y^2 \right) \left( 1 - k^2 + k^2 y^2 \right) . \]
dn(t) solves the differential equations
\[ \frac{{\text d}^2 y}{{\text d}t^2} - \left( 2 - k^2 \right) y + 2\, y^3 =0 \qquad\mbox{and} \qquad \left( \frac{{\text d} y}{{\text d}t} \right)^2 = \left( y^2 -1 \right) \left( 1 - k^2 - y^2 \right) . \]

The functions satisfy the two square relations

\begin{eqnarray*} \mbox{cn}^2 (t) + \mbox{sn}^2 (t) &=& 1 , \\ \mbox{dn}^2 (t) + k^2 \mbox{sn}^2 (t) &=& 1 . \end{eqnarray*}

These functions are doubly periodic generalizations of the trigonometric functions satisfying

\begin{eqnarray*} \mbox{sn} (t,0) &=& \sin t , \\ \mbox{cn} (t,0) &=& \cos t , \\ \mbox{dn} (t,0) &=& 1 . \end{eqnarray*}

Now we plot some of them:

■

Djang, G.-F., A modified method of iteration of the Picard type in the solution of differential equations, Journal of The Franklin Institute, 1948, Vol. 246, Issue 6, pp. 453--457. https://doi.org/10.1016/0016-0032(48)90261-0
Hu, J. and Li, W.-P., Theory of Ordinary Differential Equations, The Hong Kong University of Science and Technology, 2005.
Let's Learn, Nemo, YouTube.
Robin, W.A., Solving differential equations using modified Picard iteration, International Journal of Mathematical Education in Science and Technology, Volume 41, 2010, Issue 5, Pages 649--665; https://doi.org/10.1080/00207391003675182

Part 2: Picard Iterations