# 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[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? – J. M.'s technical difficulties Jul 27 '15 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? – user1337 Jul 27 '15 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. – John McGee Jul 27 '15 at 13:26
• I meant the recursive formula you've given; something seems missing within the integral… – J. M.'s technical difficulties Jul 27 '15 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/… – Andreas Lauschke Jul 29 '15 at 19:51

I believe that the following code does what you want

For[{n = 1, y[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? – user1337 Jul 27 '15 at 13:29
• This might be useful. – J. M.'s technical difficulties Jul 27 '15 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 – Michael E2 Jul 27 '15 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… – J. M.'s technical difficulties Jul 29 '15 at 1:22
• @J. M. Yes, this one's hard to run in GedankenMathematica. – Jens Jul 29 '15 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. – march Jul 29 '15 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. – J. M.'s technical difficulties Aug 1 '15 at 18:31
• @J M. Happy to hear that you found a good use for this... – Jens Aug 1 '15 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[x]
y[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[x_] = 1
y[n_][x_] := y[n][x] = h[y[n - 1][x]]


Then,

y[x]


outputs

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


and in the process, y[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]
(* {1, 1 + x, 1 + x + x^2 + x^3/3} *)


and in the process, y and y get defined (but not y).

• What is code Picard's interation of y'''[x]=y''[x]+y'[x]+y[x]+x,y=1,y'=2,y''=3.? – Mariusz Iwaniuk Jul 28 '15 at 7:32
• @Mariusz. Why don't you post a separate question with your particular problem, showing us what you tried? – march Jul 28 '15 at 13:58
• I upvoted this, but have also posted an answer that is lightning-fast compared to this one. – Jens Jul 29 '15 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 – Michael E2 Jun 22 '16 at 18:04