Implementing Picard's Iteration for solving ODEs

The answers by march and John McGee become very slow for larger numbers of iteration, to the extent that I had to abort the calculations when going to 7 or 8 iterations.

The reason is that Integrate appears to be trying too many unnecessary simplifications at each level, and these steps proliferate because the integrals are iterated.

The following makes the calculations much faster - by many orders of magnitudes for large iterations. I assume that the initial condition is given at x=0, so the lower integration limit is always 0.

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}
  ]

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 example from the original question:

flow[t_, f_] := f^2

picardSeries[1, flow, 5, x]

(*
==> 1 + x + x^2 + x^3 + x^4 + x^5 + (43 x^6)/45 + (
 13 x^7)/15 + (943 x^8)/1260 + (3497 x^9)/5670 + (
 27523 x^10)/56700 + (1477 x^11)/4050 + (17779 x^12)/68040 + (
 13141 x^13)/73710 + (1019 x^14)/8820 + (63283 x^15)/893025 + (
 43363 x^16)/1058400 + (1080013 x^17)/48580560 + (
 2588 x^18)/229635 + (162179 x^19)/30541455 + (16511 x^20)/7144200 + (
 207509 x^21)/225042300 + (557 x^22)/1666980 + (
 2447 x^23)/22504230 + (16927 x^24)/540101520 + (
 5309 x^25)/675126900 + x^26/595350 + (2 x^27)/6751269 + (
 13 x^28)/315059220 + x^29/236294415 + x^30/3544416225 + \
x^31/109876902975
*)

This result is obtained almost instantaneously. You can also get results for n = 7, 8 with almost no delay (I just don't want to waste space by listing them). These orders are practically out of reach for the other answers.

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

This solution works equally well for vectorial initial-value problems, i.e., the flow can be a vector function and the initial condition a vector. This is needed, e.g., if you want to apply this method to a higher-order differential equation for a scalar function by converting it to a first-order equation for a vector function (a standard technique I don't think I have to go into in detail).

Application: a two-state system

As a non-trivial example of a vectorial initial-value problem, here is the solution to a quantum time evolution of a two-state system that is initially in the state {1,0} and subjected to a periodic driving field:

Clear[e0, v0, ω, t];

$Assumptions = {e0 > 0, v0 > 0, ω > 0, t > 0};

flow1[t_, ψ_] := -I {{0, E^(-I t (e0 - ω)) v0}, 
                     {E^(I t (e0 - ω)) v0, 0}}.ψ

FullSimplify[picardSeries[{1, 0}, flow1, 3, t]]

$$\left( \begin{array}{c} 1+\frac{\text{v0}^2 \left(i t (\text{e0}-\omega )+e^{-i t (\text{e0}-\omega )}-1\right)}{(\text{e0}-\omega )^2} \\ \frac{\text{v0} \left(e^{i t (\text{e0}-\omega )} \left(2 \text{v0}^2-(\text{e0}-\omega ) \left(\text{e0}+i t \text{v0}^2-\omega \right)\right)+(\text{e0}-\omega ) \left(\text{e0}-i t \text{v0}^2-\omega \right)-2 \text{v0}^2\right)}{(\text{e0}-\omega )^3} \\ \end{array} \right)$$

This illustrates the important lesson that Simplify (let alone FullSimplify) shouldn't be called (or implicitly invoked) in the iteration step. On the other hand, using FullSimplify on the final result (which here contains many symbolic constants, too), is OK.

Edit: speed considerations

The individually labelled integration variables in each step create a slight overhead (linear in n, i.e., rather benign).

But I found that nested integrations with identically labelled integration variables run more slowly (although they produce correct results). This is very noticeable for the second example, flow1, going e.g. to n = 9. Therefore, the small overhead in creating unique variable names is justified when going to large n.

Here are some possibilities more in line with Mathematica idioms (i.e. that avoid using procedural loops).

Option 1

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

Option 2

To set the ys in the process:

Clear[y, h]
h = Function[{t}, 1 + Integrate[#^2, {x, 0, t}]][x] & ;
NestList[{1 + First@#, y[1 + First@#][x_] = h@Last@#} &, {0, 1}, 2]
y[2][x]
y[1][x]

results in

(* {{0, 1}, {1, 1 + x}, {2, 1 + x + x^2 + x^3/3}} *)
(* 1 + x + x^2 + x^3/3 *)
(* 1 + x *)

Option 3

Or we could memo-ize while using a recursive definition of y:

Clear[y, h]
h = Function[{t}, 1 + Integrate[#^2, {x, 0, t}]][x] & ;
y[0][x_] = 1
y[n_][x_] := y[n][x] = h[y[n - 1][x]]

Then,

y[2][x]

outputs

(* 1 + x + x^2 + x^3/3 *)

and in the process, y[1][x] also gets defined.

Option 4

Using FoldList:

h = (y[#2][x_] = Function[{t}, 1 + Integrate[#1^2, {x, 0, t}]][x]) &;
FoldList[h, 1, Range[2]]
(* {1, 1 + x, 1 + x + x^2 + x^3/3} *)

and in the process, y[1] and y[2] get defined (but not y[0]).

I believe that the following code does what you want

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