Charles Picard

Consider the initial value problem

\[ \frac{{\text d}y}{{\text d}x} = f(x,y), \qquad y(x_0 ) = y_0 . \]
The derivative operator \( \texttt{D} = {\text d}/{\text d}x \) is an unbounded operator because it could transfer a bounded function into unbounded function. For example, bounded and smooth function sin(ωt) is mapped by the derivative operator into unbounded function ωcos(ωt), which may exceed any number upon choosing an appropriate value of ω. This example shows that a fast oscillating function may have a very large derivative, so its finite difference approximation will be problematic. The same affect is observed when the solution of a differential equation may be unbounded at some point. Therefore, a practical approximation of such functions may not be accurate unless we know that a solution is a smooth varying function. As an alternative, we try to reformulate the initial value problem as an integral equation.

The inverse of the derivative operator \( \texttt{D} \) is not an operator in general mathematical sense because it assigns to every function infinite many outputs that are usually called the «antiderivative». In order to single out one output, we consider the differential operator on a special set of functions that have a prescribed value at one point x = x0, so f(x0) = y0. We denote \( \texttt{D} \) on such set of functions (which is not a linear space) as L. If we rewrite our initial value problem in the operator form \( L[y] = f , \) then apply its inverse L-1 to both sides, we reduce the given problem to the fixed point problem

\[ y = L^{-1} \left[ f(x,y) \right] \]
for a bounded operator L-1. If the given initial value problem is homogeneous, then Banach's or Brouwer's fixed point theorem assures us that for some bounded operators (such as L-1f for Lipschitz slope function f) the iteration \( y_{n+1} = L^{-1} \left[ f(x,y_n ) \right] , \ y_0 = y(x_0 ) \) converges to a true solution. So Picard's iteration is actually an extension and implementation of this contraction theorem. Banach’s fixed point theorem formulates an easy to check assumption under which we have existence and uniqueness and computabilty of a fixed point is guaranteed. Moreover, it presents the order of the convergence of the approximating sequence and gives hints to stop the iteration infinite steps.

