MATLAB TUTORIAL for the First Course: Picard Iterations

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 = x₀, so f(x₀) = y₀. 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 x₀ to x, we get the following integral equation

\[ 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, φ₂(x) is obtained from φ₁(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) = y₀. 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

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
Boyce, W.E. and DiPrima, R.C., Elementary Differential equations and Boundary Value Problems, 11th edition, Wiley, New York; section 2.8.
Dobrushkin, V.A., Applied Differential Equations. The Primary Course, CRC Press, Boca Raton, FL, 2015. Section 1.6.
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.
Goode, S.W. and Annin, S.W., Differential Equations and Linear Algebra, 4th edition, Prentice Hall, Upper Saddle River, NJ, 2016, Appendix 4.
Hille, E., Lectures on Ordinary Differential Equations, Addison-Wesley Pub. Co., Reading, MA, 1969, pp. 32-41.
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.
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:

Example: Consider the initial value problem for the linear equation

\[ y' = x \left( 3 - 2y \right) , \qquad y(0) =1. \]

Since the given equation is actually separable, we have

\[ \frac{{\text d}y}{3-2y} = \frac{{\text d}x}{x} \qquad \Longrightarrow \qquad \ln | 3 - 2y | = x^2 + c . \]

The its explicit solution follows immediately

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

which guarantees the validity interval to be the whole real axis.

Expanding the explicit solution into Maclaurin series, we obtain

\[ y(x) = 1 + \frac{x^2}{2} - \frac{x^4}{4} + \frac{x^6}{12} - \frac{x^8}{2\cdot 4!} + \frac{x^{10}}{2\cdot 5!} + \cdots . \]

Indeed, Mathematica confirms

Series[3/2 - 1/2 *Exp[-x^2] , {x, 0, 10}]

SeriesData[x, 0, {1, 0, Rational[1, 2], 0, Rational[-1, 4], 0, Rational[1, 12], 0, Rational[-1, 48], 0, Rational[1, 240]}, 0, 11, 1]

Now we perform Picard's iterations:

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

and so on. So we see that first terms stabilize with each iteration, and we hope to achieve the explicit solution (expansion). We check with Mathematica:

h = Function[{t}, 1 + Integrate[x*(3 - 2*#), {x, 0, t}]][x] &; NestList[h, 1, 3]

which confirms our conclusion! ■

Example: Consider the initial value problem for the Riccati equation

\[ y' = x^2 + y^2 , \qquad y(0) =1. \]

First, we find its series solution.

Clear[seriesDSolve];
seriesDSolve[ode_, y_, iter : {x_, x0_, n_}, ics_: {}] :=
Module[{dorder, coeff},
dorder = Max[0, Cases[ode, Derivative[m_][y][x] :> m, Infinity]];
(coeff = Fold[
Flatten[{#1,
Solve[#2 == 0 /. #1,
Union@Cases[#2 /. #1, Derivative[m_][y][x0] /; m >= dorder,
Infinity]]}] &,
ics,
Table[
SeriesCoefficient[ode /. Equal -> Subtract, {x, x0, k}],
{k, 0, Max[n - dorder, 0]}]
];
Series[y[x], iter] /. coeff) /; dorder > 0 ]

seriesDSolve[y'[x] == x^2 + (y[x])^2, y, {x, 0, 15}, {y[0] -> 1}]

1 + x + x^2 + (4 x^3)/3 + (7 x^4)/6 + (6 x^5)/5 + (37 x^6)/30 + ( 404 x^7)/315 + (369 x^8)/280 + (428 x^9)/315 + (1961 x^10)/1400 +
( 75092 x^11)/51975 + (1238759 x^12)/831600 + (9884 x^13)/6435 + ( 17121817 x^14)/10810800 + (115860952 x^15)/70945875 + SeriesData[ x, 0, {}, 0, 16, 1]

So we have

\begin{align*} y(x) &= 1 + x + x^2 + \frac{4}{3}\,x^3 + \frac{7}{6}\, x^4 + \frac{6}{5}\,x^5 + \frac{37}{30}\, x^6 + \frac{404}{315}\, x^7 + \frac{369}{280} \, x^8 + \frac{428}{315}\,x^9 + \frac{1961}{1400}\,x^{10} \\ & \quad + \frac{75092}{51975}\,x^{11} + \frac{1238759}{831600}\, x^{12} + \frac{9884}{6435}\,x^{13} + \frac{17121817}{10810800}\, x^{14} + \frac{115860952}{70945875}\, x^{15} + \cdots . \end{align*}

Incidentally, this procedure can be automated very easily with a computer algebra system. Here's how to do so with Mathematica:

For[{n = 1, y[0][x_] = 1}, n < 4, n++,
y[n][x_] = 1 + Integrate[y[n - 1][t]^2 + t^2, {t, 0, x}];
Print[{n, y[n][t]}]]

{1,1+t+t^3/3}
{2,1+t+t^2+(2 t^3)/3+t^4/6+(2 t^5)/15+t^7/63}
{3,1+t+t^2+(4 t^3)/3+(5 t^4)/6+(8 t^5)/15+(29 t^6)/90+(47 t^7)/315+
(41 t^8)/630+(299 t^9)/11340+(4 t^10)/525+(184 t^11)/51975+t^12/2268+(4 t^13)/12285+t^15/59535}

\begin{align*} \phi_0 &= 1, \\ \phi_1 (x) &= 1 + x + \frac{x^3}{3} , \\ \phi_2 (x) &= 1 + x + x^2 + \frac{2\,x^3}{3} + \frac{1}{6}\, x^4 + \frac{2}{15}\, x^5 + \frac{x^7}{63} , \\ \phi_3 (x) &= 1 + x + x^2 + \frac{4}{3}\, x^3 + \frac{5}{6}\, x^4 + \frac{8}{15}\, x^5 + \frac{29}{90}\, x^6 + \frac{47}{315}\, x^7 \\ &\quad + \frac{41}{630}\, x^8 + \frac{299}{11340}\, x^9 + \frac{4}{525}\, x^{10} + \frac{184}{51975}\, x^{11} + \frac{1}{2268}\,x^{12} + \frac{4}{12285}\,x^{13} + \frac{x^{15}}{59535} . \end{align*}

If you want intermediate results, too, here is a function that keeps them:

picardList[initialVector_, flow_, n_, var_] :=
Module[{time},
FoldList[
iterate[
initialVector, flow, #1, #2, time
] &,
initialVector,
Reverse[Range[n] - 1]] /. {Subscript[time, _] -> var} ]

picardList[1, flow, 4, x] // TableForm

Picard's iterations could be found using the following subroutine:

Clear[picardSeries, iterate];
iterate[initial_, flow_, psi_, n_, t_] :=
initial +
Expand@Apply[Subtract,
Assuming[{Subscript[t, n + 1] > 0},
Integrate[flow[Subscript[t, n + 1], psi],
Subscript[t, n + 1]]] /. {{Subscript[t, n + 1] ->
Subscript[t, n]}, {Subscript[t, n + 1] -> 0}}]

picardSeries[initialVector_, flow_, n_, var_] :=
Module[{time},
Fold[iterate[initialVector, flow, #1, #2, time] &, initialVector,
Reverse[Range[n] - 1]] /. {Subscript[time, 0] -> var}]

Then we apply this subroutine to our problem:

flow[t_, f_] := f^2 + t^2
picardSeries[1, flow, 5, x]

There is no way to find a pattern or estimate the radius of convergence of the above series. Now we try to find its explicit solution:

DSolve[{y'[x] == (y [x])^2 + x^2, y[0] == 1}, y[x], x]

{{y[x] -> (-x^2 BesselJ[-(3/4), x^2/2] Gamma[1/4] + x^2 BesselJ[-(5/4), x^2/2] Gamma[3/4] +
BesselJ[-(1/4), x^2/2] Gamma[3/4] - x^2 BesselJ[3/4, x^2/2] Gamma[3/ 4])/
(x (BesselJ[1/4, x^2/2] Gamma[1/4] - 2 BesselJ[-(1/4), x^2/2] Gamma[3/4]))}}

Therefore, its solution is expressed through Bessel functions:

\[ y(x) = \dfrac{-x^2 J_{-3/4} (x^2 /2)\,\Gamma (1/4) + x^2 J_{-5/4} (x^2 /2)\,\Gamma (3/4) + J_{-1/4} (x^2 /2)\,\Gamma (3/4) -x^2 J_{3/4} (x^2 /2)\,\Gamma (3/4)}{x \, J_{1/4} (x^2 /2)\,\Gamma (1/4) -2\,x\,J_{-1/4} (x^2 /2)\,\Gamma (3/4)} . \]

Obviously, the solution blows up at the points where the denominator is zero.

FindRoot[(BesselJ[1/4, x^2/2] Gamma[1/4] - 2 BesselJ[-(1/4), x^2/2] Gamma[3/4]) == 0, {x, 1}]

{x -> 0.969811}

So we expect that the solution exists in the interval (0,0.969811). Upon plotting the fourth Picard iteration and the explicit solution, we see that

y[x_] = (-x^2 BesselJ[-(3/4), x^2/2] Gamma[1/4] +
x^2 BesselJ[-(5/4), x^2/2] Gamma[3/4] +
BesselJ[-(1/4), x^2/2] Gamma[3/4] -
x^2 BesselJ[3/4, x^2/2] Gamma[3/ 4])/
(x (BesselJ[1/4, x^2/2] Gamma[1/4] - 2 BesselJ[-(1/4), x^2/2] Gamma[3/4]))
phi4[x_] =
1 + x + x^2 + (4 x^3)/3 + (7 x^4)/6 + (6 x^5)/5 + (107 x^6)/90 +
( 362 x^7)/315 + (2683 x^8)/2520 + (2689 x^9)/2835 + (92179 x^10)/ 113400 +
(6673 x^11)/9900 + (201503 x^12)/374220 + (612853 x^13)/ 1474200 +
(3774677 x^14)/12162150 + (127144841 x^15)/567567000
Plot[{y[x], phi4[x]}, {x, 0, 1}, PlotStyle -> Thick]

Picard 4th iteration and explicit solution

Plotted figure shows a good approximation on the interval about [0,0.7]. Picard's theorem guarantees the existence of solution on the interval $ [0, 1/\sqrt{2}) \approx (0,0.7) , $ which is far away from the root 0.969811, where the denominator is zero. ■

Example: Consider the initial value problem for the Riccati equation

\[ y' = x^2 - y^2 , \qquad y(0) =1. \]

First, we find its explicit solution

DSolve[{y'[x] == -(y [x])^2 + x^2, y[0] == 1}, y[x], x]

which gives a complex-valued answer. So we look for another option:

\[ y(x) = x\,\frac{I_{-3/4} \left( x^2 /2 \right) \pi \left( \sqrt{2} + \pi / \Gamma^2 (3/4) \right) - 2 \,K_{3/4} \left( x^2 /2\right) }{\pi \left( \sqrt{2} + \pi / \Gamma^2 (3/4) \right) I_{1/4} \left( x^2 /2 \right) + 2\,K_{1/4} \left( x^2 /2 \right)} . \]

Actually, we need to know about the true solution only whether it exists for all x and its Maclaurin power series expansion:

\[ y(x) = 1-x + x^2 - \frac{2}{3}\,x^3 + \frac{5}{6}\,x^4 - \frac{4}{5}\, x^5 + \frac{23}{30}\,x^6 - \frac{236}{315}\, x^7 + \frac{201}{280}\, x^8 - \frac{218}{315}\,x^9 + \frac{2803}{4200}\,x^{10} - \frac{33412}{51975}\,x^{11} + \cdots . \]

Now we find Picard's iterations. Of course, we use first Mathematica:

For[{n = 1, y[0][x_] = 1}, n < 8, n++,
y[n][x_] = 1 + Integrate[-y[n - 1][t]^2 + t^2, {t, 0, x}];
Print[{n, y[n][t]}]]

\begin{align*} \phi_0 &= 1, \\ \phi_1 (x) &= 1 + \int_0^x \left( t^2 -1^2 \right) {\text d}t = 1 - x + \frac{x^3}{3} , \\ \phi_2 (x) &= 1 + \int_0^x \left( t^2 -\phi_1^2 (t) \right) {\text d}t = 1 -x + x^2 - \frac{x^4}{6} + \frac{2\,x^5}{15} - \frac{x^7}{63} , \\ \phi_3 (x) &= 1 + \int_0^x \left( t^2 -\phi_2^2 (t) \right) {\text d}t = 1 - x + x^2 - \frac{2x^3}{3} + \frac{x^4}{2} - \frac{2x^5}{15} - \frac{x^6}{10} + \cdots \\ \phi_4 (x) &= 1 + \int_0^x \left( t^2 -\phi_3^2 (t) \right) {\text d}t = 1 - x + x^2 - \frac{2x^3}{3} + \frac{5\,x^4}{6} - \frac{2\,x^5}{3} + \frac{13\,x^6}{30} - \cdots . \end{align*}

There are some discrepancies in coefficients; however, Picard's theorem guarantees that they will eventually go away, and the series is stabilized. ■

Example: Consider the initial value problem for the autonomous equation

\[ y' = -\sqrt{1/4 + y^2} , \qquad y(0) =3. \]

Its solution exists for all x and is given by

\[ y(x) = -\frac{1}{2}\,\sinh \left( x- \mbox{arcsinh}(6) \right) . \]

Mathematica confirms:

DSolve[{y'[x] == -Sqrt[1/4 + (y[x])^2], y[0] == 3}, y[x], x]

Out[1]= {{y[x] -> -(1/2) Sinh[x - ArcSinh[6]]}}

Next, we find its series solution (using previously defined subroutine seriesDSolve).

seriesDSolve[ y'[x] == Sqrt[1/4 + (y[x])^2], y, {x, 0, 15}, {y[0] -> 3}]

3 + (Sqrt[37] x)/2 + (3 x^2)/2 + (Sqrt[37] x^3)/12 + x^4/8 +
( Sqrt[37] x^5)/240 + x^6/240 + (Sqrt[37] x^7)/10080 + x^8/13440 +
( Sqrt[37] x^9)/725760 + x^10/1209600 +
( Sqrt[37] x^11)/79833600 + x^12/159667200 +
( Sqrt[37] x^13)/12454041600 + x^14/29059430400 +
( Sqrt[37] x^15)/2615348736000 + SeriesData[x, 0, {}, 0, 16, 1]

Now we try to use Picard's iterations:

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

Mathematica failed to find the second Picard iteration. However, Maple can handle this integral providing the answer:

\[ \phi_2 (x) = 3/2-{\frac {\sqrt {37}\ln \left( 6+\sqrt {37} \right) }{148}} + {\frac{3\,\sqrt {37+37\,{x}^{2}+12\,\sqrt {37}x}\sqrt {37}}{74}}+1/4\,\sqrt {37+37\,{x}^{2}+12\,\sqrt {37}x}x+{\frac {\sqrt {37}{\rm arcsinh} \left(1/37\,\sqrt {37} \left( 6\,\sqrt {37}+37\,x \right) \right)}{148 }} . \]

So we conclude that Picard's procedure is not practical because we cannot find the next iteration. ■

Example: Consider the initial value problem

\[ \dot{y} = \sin (y) + \sin (t) , \qquad y(0) = 1. \]

Since the slope function is not a polynomial in two variables y,t, we do not expect that Picard's iteration can be performed because the explicit integration is not possible. However, chebfun can be used that is based on Chebyshev approximation to the slope function.


d = 8; t = chebfun('t',[0 d]); u0 = 1;
L = chebop(0,d); L.op = @(t,u) diff(u) - sin(u); L.lbc = u0;
yexact = L\sin(t);




This first plot shows iterates k=0,…,4, with the exact solution in red.

 




   
        
    
u = u0 + 0*t;
f = @(u,t) sin(u) + sin(t);
LW = 'linewidth'; FS = 'fontsize'; IN = 'interpret'; LT = 'latex';
hold off
ss = @(k) ['$k = ' int2str(k) '$'];
for k = 0:4
  plot(u,'b',LW,1.6), hold on, ylim([-3 10])
  text(1.015*d,u(end),ss(k),IN,LT)
  u = u0 + cumsum(f(u,t));
end
plot(uexact,'r',LW,1.6), xlabel('t',FS,10), ylabel('u',FS,10)
title('Picard iterates $k = 0,\dots,4$',FS,12,IN,LT)
    
        







A second plot shows k=5,…,9.





   
        
    
hold off
for k = 5:9
  plot(u,'b',LW,1.6), hold on, ylim([0 7])
  text(1.015*d,u(end),ss(k),IN,LT)
  u = u0 + cumsum(f(u,t));
end
plot(uexact,'r',LW,1.6), xlabel('t',FS,10), ylabel('u',FS,10)
title('Picard iterates $k = 5,\dots,9$',FS,12,IN,LT)
    
        






A third plot shows k=10,…,14.






   
        
    
hold off
for k = 10:14
  plot(u,'b',LW,1.6), hold on, ylim([1 6])
  text(1.015*d,u(end),ss(k),IN,LT)
  u = u0 + cumsum(f(u,t));
end
plot(uexact,'r',LW,1.6), xlabel('t',FS,10), ylabel('u',FS,10)
title('Picard iterates $k = 10,\dots,14$',FS,12,IN,LT)
    
        








These plots show vividly the kind of convergence one can expect from a
Picard iteration: starting at the initial condition, sweeping slowly
across the domain. There is a numerical method based on this idea,
called waveform relaxation, but one can see immediately from the
pictures that it is unlikely to be efficient when carried out over
long time intervals. Instead, standard numerical methods just march
once rather than many times from left to right, but they march with a
small discrete time step and a discrete formula of higher order. 



To see the convergence in a quantitative fashion, it is interesting plot the errors of iterates 0,…,4 as a function of t on a loglog plot. The zeroth iterate has accuracy O(t), the first has accuracy O(t²), and so on:




   
        

    
u = u0 + 0*t;
ss = @(k) ['$k = ' int2str(k) '$'];
tt = logspace(-2,log10(8),600); hold off
for k = 0:4
  errtt = abs(u(tt)-uexact(tt));
  loglog(tt,errtt,'k',LW,.7), hold on
  text(8.7,errtt(1),ss(k),IN,LT)
  u = u0 + cumsum(f(u,t));
end
xlabel('t',FS,10), ylabel('error',FS,10)
axis([1e-2 8 1e-16 1e3])
title('Errors of iterates $0,\dots,4$',FS,12,IN,LT)
    
        













   ■










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




    






    
    
Example: 
**DESCRIPTION OF PROBLEM GOES HERE**
This is a description for some MATLAB code. MATLAB is an extremely useful tool for many different areas in engineering, applied mathematics, computer science, biology, chemistry, and so much more. It is quite amazing at handling matrices, but has lots of competition with other programs such as Mathematica and Maple. Here is a code snippet plotting two lines (y vs. x and z vs. x) on the same graph. Click to view the code!



    




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.
   ▣