Preface

---

This secton demonstrates applications of Picard's iteration scheme for power series determination of solutions to the initial value problems for second and higher order differential equations.

Picard iterations

---

Picard's theorem on the existence and uniqueness of solutions of first order differential equations was discussed previously in the introductory section. Let us outline the main ideas in Picard's iteration procedure, starting with the first order differential equation Here \( y' = {\text d}y/{\text d}x \) is the derivative of function y(x). It is convenient to use the derivative operator \( \texttt{D} = {\text d}/{\text d}x , \) and rewrite Eq.\eqref{EqPicard.1} in operator form The iteration procedure is very basic mathematical technique that conceptually sets the tone for a large class of analytical and numerical methods. We are going to apply this scheme for determination of series solutions for wide class of differential equations.

Since the differential equation contains the unbounded derivative operator, it is hard to expect that any iteration procedure applied to the differential equation \eqref{EqPicard.1} will lead to an acceptable solution. Therefore, the first step is to get rid of this unbounded operator by integration: \begin{equation} \label{EqPicard.2} y(x) = y\left( x_0 \right) + \int_{x_0}^x f\left( s,y(s) \right)\,{\text d} s . \end{equation} There are some advantages in converting the given initial value problem \eqref{EqPicard.1} into an equivalent problem without derivatives. You can recognize in the right-hand side integral the inverse to the derivative operator \( \texttt{D} = {\text d}/{\text d}x , \) acting in the space of smooth function vanishing at x₀. So adding y(x₀) is needed to satisfy the initial condition. Next, the form of Eq.\eqref{EqPicard.2} tells us that we can apply a fixed point theorem to the equation y = Φ(y) for a bounded integral operator \( \Phi (y) = \int_{x_0}^x f\left( s,y(s) \right)\,{\text d} s . \) Eq.\eqref{EqPicard.2} is known as the Volterra integral equation of the second kind, which is suitable for the iteration procedure: \begin{equation} \label{EqPicard.3} \phi_{m+1} (x) = y\left( x_0 \right) + \int_{x_0}^x f\left( s,\phi_m (s) \right)\,{\text d} s , \qquad m=0,1,2,\ldots , \qquad \phi_0 = y\left( x_0 \right) . \end{equation} When the slope function f(x, y) is a Lipschitz continuous function with respect to variable y, the Picard's iteration \eqref{EqPicard.3} converges uniformly to a unique solution of the initial value problem \eqref{EqPicard.1} on some interval containing the initial point x₀.

Picard's iteration scheme can be implemented in Mathematica in many ways. Some of them are presented below.

The iteration step is called iterate, and it keeps track of the iteration order n so that it can assign a separate integration variable name at each step. As additional tricks to speed things up, I avoid the automatic simplifications for definite integrals by doing the integral as an indefinite one first, then using Subtract to apply the integration limits.

Instead of a generic simplification, I add Expand to the result in order to help the subsequent integration step recognize how to split up the integral over the current result. The last argument in iterate is the name that I intend to use as the base variable name for the integration, but it will get a subscript n attached to it.

Then iterate is called by picardSeries, using increasing subscripts starting with the outermost integration. This is implemented using Fold with a Reverse range of indices. The first input is the initial condition, followed by the function defining the flow (specifically, defined below as flow[t_, f_]). After the order n, I also specify the name of the independent variable var as the last argument. Here is the code from the following example involving the Riccati equation:

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

Then the previous code will be as follows

This solution works equally well for vector initial value problems, i.e., the flow can be a vector function and the initial condition a vector. It can be extended for higher order differential equations.

Example: Riccati equation Example: Abel equation Since every iteration step in \eqref{EqPicard.3} requires integration, it can be efficiently performed only when the slope function is a polynomial. Therefore, the general effectiveness of Picard’s method fails apparently even for relatively simple problems such as the pendulum equation. This criticisms can be overcomed by applying a well-known transformation technique or polynomial approximations (as shown by L.N. Trefethen using Chebyshev polynomials).

 

In this subsection, we specifically consider second order ordinary differential equations (ODEs for short), i.e., equations of the form The theory of second order ODEs is an extensive and highly developed part of mathematics, not least because of its numerous applications in the physical sciences. For example, Newton told us, with impressive foresight, that the radial motion of the space shuttle is governed by the equation where M is the mass of the earth and G is the gravitational constant. As another sample, one may consider a simple pendulum equation None of these equations is suitable to practical Picard's iteration because they involve nonpolynomial functions. Therefore, they can be solved with Picard's method only upon transformations into polynomial forms---it will be shown in the second part of the tutorial.

Integration twice the differential equation \eqref{EqPicard.4}, we obtain

This equation is suitable for a fixed point theorem and the Picard's iteration yields Sometimes (at least in case of linear equations) it is possible to integrate by parts and eliminate the derivative from the sign of the integral in Eq.\eqref{EqPicard.5}. Next, we illustrate the iteration procedure on several examples. Example: Airy equation Example: stiff differential equation Example: Legendre's differential equation Example: Lane--Emden equation

 

Consider a n-th order (nonlinear) differential equation subject to the initial conditions: As usual, \( y^{(k)} = \texttt{D}^k y = {\text d}^k / {\text d} x^k \) is the k-th derivative of the function y(x), with \( y' = \texttt{D}\, y = {\text d} / {\text d} x \) and \( y'' = \texttt{D}^2 y = {\text d}^2 / {\text d} x^2 . \) It is assumed that the given initial value problem \eqref{EqPicard.n1} has a unique solution for given n constants α₀, α₁, … , α_n-1.

We integrate Eq.\eqref{EqPicard.n1} repeatedly, n times, to get

\begin{equation} \label{EqPicard.n2} y(x) = \frac{1}{(n-1)!} \int_{x_0}^x \left( x- t \right)^{n-1} f(t)\,{\text d}t + \sum_{k=0}^{n-1} \frac{\alpha_k}{k!}\, x^k . \end{equation}

Here and we call the initial polynomial because it satisfies the given initial conditions in \eqref{EqPicard.n1}. Equation \eqref{EqPicard.n2} is solved now by iteration, upon choosing a suitable starting function that satisfies the initial conditions and then substituting φ₀ into the right hand side of \eqref{EqPicard.n1}, i.e. \eqref{EqPicard.n2}, which can be integrated to yield ϕ₁. Repeating the iterative process, we get the general formula to obtain ϕ_m+1 from ϕ_m The iteration process \eqref{EqPicard.n3} is terminated when the required degree of accuracy is attained. Example: Blasius differential equation

 

References

 

Return to Mathematica page
Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)
Return to the Part 7 (Boundary Value Problems)