If we integrate (that is, apply the inverse operator L-1 to the unbounded derivative operator restricted on the set of functions with a specified initial condition) both sides of the differential equation \( y' = f(x,y) \) with respect to x from x0 to x, we get the following integral equation

Function sin(100t)
\[ y(x) = y(x_0 ) + \int_{x_0}^x f\left( t, y(t) \right) {\text d} t . \]
The above (Volterra) integral equation in the right-hand side is actually the inverse operator generated by the initial value problem. Therefore, both problems, the initial value problem and its integral reformulation, are equivalent in the sense that every solution of one of them is also a solution of its counterpart and vise versa. Although this integral equation is very hard (if not impossible) to solve, it has some advantages over the initial value problem. The main feature of the integral equation is that it contains integration---a bounded operator. This should be almost obvious for everyone who took an introductory calculus course: integration (of a positive function) is an area under the curve, which cannot be infinite for any smooth and bounded function on a finite interval.

Now we deal with the fixed point problem that is more computationally friendly because we can approximate an integral with any accuracy we want. The Picard iterative process consists of constructing a sequence of functions { φn } that will get closer and closer to the desired solution. This is how the process works:

\[ \phi_0 \equiv y_0 \qquad \mbox{for all $x$}. \]
For nonzero parameter n, the recurrent formula holds
\[ \phi_{n+1} (x) = y_0 + \int_{x_0}^x f\left( t, \phi_n (t) \right) {\text d} t, \qquad n=0,1,2,\ldots . \]
This recurrence is equivalent to the sequence of initial value problems:
\[ \frac{{\text d}\phi_{n+1} (x)}{{\text d}x} = f\left( x, \phi_n (x) \right) , \qquad \phi_n (x_0 ) = y_0 , \qquad n=0,1,2,\ldots . \]
Then the limit function \( \phi (x) = \lim_{n\to \infty} \phi_n (x) \) gives the solution of the original initial value problem. However, practical implementation of Picard's iteration is possible only when the slope function is a polynomial. Otherwise, explicit integration is most likely impossible to achieve, even for a very sophisticated computer algebra system. We will show in the second part how to bypass this obstacle.
Theorem: Suppose that the slope function f satisfies the existence and uniqueness theorem ( Picard--Lindelöf) in some rectangular domain \( \displaystyle R = [x_0 -a, x_0 +a]\times [y_0 -b , y_0 +b ] : \)
  • The slope function is bounded: \( M = \max_{(x,y) \in R} \, |f(x,y)| . \)
  • The slope function is Lipschitz continuous: \( |f(x,y) - f(x,z) | \le L\,|y-z| . \) for all (x,y) from R.
Then the successive approximation
\[ \phi_n (x) = y_0 + \int_{x_0}^x f \left( t, \phi_{n-1} (t) \right) {\text d}t , \qquad \phi_0 = y_0 , \]
converges to the unique solution of the initial value problem \( y' = f(x,y) , \ y(x_0 ) = y_0 , \) in the interval \( x_0 - h \le x \le x_0 +h , \) where \( \displaystyle h = \min \left\{ a, \frac{b}{M} \right\} . \)    ⧫

If we chose the initial approximation as \( \phi_0 (x_0 = y_0 , \) then at least the initial condition \( y( x_0 ) = y_0 \) is satisfied. The next approximation is obtained by substituting the initial approximation into the right hand side of integral equation:
\[ \phi_1 (x) = y_0 + \int_{x_0}^x f \left( t, \phi_{0} (t) \right) {\text d}t = y_0 + \int_{x_0}^x f \left( t, y_{0} (t) \right) {\text d}t . \]
Similarly, φ2(x) is obtained from φ1(x):
\[ \phi_2 (x) = y_0 + \int_{x_0}^x f \left( t, \phi_{1} (t) \right) {\text d}t . \]
In general, we have
\[ \phi_{n+1} (x) = y_0 + \int_{x_0}^x f \left( t, \phi_{n} (t) \right) {\text d}t . \]
In this manner, we generate the sequence of functions \( \displaystyle \left\{ \phi_n (x) \right\} = \left\{ \phi_0 , \phi_1 (x), \phi_2 (x) , \ldots \right\} . \) Each member of this sequence satisfies the initial condition, but in general none satisfies the differential equation. If for some index m, we find that \( \displaystyle \phi_{m+1} (x) = \phi_m (x) , \) then φm(x) is a fixed point of the integral equation and it is a solution of the differential equation as well. Hence, the sequence is terminated at this point.

First, we verify that all members of the sequence exist. This follows from the inequality:

\[ \left\vert \phi_{n+1} (x) \right\vert \le \left\vert y_0 \right\vert + \max \left\vert f(x,y) \right\vert \int_{x_0}^x {\text d}t = M\,|x-x_0 | \le b \]
for \( \displaystyle \left\vert x - x_0 \right\vert \le b\,M, \) where \( \displaystyle M = \max \left\vert f(x,y \right\vert . \)

To show that the sequence converges, we represent the limit function as infinite series

\[ \phi (x) = \lim_{n\to \infty} \phi_n (x) = \phi_0 + \sum_{n\ge 1} \left[ \phi_n (x) - \phi_{n-1} (x) \right] . \]
The convergence of the sequence \( \displaystyle \left\{ \phi_n (x) \right\} \) is established by showing that the above infinite series converges. To do this, we estimate the magnitude of the general term \( \displaystyle \left\vert \phi_n (x) - \phi_{n-1} (x) \right\vert \) .

We go step by step, and consider the first term

\[ \left\vert \phi_1 (x) - \phi_0 \right\vert = \left\vert \int_{x_0}^x f(t, y_0 ) \,{\text d}t \right\vert \le M \int_{x_0}^x f(t, y_0 ) \,{\text d}t = M \left\vert x - x_0 \right\vert . \]
For the second term, we have
\[ \left\vert \phi_2 (x) - \phi_1 (x) \right\vert \le \int_{x_0}^x \left\vert f(t, \phi_1 (t) ) - f(t, \phi_0 ) \right\vert {\text d}t \le L \int_{x_0}^x M \left( t - x_0 \right) {\text d}t = \frac{1}{2}\,ML\left( x- x_0 \right)^2 , \]
where L is the Lipschitz constant for the function f(x,y), that is, \( \displaystyle \left\vert f(x, y_1 ) - f(x, y_2 ) \right\vert \le L \left\vert y_1 - y_2 \right\vert . \) For the n-th term, we have
\begin{align*} \left\vert \phi_n (x) - \phi_{n-1} (x) \right\vert &\le \int_{x_0}^x \left\vert f \left( t , \phi_{n-1} (t) \right) - f \left( \phi_{n-2} (t) \right) \right\vert {\text d}t \\ &\le L \int_{x_0}^x \left\vert \phi_{n-1} (t) - \phi_{n-2} (t) \right\vert {\text d}t \le L \int_{x_0}^x \frac{M\,L^{n-2} \left( t- x_0 \right)^{n-1}}{(n-1)!} \, {\text d}t \\ &= \frac{M\,L^{n-1} \left\vert x - x_0 \right\vert^n}{n!} \le \frac{M\,L^{n-1} h^n}{n!} , \qquad h= \max \left\vert x - x_0 \right\vert . \end{align*}
Substituting these formulas into the finite sum
\[ \phi_n (x) = \phi_0 + \sum_{k=1}^n \left[ \phi_k (x) - \phi_{k-1} (x) \right] , \]
we obtain
\begin{align*} \left\vert \phi_n (x) - \phi_{n-1} (x) \right\vert &\le M \left\vert x - x_0 \right\vert + \frac{M\,L \left\vert x-x_0 \right\vert^2}{2} + \cdots + \frac{M\,L^{n-1} \left\vert x-x_0 \right\vert^n}{n!} . \end{align*}
When n approaches infinity, the sum
\[ M \left\vert x - x_0 \right\vert + \frac{M\,L \left\vert x-x_0 \right\vert^2}{2} + \cdots + \frac{M\,L^{n-1} \left\vert x-x_0 \right\vert^n}{n!} \]
approaches \( \displaystyle e^{L \left\vert x - x_0 \right\vert} -1 . \) Therefore, we have
\[ \left\vert \phi_n (x) - \phi_{n-1} (x) \right\vert \le \frac{M}{L} \left[ e^{L \left\vert x - x_0 \right\vert} -1 \right] \]
for any n. Hence, by the Weierstrass M-test, the series \( \displaystyle \phi (x) = \lim_{n\to \infty} \phi_n (x) = \phi_0 + \sum_{n\ge 1} \left[ \phi_n (x) - \phi_{n-1} (x) \right] \) converges absolutely and uniformly on the interval \( \displaystyle \left\vert x - x_0 \right\vert \le h . \) Then it follows that the limit function \( \displaystyle \phi (x) = \lim_{n\to \infty} \phi_n (x) \) of the sequence \( \displaystyle \left\{ \phi_n (x) \right\} \) is a continuous function on the interval \( \displaystyle \left\vert x - x_0 \right\vert \le h . \)

Next we prove that the limit function is a solution of the given initial value problem. Allowing n to approach infinity on both sides of

\[ \phi_{n+1} (x) = \phi_0 + \int_{x_0}^x f \left( t, \phi_n (t) \right) {\text d}t , \]
we get
\[ \phi (x) = \lim_{n\to\infty} \,\phi_{n+1} (x) = \phi_0 + \lim_{n \to \infty} \int_{x_0}^x f \left( t, \phi_n (t) \right) {\text d}t , \]
Since the slope function satisfies the Lipschitz condition
\[ \left\vert f \left( t, \phi_n (t) \right) - f \left( t, \phi_{n-1} (t) \right) \right\vert \le L \left\vert \phi_n (t) - \phi_{n-1} (t) \right\vert \]
on the interval \( \displaystyle \left\vert t - x_0 \right\vert \le h , \) it follows that the sequence \( \displaystyle f \left( t, \phi_n (t) \right) \) also converges uniformly to \( \displaystyle f \left( t, \phi (t) \right) \) on this interval. So we can interchange integration with the limit operation on the right-hand side to obtain
\[ \phi (x) = y_0 + \int_{x_0}^x \lim_{n\to +\infty} \, f \left( t, \phi_n (t) \right) {\text d}t = y_0 + \int_{x_0}^x f \left( t, \phi (t) \right) {\text d}t . \]
Thus, the limit function is a solution of the Volterra integral equation that is equivalent to the given initial value problem. Note that the limit function automatically satisfies the initial condition y(0) = y0. Also, differentiating both sides of the latter with respect to x and noting that that the right-hand side is a differentiable function of the upper limit, we find
\[ \phi' (x) = f(x, \phi (x)) . \]

 

More details for the existence and uniqueness of the initial value problem can be found in the following textbooks
  1. Birkhoff, G. and Rota, G.-C., Ordinary Differential equations, 4rd Ed., John Wiley and Sons, New York, NY, 1989, Chapter 1, p. 23. https://doi.org/10.1137/1005043
  2. Boyce, W.E. and DiPrima, R.C., Elementary Differential equations and Boundary Value Problems, 11th edition, Wiley, New York; section 2.8.
  3. Dobrushkin, V.A., Applied Differential Equations. The Primary Course, CRC Press, Boca Raton, FL, 2015. Section 1.6.
  4. Edwards, C.H. and Penney, D.E., Differential Equations and Boundary Value Problems: Computing and Modeling, 5th edition, 2014, Prentice Hall, Upper Saddle River, NJ, Appendix A.1.
  5. Goode, S.W. and Annin, S.W., Differential Equations and Linear Algebra, 4th edition, Prentice Hall, Upper Saddle River, NJ, 2016, Appendix 4.
  6. Hille, E., Lectures on Ordinary Differential Equations, Addison-Wesley Pub. Co., Reading, MA, 1969, pp. 32-41.
  7. Nagle, R.K., Saff, E.K., and Snider, A.D., Fundamentals of Differential Equations and Boundary Value Problems, 8th edition, Addison-Wesley, Boston, MA, 2019, Chapter 13, Sections 1 and 2.
  8. Wouk, A., On the Cauchy-Picard Method, The American Mathematical Monthly, Vol. 70, No. 2. (Feb., 1963), pp. 158-162.
   ▣

Although this approach is most often connected with the names of Charles Emile Picard, Giuseppe Peano, Ernst Lindelöf, Rudolph Lipschitz, and Augustin Cauchy, it was first published by the French mathematician Joseph Liouville in 1838 for solving second order homogeneous linear equations. About fifty years later, Picard developed a much more generalized form, which placed the concept on a more rigorous mathematical foundation, and used this approach to solve boundary value problems described by second order ordinary differential equations.

Note that Picard's iteration procedure, if it could be performed, provides an explicit solution to the initial value problem. This method is not for practical applications mostly for two reasons: finding next iterate may be impossible, and, when it works, Picard's iteration contains repetitions. Therefore, it could be used sometimes for approximations only. Moreover, the interval of existence (-h, h) for the considered initial value problem is in most cases far away from the validity interval. There are known some improvements in this direction (for instance, S.M. Lozinsky theorem), but we are far away from complete determination of the validity interval based on Picard's iteration procedure. Also, there are a number of counterexamples in which Picard iteration is guaranteed not to converge, even if the starting function is arbitrarily close to the real solution.

Here are some versions of Picard's iteration procedure for matlab:

================== to be removed ============

Two n-by-n matrices A and B are called similar if there exists an invertible n-by-n matrix S such that
Theorem: If λ is an eigenvalue of a square matrix A, then its algebraic multiplicity is at least as large as its geometric multiplicity.    ▣