Picard Iteration

Differential Equations

The Picard iteration is an iterative procedure for solving differential equations. It is most often used in proofs of the existence of a differential equation's solution and is not very numerically performant. There are some approaches to remedy this which we will explore.

Introduction

Given a differential equation of the following form

\frac{d y}{d x} = f (x, y (x)) y (x_{0}) = y_{0}

We may write this in integral form

y (x) - y (x_{0}) = \int_{x_{0}}^{x} f (t, y (t)) d t

Rearranging the formula:

y (x) = y (x_{0}) + \int_{x_{0}}^{x} f (t, y (t)) d t

We can conveniently change this into a fixed point iteration.

y_{n + 1} = y_{0} + \int_{x_{0}}^{x} f (t, y_{n} (t)) d t

An Application

Let's apply this to a simple differential equation. I use this same example on the Adomian Decomposition page, but I don't remind repeating it here.

y_{0} = 1

y_{1} = 1 + \int_{0}^{x} d x = 1 + x

y_{2} = 1 + \int_{0}^{x} (1 + x) d x = 1 + x + \frac{x ^{2}}{2 !}

y_{3} = 1 + \int_{0}^{x} (1 + x + \frac{x ^{2}}{2}) d x = 1 + x + \frac{x ^{2}}{2 !} + \frac{x ^{3}}{3 !}

We recognize this as a Taylor series.

y_{n} = k = 0 \sum n \frac{x ^{k}}{k !}

So this is easy enough for a simple linear differential equation, but if we make the linear equation just a little bit more difficult we start to run into issues. It takes many steps in order to get comparable convergence with a first-order euler method. In fact, even for some linear problems it is essentially intractable (symbolically) and we may have to resort to numerical integration.

Symbolic

Using a computer algebra system, one can create an iterative symbolic solution. The trouble with this method is that when nonlinearities are introduced some integrals become intractable. It is possible, however to Taylor expand the nonlinearity and proceed, but the radius of convergence may be finite.

Pade Approximant

Another modification to the picard iteration is you can use a Pade approximant to extrapolate out from a local region. This approximant can then be re-introduced into the iterative scheme using the chebyshev basis.

Chebyshev

Using a Chebyshevv basis it is possible to define an iterative procedure and since all of the integrations are carried out numerically it is simple to implement. Additionally, derivatives may be calculated using matrix operators and their inversions. This is a convenient representation because we have access to both the derivatives and integrals.

Spline

Splines are another possible solution. Instead of using a Chebyshev basis one can use low-order splines. These splines are easy to integrate and also admit simple derivative representations. Additionally, if cubic splines are used then the integration is exact when using Simpson's rule.

Cumulative Summation

One may also resort to a fully numerical solution using Simpson or Trapezoidal methods. The reason for using these methods is that they may be easily transformed into cumulative forms. Methods which use a gaussian quadrature cannot easily calculate cumulative forms, making the anti-derivative difficult to calculate.

Anderson Acceleration

Until now, no explicit use of the residual has been mentioned. Anderson Acceleration uses a regression of the past few iterates to create a linear combination which has the lowest residual error, theoretically accelerating convergence.

Gradient-based Method

Newton-Raphson uses the gradient of the residual to optimize the iterations.