# Implementing Picard's Iteration for solving ODEs

Picard's Iteration is a way of solving the IVP

$$y'(x)=f(x,y(x)), \quad y(x_0)=y_0$$

It consists of defining the following sequence of functions recursively:

$$y_0(x):=y_0 \\ y_{n}(x):=y_0+\int_{x_0}^x f(t,y_{n-1}(t)) \mathrm dt.$$

I've tried implementing it in Mathematica, for the particular problem

$$y'=y^2,y(0)=1$$

as follows

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


which doesn't work at all. What am I doing wrong here, and how can I fix it?

Thank you!

• Your formula for Picard-Lindelöf does not look right. Are you sure you didn't forget, say, a subscript? Jul 27, 2015 at 13:22
• @J. M. instead of subscripting my functions as $\{y_n\}$ I've decided to use an array of functions $\{y[n] \}$. Is that the reason why it fails? Jul 27, 2015 at 13:23
• I don't think you want to use := (delayed assignment) in this case. Change your assignments of y[n] to = and start at n=1. Jul 27, 2015 at 13:26
• I meant the recursive formula you've given; something seems missing within the integral… Jul 27, 2015 at 13:26
• FWIW, there was an article about Picard iteration in TMJ about ten years ago, I've found it online: mathematica-journal.com/issue/v10i1/contents/Picard/… Jul 29, 2015 at 19:51

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

• This seems to work beautifully, thanks! I'm used to adding the : for unknown reasons. Could you please refer me to a place where I could learn what it means? Jul 27, 2015 at 13:29
• This might be useful. Jul 27, 2015 at 13:34
• @user1337 In the documentation center, see tutorial/ImmediateAndDelayedDefinitions, and on the site, look at mathematica.stackexchange.com/search?q=Set+SetDelayed+is%3Aq Jul 27, 2015 at 16:48

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.

• Okay, this code got me excited. I will have to try this the next time I see a computer… Jul 29, 2015 at 1:22
• @J. M. Yes, this one's hard to run in GedankenMathematica.
– Jens
Jul 29, 2015 at 1:25
• This is excellent. I always forget how Integrate tries to be clever (and I did wonder why mine would slow down so much once n is 8 or so). I guess my solutions are more syntactic fluff and less useful for the actual problem. Jul 29, 2015 at 2:26
• Just tried this out. I love it! This is more efficient than the method I had for generating series coefficients of the first Painlevé transcendent. Aug 1, 2015 at 18:31
• @J M. Happy to hear that you found a good use for this...
– Jens
Aug 1, 2015 at 19:16

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]).

• What is code Picard's interation of y'''[x]=y''[x]+y'[x]+y[x]+x,y[0]=1,y'[0]=2,y''[0]=3.? Jul 28, 2015 at 7:32
• @Mariusz. Why don't you post a separate question with your particular problem, showing us what you tried? Jul 28, 2015 at 13:58
• I upvoted this, but have also posted an answer that is lightning-fast compared to this one.
– Jens
Jul 29, 2015 at 0:07

for example for function y'=f(x,y)=x-y , I think the easiest way is...

f[x_, y_] := x - y

yn[x_,n_,x0_,y0_]:=If[n>0, y0+Integrate[f[s,yn[s,n-1,x0,y0]],{s,x0,x}], y0]

example (y(0)=1 and 3 iterations):

yn[x, 3, 0, 1]

1 - x + x^2 - x^3/3 + x^4/24

• You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this meta Q&A helpful Jun 22, 2016 at 18:04